Variant calling and analysis
The main steps comprising variant calling and analysis are
- mapping short reads
- calling raw variants
- adding filters (really more like 'tagging' to identify raw variants that are really variants and not technical errors) and some annotations to variants
- annotating effect(s) of variants on genes (like if they change protein sequence)
Which mapper and variant caller works best? No simple answer (of course), but see
- Hwang at al., 2015. Systematic comparison of variant calling pipelines using gold standard personal exome variants.
- Pirooznia et al., 2014. Validation and assessment of variant calling pipelines for next-generation sequencing.
- Liu et al., 2013. Variant callers for next-generation sequencing data: a comparison study.
- Note: GATK is optimized for large human datasets, whereas GATK and samtools may perform similarly with other species and smaller-scale experiments.
Map short reads
After quality control and/or filtering of reads...
Single reads with bowtie2:
bowtie2 -x /nfs/genomes/sgd_2010/bowtie/sacCer3 A_reads.fq -S A_reads.bt2.sam
Single reads with bwa:
#assuming the index file name is sacCer3.fa: use the name that was used when creating the index bwa aln /nfs/genomes/sgd_2010/bwa/sacCer3.fa A_reads.fq > A_reads.sai bwa samse /nfs/genomes/sgd_2010/bwa/sacCer3.fa A_reads.sai A_reads.fq > A_reads.bwa.sam
Paired-end reads with bowtie2:
bowtie2 -x /nfs/genomes/sgd_2010/bowtie/sacCer3 -1 A_reads.1.fq -2 A_reads.2.fq -S A_reads.1+2.bt2.sam
Paired-end reads with bwa (for reads of at least 70 nt):
bwa mem /nfs/genomes/sgd_2010/bwa/sacCer3.fa A_reads.1.fq A_reads.2.fq > A_reads.1+2.bwa.sam
Convert to BAM, sort, and index [with a custom BaRC script that uses samtools]:
/nfs/BaRC_Public/BaRC_code/Perl/SAM_to_BAM_sort_index/SAM_to_BAM_sort_index.pl A_reads.bwa.sam
Get uniquely mapping reads (Is this recommended or not?)
bsub samtools view -b -q 5 -o A_reads.bt2.unique.bam A_reads.bt2.bam
Remove duplicated reads assumed to be PCR artifacts. This should be done even if duplication rate is low.
# paired-end reads samtools rmdup A_reads.bt2.sorted.bam A_reads.bt2.sorted.noDups.bam # single-end reads samtools rmdup -s A_reads.bt2.sorted.bam A_reads.bt2.sorted.noDups.bam
Call raw variants with mpileup+bcftools
Call variants (one sample vs. reference) with samtools' mpileup+bcftools (see the samtools' variant calling workflow for more details). In our experience, "-B" (disable BAQ) or "-E" (recalculate BAQ) works better than the default method, which can remove some obvious variants.
samtools mpileup -d100000 -t DP,AD,ADF,ADR,SP,DP4,INFO/AD -Buf /nfs/genomes/sgd_2010/bwa/sacCer3.fa A_reads.bt2.sorted_unique.bam | bcftools call -O b -v -c - >| A_reads.bt2.sorted_unique.raw.bcf #in bcftools call: -m for multi-allelic, -c for bi-alleleic
Call variants (multiple sample vs. reference) using a set of BAM files
samtools mpileup -d100000 -t DP,AD,ADF,ADR,SP,DP4,INFO/AD -Buf /nfs/genomes/sgd_2010/bwa/sacCer3.fa *_reads.bt2.sorted_unique.bam | bcftools call -O b -v -c - >| ALL_reads.bt2.sorted_unique.raw.bcf
Add filters and annotations to raw variants
This step uses vcf-annotate from the VCFtools suite
Annotate variants by adding tags ("filters" but all variants are kept) to each variant, using all default filters.
bcftools view A_reads.bt2.sorted_unique.raw.bcf | vcf-annotate -f + > A_reads.bt2.sorted_unique.withTags.vcf
Generate variant statistics for each sample (with 'bcftools stats') and plot them.
bcftools stats -F <ref.fa> -s - <results.vcf> > <results.vcf.stats.txt> plot-vcfstats -p vcfstats <results.vcf.stats.txt>
Prepare file of known SNPs for use with vcf-annotate. Start with tab-delimited file (ex: SNP137.bed) that looks like
chr1 1360 1361 rs000000001
bgzip SNP137.bed tabix -p bed SNP137.bed.gz
Annotate variants by adding tags, more analysis, and any SNPdb overlaps
bcftools view A_reads.bt2.sorted_unique.raw.bcf | vcf-annotate -f +/d=10 --fill-HWE --fill-type -n -a SNP137.bed.gz -c CHROM,FROM,TO,INFO/SNP_ID -d key=INFO,ID=SNP_ID,Number=1,Type=Integer,Description='SNP137 sites' > A_reads.bt2.sorted_unique.filtered.vcf
In the output VCF file, any tags will appear in the FILTER field, and SNPdb overlaps will appear in the INFO field.
Overlap with annotations can also be identified with intersectBed (where annotation will appear in new fields of the output VCF file):
intersectBed -wao -split -a A_reads.bt2.sorted_unique.raw.vcf -b SNP137.bed > A_reads.bt2.sorted_unique.annotated.vcf
Annotate effect(s) of variants on genes
- snpEff
Use snpEff (assuming snpEff has gene+protein annotations for your genome):
snpEff -c /usr/local/share/snpEff/snpEff.config -s A_snpEff.html SacCer_Apr2011.18 A_reads.bt2.sorted_unique.filtered.vcf > A_reads.bt2.sorted_unique.filtered.snpEff.vcf
Get information only for variants that overlap protein-coding exons:
snpEff -c /usr/local/share/snpEff/snpEff.config -no-downstream -no-intergenic -no-intron -no-upstream -no-utr -s A_snpEff.html SacCer_Apr2011.18 A_reads.bt2.sorted_unique.filtered.vcf > A_reads.bt2.sorted_unique.filtered.snpEff.vcf
- Ensembl's Variant Effect Predictor (VEP)
For large number of variants, download the annotation and use --cache or --offline instead of --database, see documentation
#default human, set by --species vep -i chr1_subset.vcf --database
See Interpreting VCF files for more information.
Get read counts and statistics at desired sites
Sometimes one may have an interesting set of genome locations that may not turn up in a variant output file (if they are variants in that sample, for example).
If you need variant-type metrics for those sites, you can use a command like
samtools mpileup -uv -t DP,DPR,DV,DP4,INFO/DPR,SP -l Desired_sites.bed -f refGenome.fa Sample.bam | bcftools call -c - > Desired_sites_info.vcf
If you want counts for all 4 nucleotides (not just reference and alternate alleles), one way is with bam-readcount. Note that you'll only get a maximum coverage of ~8000 reads unless you include location (ex: chr19:50-100) as the final argument for 'bam-readcount'. The file of coordinates of your desired sites must include chr, start, end but also be 1-based (so not bed format).
bam-readcount -w 1 -l Desired_sites.txt -f refGenome.fa Sample.bam chr1 > Desired_sites.nt_counts_etc.txt # Extract just nt counts (fields: chr, position, reference_base, depth, As, Cs, Gs, Ts) perl -pe 's/:/\t/g' Desired_sites.nt_counts_etc.txt | cut -f1-4,20,34,48,62 > Desired_sites.nt_counts_only.txt
Create a consensus sequence by applying VCF variants to a reference genome file
bgzip All_samples.variants.withTags.vcf tabix -p vcf All_samples.variants.withTags.vcf.gz samtools faidx Ref_genome.fa chr1:1-1000 | bcftools consensus All_samples.variants.withTags.vcf.gz > chr1_1_1000.consensus.fa
