Sequence retrieval and variant calling

Sequence retrieval and variant calling

1. Obtain accession numbers for 2773 clinicalstrains of Mtb from 10 previous studies(Supplementary file 1), representing all 4 lineagesfrom 9 countries.
    
2. All the sequenced data of clinical strains is given at the EBA and NCBI. Use SRA toolkit (NCBI) to retrieve fastq data from NCBI databases (https://www.ncbi.nlm.nih.gov/). Sort sequencesbased on the paired end and singleend fastq files.
   
   mkdir -p /path/to/dump/folder 
   module load sratoolkit
   fasterq-dump SRR_ accession number -O /path/to/dump/folder
    
3. For downloading files from EBA (https://www.ebi.ac.uk/), use ena_file_dowloader.jar program. All accessions should be stored in a file (accession.txt) and provide file to the downloader as follows:
   
   java -jar ena-file-downloader.jar accession.txt
    
4. Download       the        trimommatic       software       using        the        following       link http://www.usadellab.org/cms/index.php?page=trimmomatic.
    
5. Quality filter FastQC files using Trimmomatic software (Bolger et al., 2014). Remove adapter and E. coli sequence contaminations, analyze the sequences at a sliding window of 4 bp, and those with an average phred value of 15 were clipped. Set parameters for the trimmomatic- based QC as (Leading:3, slidingwindow:4:15 trailing:3 minlen:60). Remove reads with a filtered lengthof <60 bp. 
   
   For fastqc based QC:
   fastqc -f fastq /path/to/Rawsample/*_1.fastq.gz
   fastqc -f fastq /path/to/Rawsample/*_2.fastq.gz
   
   For Trimmomatic based QC:
   java -jar trimmomatic-0.36.jar PE -phred33 -threads 40 file-R1 file_R2 file_R1_PE.fastq.gz file_R1_SE.fastq.gz file_R2_PE.fastq.gz file_R2_SE.fastq.gz ILLUMINACLIP:TruSeq3-PE- 2.fa:2:30:10 LEADING:3SLIDINGWINDOW:4:15 TRAILING:3 MINLEN:60
    
6. Download the H37Rv reference genomeand corresponding gff file from NCBI database(https://www.ncbi.nlm.nih.gov/datasets/genome/) with accession number ASM19595.
    
7. Create indexing file for MTB genome. Perform this using the Burrows Wheeler Alignment (BWA) mem algorithm (Li and Durbin,2009). A detailedmanual can be checked at the websitehttps://bio-bwa.sourceforge.net/bwa.shtml.
   
   mkdir /path/to/reference_directory 
   cd /path/to/reference_directory
   cp path/to/downloaded/reference_file 
   bwa index REF.fa
    
8. Map the filteredSE/PE reads on the H37Rv reference genomeusing the BurrowsWheeler Alignment (BWA) mem algorithm (Li and Durbin,2009).
   
   bwa mem -M -t 40 -R "@RG\tID:FLOWCELL1.LANE1\tPL:ILLUMINA\tLB:LIB-$BASE\tSM:$BASE" REF.fa file_R1_PE.fastq.gz file_R2_PE.fastq.gz > file.samlignment
    
9. Convert the generated SAM file into .bam file using samtools(Li et al, 2009).
   
   Samtools view -Sb /path/to/Alignment.sam > /path/to/Alignment.bam Samtools sort /path/to/Alignment.bam /path/to/Alignment.srt.bam
    
10. Determine the mappingquality of the sorted sam file using qualimap software.
    Mkdir /path/to/quality_check/file
    qualimap bamqc -bam file.srt.bam --java-mem-size=15G -outdir file.qualimap.srtbam -outfile file.qualimap_srt.bam.pdf
     
11. Quality check improvement and realignments were done using GATK pipeline (Van der Auwera et al., 2013). A detailed manual of the GATK software is already published at the websitehttps://gatk.broadinstitute.org/hc/en-us. Variants (SNP and small InDels)were predicted in a batch mode with Platypus or GATK pipeline (Rimmer et al., 2014) using the followingcommand.
    
    mkdir /path/to/combined_Annotation
    gatk.jar GenotypeGVCFs -R REF --dbsnp known_file.vcf --variant file.recal.g.vcf --output file.combined_annotated.vcf
     
12. Annotation of the identified SNPs using SnpEff toolbox (Cingolani et al., 2012). [http://pcingola.github.io/SnpEff/]
    
    Java -jar snpEff.jar mycobacterium file.combined.filtered.vcf >/path/to/combined_Annotation/file.combined_annotated.vcf
    /path/to/Tab_annotated/
    
    cat /path/to/combined_Annotation/file.combined_annotated.vcf | ./vcfEffOnePerLine.pl | java-Xmx40g -jar SnpSift.jar extractFields -s \",\" -e \".\" - CHROM POS REF ALT\"ANN[*].GENE\" \"ANN [*].EFFECT\" \"ANN[*].IMPACT\" \"ANN[*].BIOTYPE\"\"ANN[*].HGVS_P\" | grep -E -v \’(^#| downstream_gene_variant|upstream_gene_variant)\’ >/path/to/Tab_annotated/ file.combined_annotated_tabdelim.txt
     
13. Finally, ~160,000 SNPs from 2773 accessions with a variant dataset with an MAF of >10% selectedfor the final analysis. The. vcf file converted to. hmp file with Tassel (Bradbury et al., 2007). SNPs with read depth <5 or mapping quality of <30 marked as missing. SNPs with missing calls in >40% of the accessions removed from the analysis.
    
    vcftools --vcf file.combined_annotated.vcf --max-missing 0.2 --mac 10 --minQ 30 --recode --
    recode-INFO-all --out file.combined_mac10_max_missing_0.2.hmp