workflow_eb.sh

NUMPROC=30
MEMORY=225
PROJ=${PWD##*/}

#Trim galore: 
#trim adapter seqs, overrep seqs, qual < 20 (also dedup w/ erika's suggestion?)
if [ ! -d "trimming" ]; then mkdir trimming; fi
if [ ! -d "assembly" ]; then mkdir assembly; fi
if [ ! -d "quast" ]; then mkdir quast; fi
if [ ! -d "mapping" ]; then mkdir mapping; fi
if [ ! -d "variantcalling" ]; then mkdir variantcalling; fi
      

find  ./raw  -name "*_R1*.f*.gz" | cut -d "_" -f1 | parallel --dryrun -j $NUMPROC trim_galore --illumina --fastqc -o trimming/ {}\_R1.fastq.gz {}\_R2.fastq.gz 

find ./raw -name "*R1*.gz"* | cut -d "_" -f1 | cut -d "/" -f3 | sort | uniq > popslist
 

if test -n "$(find ./raw -maxdepth 1 -name '*R2*' -print -quit)"; then
    echo reverse reads found. proceeding with paired end assembly
    #setspades.R
    echo "library(here)" > spades_yaml.R
    echo "workingdir <- paste0(here(),\"/\")" >> spades_yaml.R
    echo "popslist <- read.table(paste0(workingdir,\"/popslist\"))" >> spades_yaml.R
    echo "setwd(workingdir)" >> spades_yaml.R
    echo "write(paste0('[ \n',  '   { \n', '     orientation: \"fr\", \n', '     type: \"paired-end\", \n', '     right reads: ['), file = \"libraries.yaml\", append = F)" >> spades_yaml.R
    echo "for(i in 1:nrow(popslist)){write(paste('       \"', workingdir, \"trimming/\", popslist[i,1], '_R1_trimmed.fq.gz\",', sep = ''), file = \"libraries.yaml\", append = T)}" >> spades_yaml.R
    echo "write(paste0('       ], \n', '       left reads: ['), file = \"libraries.yaml\", append = T)" >> spades_yaml.R
    echo "for(i in 1:nrow(popslist)){write(paste0('       \"', workingdir,\"trimming/\", popslist[i,1], '_R2_trimmed.fq.gz\",'), file = \"libraries.yaml\", append = T)}" >> spades_yaml.R                   
    echo "write(paste0('       ] \n', '   } \n', ']'), file = \"libraries.yaml\", append = T)" >> spades_yaml.R


    chmod 755 spades_yaml.R
    Rscript spades_yaml.R

    #paired-end assembly
    spades.py --dataset libraries.yaml -k 71,81,91,99,121,127 --careful -t 30 -m 300 -o ./assembly/$PROJ_usr_kmer
    spades.py --dataset libraries.yaml --careful -t $NUMPROC -m $MEMORY -o ./assembly/$PROJ_default_kmer
else
    echo no reverse reads found. proceeding with single end assebmly
    #single end assembly
    pools=(` find ./trimming -name "*R1*.gz" | sort | uniq | cat `) 
    libnum=( `seq 1 "${#pools[@]}"` )
    libnum=("${libnum[@]/#/--s}")
    unset reads
    for (( i=0; i<${#libnum[*]}; ++i)); do reads+=( ${libnum[$i]} ${pools[$i]} ); done
    #kmers="$(seq -s ',' 21 2 127)"
    #spades.py "${reads[@]}" -k "$kmers" --careful -t $NUMPROC -m $MEMORY -o ./assembly/ipyRAD_all_kmer_output

    spades.py "${reads[@]}" -k 71,81,91,121,127 --careful -t $NUMPROC -m $MEMORY -o ./assembly/$PROJ_usr_kmer
    spades.py "${reads[@]}" --careful -t $NUMPROC -m $MEMORY -o ./assembly/$PROJ_default_kmer
fi

#####QUAST#####
#per KMER
for genome in `find ./assembly -name "*usr*" -type d | cut -d "/" -f3`; do
    KMERS=(` find ./assembly/$genome/ -maxdepth 1 -name "K*" | cut -d "/" -f4 | sort | uniq `)
    cp ./assembly/${genome}/contigs.fasta ./quast/"${genome}_contigs.fasta"
         for KMER in `find ./assembly/$genome/ -maxdepth 1 -name "K*" | cut -d "/" -f4 | sort | uniq`; do
             cp ./assembly/${genome}/${KMER}/final_contigs.fasta ./quast/${genome}_${KMER}.fasta
        done
    by_kmer=("${KMERS[@]/#/./quast/${genome}_}")
    by_kmer_files=("${by_kmer[@]/%/.fasta}")
    quast.py -e -t $NUMPROC -m $MEMORY ./quast/${genome}_contigs.fasta "${by_kmer_files[@]}" -o ./quast/${genome}_results/
done

#all default KMER selection
all_user=(` find ./quast/*usr*contigs.fasta  -type f `)
quast.py -e -t $NUMPROC -m $MEMORY  "${all_user[@]}" -o ./quast/usr_spades/

#spades default settings
#all default KMER selection
all_default=(` find ./quast/*default*contigs.fasta  -type f `)
quast.py -e -t $NUMPROC -m $MEMORY "${all_default[@]}" -o ./quast/default_spades/ 

#########STOP######
#before proceeding: evaluate your quast results. Choose desired de novo, copy to ./assemly with a name to include chosen.fasta

#THECHOSEN=./assembly/<nameofchosenreference>
THECHOSEN="./assembly/*chosen*"

LIBS=(` find  ./trimming  -name "*trimmed.fq.gz" | cut -d "/" -f3 `)
POPS=(` find  ./trimming  -name "*trimmed.fq.gz" | cut -d "_" -f1 | cut -d "/" -f3 | sort | uniq `)
optA=1
optB=4
optO=6
SPLITS=2

cp $THECHOSEN ./mapping/reference.fasta
cp ./trimming/*_trimmed.fq* ./mapping/
for i in "${POPS[@]}"; do mv ./mapping/${i}_R1__trimmed.fq.gz ./mapping/${i}.fq.gz; done

find ./mapping -name "*.fq.gz" -o -name "*.fasta" | parallel -j $NUMPROC bwa index {}

#PAIRED
for i in "${POPS[@]}"; do
    echo $i | parallel -j 1 bwa mem reference.fasta  {}.F.fq.gz {}.R.fq.gz -L 20,5 -t $NUMPROC -M -T 35 -A $optA -B $optB -O $optO -R \"@RG\\tID:{}\\tSM:{}\\tPL:Illumina\" 2> bwa.$i.log | mawk '$6 !~/[2-9].[SH]/ && $6 !~ /[1-9][0-9].[SH]/' | samtools view -@ $NUMPROC -q 1 -SbT reference.fasta - > $i.bam 2>$i.bam.log
    samtools sort -@ $NUMPROC -O BAM -o $i.sort.bam $i.bam
done

#Single end
for i in "${POPS[@]}"; do
    echo $i | parallel -j $SPLITS bwa mem ./mapping/reference.fasta ./mapping/{}.fq.gz -L 20,5 -t $(( $NUMPROC / $SPLITS )) -M -T 35 -A $optA -B $optB -O $optO -R \"@RG\\tID:{}\\tSM:{}\\tPL:Illumina\" 2> ./mapping/bwa.$i.log | mawk '$6 !~/[2-9].[SH]/ && $6 !~ /[1-9][0-9].[SH]/' | samtools view -@ $NUMPROC -q 1 -SbT ./mapping/reference.fasta - > ./mapping/$i.bam 2>./mapping/$i.bam.log
    samtools sort -@ $NUMPROC -O BAM -o ./mapping/$i.sort.bam ./mapping/$i.bam
done

find ./mapping -name "*.sort.bam" | parallel -j $NUMPROC samtools index {}
BAMS=("${POPS[@]/#/./mapping/}")
printf "%s\n" "${BAMS[@]/%/.sort.bam}" > ./variantcalling/bams.list
samtools merge -@ $$NUMPROC -b ./variantcalling/bams.list -f ./variantcalling/merged_pools.bam
samtools index ./variantcalling/merged_pools.bam

bamToBed -i ./variantcalling/merged_pools.bam | bedtools merge -i - > ./variantcalling/mapped.bed
bedtools coverage -b ./variantcalling/merged_pools.bam -a ./variantcalling/mapped.bed -counts -sorted > ./variantcalling/cov.stats

DP=$(mawk '{print $4}' ./variantcalling/cov.stats | sort -rn | perl -e '$d=.001;@l=<>;print $l[int($d*@l)]')
CC=$( mawk -v x=$DP '$4 < x' ./variantcalling/cov.stats | mawk '{len=$3-$2;lc=len*$4;tl=tl+lc} END {OFMT = "%.0f";print tl/"'$NUMPROC'"}')

mawk -v x=$DP '$4 < x' ./variantcalling/cov.stats | sort -V -k1,1 -k2,2 | mawk -v cutoff=$CC 'BEGIN{i=1} 
{
len=$3-$2;lc=len*$4;cov = cov + lc
if ( cov < cutoff) {x="./variantcalling/mapped."i".bed";print $1"\t"$2"\t"$3 > x}
else {i=i+1; x="./variantcalling/mapped."i".bed"; print $1"\t"$2"\t"$3 > x; cov=0}
}' 

mawk -v OFS='\t' {'print $1,$2'} ./mapping/reference.fasta.fai > ./variantcalling/genome.file
for i in "${POPS[@]}"; do
    echo $i | parallel -j $NUMPROC "bedtools coverage -b ./mapping/{}.sort.bam -a ./variantcalling/mapped.bed -counts -sorted -g ./variantcalling/genome.file > ./variantcalling/{}.cov.stats"
done
cat ./variantcalling/*.cov.stats | sort -k1,1 -k2,2n | bedtools merge -i - -c 4 -o sum > ./variantcalling/cov.stats

DP=$(mawk '{print $4}' ./variantcalling/cov.stats | sort -rn | perl -e '$d=.001;@l=<>;print $l[int($d*@l)]')
CC=$( mawk -v x=$DP '$4 < x' ./variantcalling/cov.stats | mawk '{len=$3-$2;lc=len*$4;tl=tl+lc} END {OFMT = "%.0f";print tl/"'$(( $NUMPROC * 2 ))'"}')

find ./trimming/ -name "*.fq.gz" | cat | parallel -j $NUMPROC "gunzip -c {} | head -2 | tail -1 >> ./variantcalling/lengths.txt"
MaxLen=$(mawk '{ print length() | "sort -rn" }' ./variantcalling/lengths.txt | head -1)
MaxLen2=$(( $MaxLen / 2 ))
ML1=$(( $MaxLen2 + 1 ))
MaxCov=$(mawk '{print $4}' ./variantcalling/cov.stats | sort -rn | head -1)
MaxCov2=$(( $MaxCov * 0.8 ))

mawk '$4 > $MaxCov2' ./variantcalling/cov.stats > ./variantcalling/cov.high.stats
mawk '$4 <= $MaxCov2' ./variantcalling/cov.stats > ./variantcalling/cov.low.stats

#split high coverage intervals into smaller, 1/2 read-length sized intervals

bedtools makewindows -b ./variantcalling/cov.high.stats -w $MaxLen2 -s $ML1 > ./variantcalling/temp1.bed
bedtools intersect -a ./variantcalling/cov.stats -b ./variantcalling/temp1.bed > ./variantcalling/temp.cov.stats

cat ./variantcalling/temp.cov.stats ./variantcalling/cov.low.stats > ./variantcalling/cov.split.stats
#rm ./variantcalling/cov.high.stats ./variantcalling/cov.low.stats ./variantcalling/temp1.bed ./variantcalling/temp.cov.stats

TT=$(( $MaxLen2 * 1000000 ))

mawk -v x=$DP '$4 < x' ./variantcalling/cov.split.stats |sort -V -k1,1 -k2,2 | mawk -v cutoff=$CC -v tt=$TT 'BEGIN{i=1}
{len=$3-$2;lc=len*$4;cov = cov + lc
if (NR == 1 && lc > tt) {x="./variantcalling/mapped."i".bed";print $1"\t"$2"\t"$3 > x; i=i+1; e=1}
else if ( cov < cutoff && lc < tt) {x="./variantcalling/mapped."i".bed";print $1"\t"$2"\t"$3 > x; e=0}
else if (lc > tt && e > 0 ) {x="./variantcalling/mapped."i".bed"; print $1"\t"$2"\t"$3 > x; cov=0;i=i+1; e=1}
else if (lc > tt && e < 1 ) {i=i+1; x="./variantcalling/mapped."i".bed"; print $1"\t"$2"\t"$3 > x; cov=0;i=i+1;e=1}
else if (cov > cutoff && lc < tt ) {i=i+1; x="./variantcalling/mapped."i".bed"; print $1"\t"$2"\t"$3 > x; cov=lc;e=0}
}'

paste popslist popslist > popmap

    call_genos(){
    samtools view -@ $NUMPROC -b -1 -L ./variantcalling/mapped.$1.bed -o ./variantcalling/split.$1.bam ./variantcalling/merged_pools.bam
    samtools index ./variantcalling/split.$1.bam
    freebayes -b ./variantcalling/split.$1.bam -t ./variantcalling/mapped.$1.bed -v ./variantcalling/raw.$1.vcf -f ./mapping/reference.fasta -E 3 -F 0.05 -C 1 --min-coverage 10 --min-repeat-entropy 1 --populations popmap -V --pooled-continuous
    #vcflib vcffilter -f "QUAL > 10" raw.$1.vcf > qualfil.$1.vcf
    rm ./variantcalling/split.$1.bam*
    }
    export -f call_genos

    ulimit -s unlimited
    ls ./variantcalling/mapped.*.bed | sed 's/mapped.//g' | sed 's/.bed//g' | cut -d "/" -f3 | shuf | parallel --env call_genos -j $NUMPROC --no-notice call_genos {}
    #rm ./variantcalling/mapped.*.bed 

    rename -n -e 's/\d+/sprintf("%02d",$&)/e' -- ./variantcalling/*.vcf
    vcflib vcfcombine ./variantcalling/raw.*.vcf > ./variantcalling/TotalRawSNPs.vcf
    #
    ls

if [ ! -d "./variantcalling/raw.vcf" ]; then
     mkdir ./variantcalling/raw.vcf
fi

    mv ./variantcalling/raw.*.vcf ./variantcalling/raw.vcf
    ls
#if output is polyploid for non GT calls (ie ././.)
grep -v "^#" ./variantcalling/TotalRawSNPs.vcf | cut -f 10- | cut -d ':' -f 1 | sort | uniq
sed 's/\.\/\.\/\./\.\/\./' ./variantcalling/TotalRawSNPs.vcf > ./variantcalling/TotalRawSNPs_fix.vcf
rm ./variantcalling/TotalRawSNPs.vcf; mv ./variantcalling/TotalRawSNPs_fix.vcf ./variantcalling/TotalRawSNPs.vcf