philippbayer

First, to download EDTA:

module load singularity
singularity pull EDTA.sif docker://quay.io/biocontainers/edta:1.9.4--0

That'll make a new file called EDTA.sif containing everything in the EDTA v1.9.4 container.

Then we have a problem: Pawsey allows only 1 million files per user and running several EDTA runs for several genomes at once will hit that limit.

philippbayer

```{r setup}
library(tidyverse)
library(ggrepel)
```

```{r}
df <- readxl::read_xlsx('./Covid_vs_State.xlsx')
head(df)
```

philippbayer

module load singularity

# we have to transform the docker image into a singularity thingy. Be careful to use the GPU docker image 0.9.0-gpu, not 0.9.0
# the following command doesn’t work on copyQ?? I get `mksquashfs not found`, works fine on other systems

singularity pull $MYGROUP/singularity/myRepository/deepvariant.sif docker://google/deepvariant:0.9.0-gpu
# you may have to make the above folders

Then, using the example data described here: https://github.com/google/deepvariant/blob/r0.9/docs/deepvariant-quick-start.md

philippbayer

I had a ddRAD project with two restriction enzymes not supported by Stacks - Hinfl and HpyCH4IV. This is how I added them to Stacks:

1. Download the code, open src/renz.cc
2. Check NEB for the restriction enzyme sequences: https://www.nebiolabs.com.au/products/r0155-hinfi#Product%20Information and https://www.nebiolabs.com.au/products/r0619-hpych4iv#Product%20Information
3. Understand how Stacks encode restriction enzymes - Hinfl is G CUT ANTC, which in the comments of renz.cc is encoded as G/ANTC
4. Choose the longer piece (in this case, ANTC) and translate that into all four possibilities AGTC AATC ATTC ACTC, and the reverse complement GACT GATT GAAT GAGT. For others like A/CGT it's easier, take the longer piece and its reverse-complement, so it's just CGT and ACG
5. Add that to the code, and don't forget your semicolons (I think the first letter of the RE name has to be lower case??)

philippbayer

import subprocess
from Bio.Seq import Seq
from Bio import SeqIO

trimal = '/path/to/trimal'
muscle = '/path/to/muscle3.8.31_i86linux64'
# this is the fasta containing the single reference gene we use for trimming all of our seqeunces 
ref_gene = '/path/to/single/ref.fasta'

all_out = open('file_to_fix.Trimmed.fasta', 'w')

philippbayer

library(tidyverse)
library(lubridate)
library(viridis)
library(ggridges)
# fix typos

a <- read_csv('./biiiiirdsssss.csv')
a$date <- gsub("16/01/017",  "16/01/2017", a$date)
a$date <- gsub("26/07/2026",  "26/07/2016", a$date)
a$date <- dmy(a$date)

philippbayer

# assuming that all results tables of BLINK are in the current folder, and assuming that all results tables were gzipped
# usage: python makeSummaryTable.py > Summary_Table.tsv

# this script will make one large summary table in the format rMVP requires, SNP, chrom, position, and then x columns for x phenotypes - each cell one p-value
# use grep -v to kick out regions of the genome you don't want to plot (unplaced contigs etc.)

import glob
import gzip
from collections import OrderedDict
all_phenos = ['SNP','chr','pos']

philippbayer

library(tidyverse)

library(lubridate) # for the days() function
library(fuzzyjoin) # for the fuzzy_left_join() function

climate <- read_tsv('./timetable.csv')
#  date         temperature
#1 2018-11-21   100 
#2 2018-11-22    80 

philippbayer

'''

This will print something like 
File	Subset	TN	TCNL	NL	CN	N	TCN	TNL	CNL
assembly.fa_chopped_out.xml_txt	complete	0	0	0	0	0	0	29	17
assembly.fa_chopped_out.xml_txt	complete (pseudogene)	0	0	0	0	0	0	16	7
assembly.fa_chopped_out.xml_txt	partial	3	0	1	0	4	1	0	0
assembly.fa_chopped_out.xml_txt	partial (pseudogene)	1	0	0	5	1	2	2	2

philippbayer