#!/bin/bash

# Author: Richard L Tillett

# open-source tool alignment, variant, and consensus methods

# Alignment performed using bowtie2 with shown settings
# Our exact raw FASTQ files could not be included in this archive
# due to the likely presence of human RNA and the consequent privacy
# concerns.
# Human-depleted FASTQs can be obtained from NCBI instead.

## alignment to reference sequence
# bowtie2 --local --no-unal \
# -x wuhan.fasta \
# -1 sampleA_R1.fastq.gz \
# -2 sampleA_R2.fastq.gz \
# -S alignedA.sam
# gzip alignedA.sam
# samtools sort -o alignedA.bam alignedA.sam.gz
# samtools index alignedA.bam

# mark dupuplicates

picard MarkDuplicates \
      I=alignedA.bam \
      O=dedupedA.bam \
      M=dedupedA.log \
      REMOVE_DUPLICATES=1

# (repeat above steps for sample B.)

# use freebayes with the bamaddrg feeding it the samples

bamaddrg \
-b dedupedA.bam -s CaseA -r CaseA \
-b dedupedB.bam -s CaseB -r CaseB | \
freebayes -p 1 \
-f wuhan.fasta \
-F 0.7 -C 5 -c > cases_only.vcf

# split vcf files per sample with vcflib tools

vcfkeepsamples cases_only.vcf CaseA | vcffixup - | vcffilter -f "AC > 0" > caseA.vcf
vcfkeepsamples cases_only.vcf CaseB | vcffixup - | vcffilter -f "AC > 0" > caseB.vcf

# generate coverage maps with bbmap's pileup.sh

pileup.sh \
in=dedupedA.bam \
ref=wuhan.fasta \
out=caseA.cov \
basecov=caseA.basecov 

# repeat above for case B

# Masking low-coverage regions
#
# replace sequences on the reference with the letter N where coverage
# is below threshold in the sample with bbmap's applyvariants.sh
# note: this step ought to be able to write the variants as well, but
# this functionality is buggy in version 38.86 of bbmap. The application 
# of unknown Ns does still work as intended, though. To achieve this
# we supply `applyvariants.sh` with a blank empty vcf file.

touch empty.vcf

applyvariants.sh \
in=wuhan.fasta \
vcf=empty.vcf \
basecov=caseA.basecov \
out=caseAs_NNN.fasta \
mincov=5 \
-da

# repeat above for case B

# Substitute variant alleles from each sample's vcf
#
# finally, apply the variants using bcftools consensus.
# There are two preparatory steps that precede this

bgzip caseA.vcf
tabix -p vcf caseA.vcf.gz

bcftools consensus -f caseAs_NNN.fasta \
-p CaseA -o caseA_sequence.fasta \
caseA.vcf.gz

# repeat above for case B
