Finding SNP Sites in the Whole Genome of a Plant Using GATK

-
As I have discussed in the previous blog, we already have aligned SAM (.sam) file to be further utilized for analyzing genomic variations. In this post, I am going to instruct the steps involved in the analysis of genomic variations among the closely related organisms. In particular, I am presenting the workflow pipeline for the GATK (Genome Analysis Tool Kit) HaplotypeCaller for detecting SNP (Single Nucleotide Polymorphism) sites and their annotation by snpEff. 

Note: In order to carry out SNP sites detection by GATK, there should be SAM file of the whole genome sequence of an organism, preinstalled GATK software version 4.1.0.1 and above (McKenna et al., 2010) with its accessory packages, Samtools, and snpEff  4.3t (Cingolani et al., 2012) in Linux server.

Image 1. a screenshot of a home page of gatk website (source: https://software.broadinstitute.org/gatk/)

In brief
1. get the SAM file of the whole genome sequence of the sample organism;
2. do SNP call with HaplotypeCaller in GATK4.1.0.1;
3. perform joint genotyping of HaplotypeCaller operated files (for multiple numbers of whole genome sequences);
4. filter the variants after they have been selected;
5. annotate the SNP and INDEL sites with snpEff.

Image 2. a screenshot of a page of snpEff on sourceforge web repository (source: http://snpeff.sourceforge.net/)

In detail
1. Create sorted BAM file from aligned SAM file
In this step, we are going to create a sorted BAM file, which will have a sequence sorted according to coordinates of the reference genome, using the following command. 

java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar SortSam \ 
    --INPUT=VRiparia.sam \ 
    --OUTPUT=VRiparia_sorted.bam \ 
    --SORT_ORDER=coordinate
--TMP_DIR=/B42T/mukesh/tmp/

Here, VRiparia.sam is an aligned SAM file of Vitis riparia whole genome sequence (the sample sequence) and TMP_DIR is a temporary directory to manage memory problems in Linux.

2. Create markduplicated BAM file or to mark duplicate reads in sorted BAM file
java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar MarkDuplicates \ 
    --INPUT=VRiparia_sorted.bam \ 
    --OUTPUT=VRiparia_dedup_reads.bam \
    --METRICS_FILE=VRiparia_dedup_metrics.txt
--TMP_DIR=/B42T/mukesh/tmp/

OR

java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar MarkDuplicates \
--INPUT=VRiparia_sorted.bam\                                  
--OUTPUT=VRiparia_markduplicates.bam \
--METRICS_FILE=VRiparia_markduplicates_metrics.txt \ 
--OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 \              
--CREATE_INDEX=true \                                       
--TMP_DIR=/tmp

3. Create index of reference genome
samtools faidx Vitis_vinifera.fa

4. Create a sequence dictionary for a reference genome 
java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar CreateSequenceDictionary \       
    --REFERENCE=/B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \ 
    --OUTPUT=Vitis_vinifera.fa.dict
5. Add and replace read groups in markduplicated .bam file
java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar AddOrReplaceReadGroups \ 
    --INPUT=/B42T/mukesh/DNA_Seq_Work/VRiparia_dedup_reads.bam \ 
    --OUTPUT=/B42T/mukesh/DNA_Seq_Work/VRiparia_addRG.bam \ 
    --RGID=HWI-D00525a \                                       
    --RGLB= VRiparia_11 \ 
    --RGPL=illumina \ 
    --RGPU=C6DFYANXX \ 
    --RGSM=SRP161488 \ 

6. Create bam index of read group added BAM file
java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar BuildBamIndex \ 
    --INPUT=VRiparia_addRG.bam 

7. to index the .vcf file with --known-sites of SNPs [when no --known-sites, may be skipped]
java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar IndexFeatureFile \
-F /B42T/mukesh/DNA_Seq_Work/VvSNPs_EnsemblPlants/vitis_vinifera.vcf

8. Do base quality score recalibration [when no --known-sites, may be skipped] 
 java -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar BaseRecalibrator \
   -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \
   --input /B42T/mukesh/DNA_Seq_Work/VRiparia__addRG.bam \
   --known-sites /B42T/mukesh/DNA_Seq_Work/VvSNPs_EnsemblPlants/vitis_vinifera.vcf \
   --output VRiparia_recal_data.table

9. Run HaplotypeCaller for Short Variants Discovery   
 java -Xmx64G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar HaplotypeCaller -nct 32 -I \
   -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \ 
   --input /B42T/mukesh/FastqFiles/VRiparia_addRG.bam\
   --emit-ref-confidence GVCF
   --output VRiparia_raw_SNPs_INDELS.g.vcf 

After we get the GVCF (or .g.vcf) file with short variants like SNPs and INDELs, we can follow two ways to perform joint genotyping:
a) by giving commands in separate steps for combining g.vcf files and GenotypeGVCFs: 

As we did not provide commands for combining g.vcf files and their joint genotyping in separate steps, the details in step 10 are missing.

10. Combine all the g.vcf files in one 



11. Perform GenotypeGVCFs     
   java -Xmx8G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar GenotypeGVCFs \
   -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \
   --variant Vitis_raw_vcfs.list \
   --output Vitis_raw_variants_output.vcf

          OR

   java -Xmx8G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar GenotypeGVCFs \
   -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \
   --variant VRupestris.g.vcf \
   --output VRupestris_raw_SNPs_INDELs.vcf;

b) by writing bash shell having command for both the combine g.vcf and GenotypeGVCFs and executing it:
   
#!/bin/bash
filename=$1
GATK4=/B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar
REFERENCE=/B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa
LISTFILE=$PWD/raw_gvcfs_files.txt
CVCFFILE=${filename}.g.vcf
GVCFFILE=${filename}.vcf
for line in `cat $LISTFILE`
do
    {
        i=`expr $i + 1`
        variant="-V $line $variant"
    }
done
echo $variant
java -Xmx24G -jar $GATK4 CombineGVCFs -R $REFERENCE  $variant -O $CVCFFILE
java -Xmx24G -jar $GATK4 GenotypeGVCFs -R $REFERENCE  -V $CVCFFILE -O $GVCFFILE
mv $CVCFFILE $PWD/tmp
bgzip $GVCFFILE
tabix -p vcf $GVCFFILE.gz
echo "$VCFFILE successfully created."


12. SelectVariants

To select SNPs
java -Xmx64G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar SelectVariants \ 
    -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \
    -V /B42T/mukesh/DNA_Seq_Work/raw_gvcfs_files.txt.vcf.gz \ 
    --select-type-to-include SNP \ 
    --output Vitis_11_raw_snps.vcf 
To select INDELs
java -Xmx64G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar SelectVariants \ 
    -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \ 
    -V /B42T/mukesh/DNA_Seq_Work/raw_gvcfs_files.txt.vcf.gz \ 
    --select-type-to-include INDEL \
    --select-type-to-include MIXED \
    --output Vitis_11_raw_INDEL_MIXED.vcf

After we selected variants from the VCF file, now, we can perform VQSR (Variant Quality Score Recalibration). Since there were no enough SNP sites for grapevines in pubicly available databases we skipped this step, instead, we did hard filtering of variants. So, I am mentioning the HardFilter command here. 
    
13. Do HardFilter the Variants
java -Xmx64G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar VariantFiltration \
-R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa 
--variant /B42T/mukesh/DNA_Seq_Work/Vitis_11_raw_snps.vcf 
--filter-expression "QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0 || MQRankSum < -10.0 || ReadPosRankSum < -8.0 || QUAL < 30" --filter-name "my_snp_filter" 
--output Vitis_11_Filtered_snps.vcf

java -Xmx64G -jar /B42T/mukesh/gatk-4.1.0.0/gatk-package-4.1.0.0-local.jar VariantFiltration \ 
    -R /B42T/mukesh/DNA_Seq_Work/index/Vitis_vinifera.fa \ 
    --variant /B42T/mukesh/DNA_Seq_Work/Vitis_11_raw_INDEL_MIXED.vcf \ 
    --filter-expression "QD < 2.0 || FS > 100.0 || ReadPosRankSum < -20.0 || InbreedingCoeff < -0.8" \ 
    --filter-name "my_indel_filter" \ 
    --output Vitis_11_Filtered_INDEL_MIXED.vcf
14. Conduct SNPs and INDELs annotation using snpEff programme
To annotate SNPs only
java -Xmx4g -jar /B42T/mukesh/snpEff/snpEff.jar Vitis_vinifera Vitis_11_Filtered_snps.vcf > Vitis_11_Ann_SNPs.vcf 

To annotate INDEL and SNPs mixed
java -Xmx8g -jar /B42T/mukesh/snpEff/snpEff.jar Vitis_vinifera Vitis_11_Filtered_INDEL_MIXED.vcf > Vitis_11_Ann_INDEL_MIXED.vcf

References
Cingolani, P., Platts, A., Wang, L. L., Coon, M., Nguyen, T., Wang, L., … Ruden, D. M. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly, 6(2), 80–92. https://doi.org/10.4161/fly.19695

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., … DePristo, M. A. (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation  DNA sequencing data. Genome Research, 20(9), 1297–1303. https://doi.org/10.1101/gr.107524.110

Comments

Post a Comment

Popular posts from this blog

How we performed the estimation of total flavonoid content ?

-
Image
This time I am sharing a method for the estimation of total flavonoid content in our experiments. We Flowchart. the process of estimating total flavonoids content. followed a previously done standard protocol  (Quettier  et al.  2000 and Stankovic 2011)  with slight modification.  Firstly, we  mixed  e ach plant extract (2 ml; 0.5 mg/ml) separately  with aluminium chloride solution (1 ml; 2 %). Then, we allowed the reaction mixture to stand for an hour at room temperature. Again, we measured the absorbance of the mixture at 415 nm using the UV-spectrophotometer. We  obtained  calibration curve (see Fig 1.) using quercetin solution series (25, 50 75, and 100 µg/ml) as a standard. Here, we used the reagent mixture with the same volume of ethanol instead of plant extract solution as a blank. Thus, the obtained curve was calibrated to estimate the concentration of total flavonoid in each plant extract. The total ...
Read more

How to Align Fastq Files of Sample Genomic DNA Sequence to Reference Genome

-
Image
Since the aim of this blog is to share the information on various aspects of plants and plant sciences, this time, I am blogging about the bioinformatics tools that help study about the genome of organisms including plants.    The application of  BWA  software or any other software that allows alignment of fastq (.fq) files of the sample genomic DNA sequence to  reference genome  is an essential step before carrying out further investigations such as the genome analysis. In this post, I am going to present the instruction for the alignment of quality trimmed fastq (.fq) files of a sample genome to the reference genome using BWA (Burrows-Wheeler Aligner) software. This instruction covers installing of BWA software, indexing the reference genome, quality trimming the raw fastq files, and aligning the quality trimmed  fastq  files to reference genome to get  SAM  (Sequence Alignment Map) file of a sample genome.  Image 1. The vi...
Read more

The Rapid Method of Estimating Antioxidant Activity and Total Polyphenols Content of Sample Plant Extracts

-
Image
Estimation of total polyphenols content Figure1. The appearance of bluish-black colour after the addition of reagents on sample plant extracts during the estimation of total polyphenols content.  We estimated  total polyphenol content  by a standard protocol (Singleton et al . 1999 and Stankovic 2011) with slight modification. Briefly, we  separately mixed  each plant extract (0.5 ml; 1 mg/ml)  with the Folin-Ciocalteu phenol reagent (2.5 ml; 10%) and aqueous sodium bicarbonate solution (2.5 ml; 7.5%). Then, the reaction mixture was allowed to stand for about 45 minutes and the absorbance was measured at 765 nm using the UV- spectrophotometer (Thermo Fisher Scientific, Genesystem -10.5).The calibration curve was plotted using gallic acid as standard in ethanol (absolute) using the concentration series of 25, 50, 75, and 100 µg/ml. Here, the reagent mixture with the  same volume of ethanol instead of plant extract solution was used ...
Read more