Sequence retrieval and variant calling

Obtain accession numbers for 2773 clinicalstrains of Mtb from 10 previous studies(Supplementary file 1), representing all 4 lineagesfrom 9 countries.
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
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
Download the trimommatic software using the following link http://www.usadellab.org/cms/index.php?page=trimmomatic.
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
Download the H37Rv reference genomeand corresponding gff file from NCBI database(https://www.ncbi.nlm.nih.gov/datasets/genome/) with accession number ASM19595.
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
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
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
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
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
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
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