...
If you happen to be working with a model organism with extensive external data (ESPECIALLY HUMAN), then there are even more sophisticated tools like the Broad Institute's GATK that can improve both sensitivity and specificity of your variant calls.
Example: The CEU Trio from the 1000 Genomes Project
Many example datasets are available from the 1000 genomes project specifically for method evaluation and training. We'll explore a trio (mom, dad, child). Their accession numbers are NA12892, NA12891, and NA12878 respectively. To make the exercise run more quickly, we'll focus on data only from chromosome 20.
All the data we'll use is located here:
| Code Block | ||
|---|---|---|
| ||
$BI/ngs_course/human_variation
|
This directory contains raw data (the .fastq files), mapped data (the .bam files) and variant calls (the .vcf files). It also contains the subdirectory ref with special references.
The 1000 Genomes project is really oriented to producing .vcf filesKeep in mind that this type of trio (or familial) analysis has been exceptionally powerful for identifying rare childhood diseases. The most prominent publication in this area is this first example of whole exome sequencing saving a life. There are many other publications since and some review articles such as this one. Familial analysis is critical for rare, autosomal dominant diseases because, almost by definition, the mutations may be "private" to each individual so we can't look across big populations to find one single causative mutation. But within families, we can identify bad private mutations in single genes or pathways and then look across populations to find commonality at the gene or pathway level to explain a phenotype.
Example: The CEU Trio from the 1000 Genomes Project
Many example datasets are available from the 1000 genomes project specifically for method evaluation and training. We'll explore a trio (mom, dad, child). Their accession numbers are NA12892, NA12891, and NA12878 respectively. To make the exercise run more quickly, we'll focus on data only from chromosome 20.
All the data we'll use is located here:
| Code Block | ||
|---|---|---|
| ||
$BI/ngs_course/human_variation
|
This directory contains raw data (the .fastq files), mapped data (the .bam files) and variant calls (the .vcf files). It also contains the subdirectory ref with special references.
The 1000 Genomes project is really oriented to producing .vcf files; the file "ceu20.vcf" contains all the latest genotypes from this trio based on abundant data from the project.
...
We'll return to this example data shortly to demonstrate a much more involved tool, GATK, to do the same steps.
Note if you're trying this on Stampede: The $BI directory is not accessible from compute nodes on Stampede so you will need to make a copy of your data on $SCRATCH and update file locations accordingly to get this demo to run.
Optional: Mapping Exercise
...
| Expand | |||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| |||||||||||||||||||||||||||||||||||||||||||||
Let's use Anna Battenhouse's shell script Don't forget that for this to work, you need to have appended
Move into your scratch directory and then try to figure out how to create and
Note that the input is paired-end data.
|
Single-sample variant calling with samtools
We would normally use the BAM file from the previous mapping step to call variants in this raw data. However, for the purposes of this course we will use the actual BAM file provided by the 1000 Genomes Project (from which the .fastq file above was derived, leading to some oddities in it).$BI/ngs_course/human_variation/NA12878.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam
...
Here are the commands, piped together (ONLY run this directly if you are in an idev session - NOT on a head node!):
| Code Block | ||
|---|---|---|
| ||
samtools mpileup -uf $BI/ref_genome/fasta/ucsc/ucsc.hg19.fasta/hs37d5.fa \ $BI/ngs_course/human_variation/NA12878.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ | bcftools view -vcg - > test.raw.vcf |
...
or via qsub:
| Code Block |
|---|
launcher_creator.py -n samtools_test -b "samtools mpileup -uf ref/hs37d5.fa |
...
Exercise:
Write a modified version of the command above to use the proper reference into a commands file, create the job submission script with launcher_creator.py and submit the job.
| Expand | Solution | Solution | |||
|---|---|---|---|---|---|
As we did for mapping, we need to submit these to the Lonestar queue:
| |||||
| Expand | |||||
| Output | Output | | Code Block | NA12878.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam | bcftools view -vcg - > test.raw.vcf" -t 01:00:00 -q development -a CCBB -m samtools
qsub samtools_test.sge |
Note that the reference chosen for mpileup must be exactly the same as the reference used to create the bam file. The 1000 genomes project has created it's own reference and so the command listed above won't work - we have to use the 1000 genomes reference which is $BI/ngs_course/human_variation/ref/hs37d5.fa. We could have chosen another mapper if we'd wanted to though.
Exercise:
Write a modified version of the command above to use the proper reference into a commands file, create the job submission script with launcher_creator.py and submit the job.
| Expand | |||||
|---|---|---|---|---|---|
| |||||
As we did for mapping, we need to submit these to the Lonestar queue:
|
| Expand | ||||
|---|---|---|---|---|
| ||||
|
After your single-sample variant calling job completes
| Expand | ||||
|---|---|---|---|---|
| ||||
If you don't want to wait, CHANGE TO A NEW DIRECTORY and do this: |
...
You can also get some quick stats with some linux one-liners on this page; there are more thorough analysis programs built to work with vcf's.
...
To instruct samtools to call variants across many samples, you must simply give it mapped data with each sample tagged separately. Samtools allows two methods to do this:
By providing separate bam files for each sample, like this:
Code Block title samtools multi-sample variants: separate bam files samtools mpileup -uf hs37d5.fa \ NA12878.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ NA12891.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ NA12892.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ | bcftools view -vcg - > all.samtools.vcfBy providing one or more bam files, each containing mapped reads from multiple samples tagged with unique samtools
@RGtags.Code Block title samtools multi-sample variants: one or more bam files using @RG samtools mpileup -uf hs37d5.fa all.bam | bcftools view -vcg - > all.raw.vcf
The output file from the first option is in $BI/ngs_course/human_variation
CAVEAT: If you intend to use the second option, you must make sure you tag your reads with the right RG tag; this can easily be done during the samse or sampe stage of mapping with bwa with the -r option, using samtools merge, with picard tools AddOrReplaceReadGroup command, or with your own perl/python/bash commands. Note that our align_bwa.sh script takes care of this detail for us. Samtools allows two methods to do this:
By providing separate bam files for each sample, like this:
Code Block title samtools multi-sample variants: separate bam files samtools mpileup -uf ref/hs37d5.fa \ NA12878.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ NA12891.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ NA12892.chrom20.ILLUMINA.bwa.CEU.exome.20111114.bam \ | bcftools view -vcg - > all.samtools.vcfBy providing one or more bam files, each containing mapped reads from multiple samples tagged with unique samtools
@RGtags.Code Block title samtools multi-sample variants: one or more bam files using @RG samtools mpileup -uf hs37d5.fa all.bam | bcftools view -vcg - > all.raw.vcf
The output file from the first option is in $BI/ngs_course/human_variation
CAVEAT: If you intend to use the second option, you must make sure you tag your reads with the right RG tag; this can easily be done during the samse or sampe stage of mapping with bwa with the -r option, using samtools merge, with picard tools AddOrReplaceReadGroup command, or with your own perl/python/bash commands. Note that our align_bwa.sh script takes care of this detail for us.
Identify the lineage
If genetics works, you should be able to identify the child based strictly on the genotypes. Can you do it?
| Expand | ||
|---|---|---|
| ||
You're trying to find the genotypes in the |
| Expand | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| |||||||||
This linux one-liner should give you a snapshot of data sufficient to figure it out:
|
GATK
GATK deserves it's own page which is here, but we've already run it and will now look at some differences in the SNP calls.
Look at the differences between single- and multiple-sample SNP calls, and between Samtools/Bcftools SNP calls and GATK SNP calls
How many SNPs were called in each case?
| Expand | ||
|---|---|---|
| ||
for file in `ls *.vcf`; do echo "File: $file `cat $file | grep -v '^#' | wc -l`"; done |
What's the overlap between all the single-sample SNP calls aggregated together and the multi-sample SNP calls?
| Expand | ||
|---|---|---|
|
In theory, GATK does a much better job ruling out false positives. But are there some SNPs GATK calls with high confidence that Samtools doesn't call at all?
| Expand | ||
|---|---|---|
| ||
What's going on here: 1) Yank out all the variant calls (comments start with '#') and add a string of "AAAAAAA" to them to make "sort" do it's job in a way that "join" will later like 2) Join the two files using their chromosome position as the join field, but also include any lines from the GATK file that DON'T match the samtools file. Use "awk" to figure out which ones came only from GATK (they are missing a bunch of fields from the samtools variant calls), look only for those that GATK has labeled "PASS" and write them to a file. 3) Sort the resultant file on the variant quality value - take the top 10 lines. You will note that many of these are complex variants, particularly insertions, so it's not too surprising that GATK does better. But here's a SNP that GATK does much better on: chr21:34278313 It has an interesting quantitative signature though... you might want to look at it in IGV. |
Other notes on bcftools
bcftools has many other sub-commands such as performing association tests and estimating allele frequency spectrums.
You can't get gene-context information from bcftools - you have to shift to variant annotation tools to do that.
Side-note on samtools mpileup
If you don't need a variant caller that uses a full Bayesian inference engine with robust filtering, protections, etc. but just want to look for weird base-artifacts within your mapped data, check out samtools mpileup output directly.
GATK: A much more sophisticated option
http://gatkforums.broadinstitute.org/discussion/44/base-quality-score-recalibration-bqsr