3,045
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchSeqShop: Variant Calling and Filtering for SNPs Practical, June 2014 (view source)
Revision as of 14:37, 13 November 2014
, 14:37, 13 November 2014no edit summary
'''Note:''' the latest version of this practical is available at: [[SeqShop: Variant Calling and Filtering for SNPs Practical]]
* The ones here is the original one from the June workshop (updated to be run from elsewhere)
==Introduction==
See the [[Media:SeqShop - GotCloud snpcall.pdf|introductory slides]] for an intro to this tutorial.
== Goals of This Session ==
== Goals of This Session ==
* What we want to learn
* What we want to learn
** How to examine the variants at particular genomic positions
** How to examine the variants at particular genomic positions
** How to evaluate the quality of SNP calls
** How to evaluate the quality of SNP calls
== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
{{SeqShopLogin}}
== Setup your run environment==
=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in. It will setup some environment variables to point you to:
This will setup some environment variables to point you to
* GotCloud program
* [[GotCloud]] program
* Tutorial input files
* Tutorial input files
* Setup an output directory
* Setup an output directory
* You won't see any output after running <code>source</code>
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /home/mktrost/seqshop/setup.txt
and press 'q' to finish.
<div class="mw-collapsible mw-collapsed" style="width:200px">
<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
View setup.txt
</div>
</div>
</div>
</div>
</div>
</div>
== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">
This tutorial builds on the alignment tutorial, if you have not already, please first run that tutorial: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical, June 2014|Alignment Tutorial]]
{{SeqShopRemoteEnv}}
== Examining GotCloud SnpCall Input files ==
== Examining GotCloud SnpCall Input files ==
Per sample BAM files contain sequence reads that are mapped to positions in the genome.
Per sample BAM files contain sequence reads that are mapped to positions in the genome.
For a reminder on how to look at/read BAM files, see: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#BAM_Files|SeqShop Aligment: BAM Files]]
For a reminder on how to look at/read BAM files, see: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical, June 2014#BAM_Files|SeqShop Aligment: BAM Files]]
For this tutorial, we will use the 4 BAMs produced in the [[SeqShop: Sequence Mapping and Assembly Practical, June 2014]] as well as with 58 BAMs that were pre-aligned to that 1MB region of chromosome 22.
=== Reference Files ===
=== Reference Files ===
Reference files can be downloaded with GotCloud or from other sources
Reference files can be downloaded with GotCloud or from other sources.
* For this practical, I already downloaded them for you.
* See [[GotCloud: Genetic Reference and Resource Files]] for more information on downloading/generating reference files
* See [[GotCloud: Genetic Reference and Resource Files]] for more information on downloading/generating reference files
#* Contains the reference base for each position of each chromosome
#* Contains the reference base for each position of each chromosome
#** Used to compare bases in sequence reads to the reference positions they mapped to
#** Used to compare bases in sequence reads to the reference positions they mapped to
#** Identify SNPs
#** Used to identify SNPs/variations in the sequence reads
#* Additional information on the FASTA format: http://en.wikipedia.org/wiki/FASTA_format
#* Additional information on the FASTA format: http://en.wikipedia.org/wiki/FASTA_format
# VCF (variant call format) files with chromosomes/positions
# VCF (variant call format) files with chromosomes/positions
#* dbsnp - used for generating summary statistics
#* dbsnp - used for generating summary statistics
We looked at them yesterday, but you can take another look at the chromosome 22 reference files included for this tutorial:
ls ${SS}/ref22
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:RefDir.png|700px]]
</div>
</div>
</ul>
=== GotCloud BAM Index File ===
=== GotCloud BAM Index File ===
The BAM index file points GotCloud to the BAM files
The BAM index file points GotCloud to the BAM files
* Alignment pipeline generates for you
* generated by the alignment pipeline
Look at the BAM index file the alignment pipeline generated
Look at the BAM index file the alignment pipeline generated
The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.
The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.
* The other 58 samples were already aligned:
* The other 58 samples were already aligned:
ls ${IN}/bams
ls ${SS}/bams
Look at the BAM index for those BAMs:
Look at the BAM index for those BAMs:
less ${IN}/bams/bam.index
less ${SS}/bams/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
Remember, use <code>'q'</code> to exit out of <code>less</code>
;Do you notice a difference between this index and yours?
;Do you notice a difference between this index and yours?
<ul>
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<div class="mw-collapsible mw-collapsed" style="width:550px">
<li>Answer:</li>
<li>Answer:</li>
<div class="mw-collapsible-content">
<div class="mw-collapsible-content">
<li>It doesn't have a full path to the BAM file, while your index has /home/...</li>
<li>It doesn't have a full path to the BAM file, while your index has /home/...</li>
[[File:Bamindex1.png|300px]]
[[File:Bamindex1.png|300px]]
<li>That's ok, <code>gotcloud.conf</code> contains the path to those BAMs</li>
<li>That's ok, we will use the <code>--base_prefix ${SS}</code> command-line option to prefix the BAM paths</li>
<li>Alternatively, we could have set BAM_INDEX in <code>gotcloud.conf</code> to the path to the BAMs
<pre>BAM_INDEX = /home/username/seqshop/example</pre> </li>
<ul>
<li>NOTE: the conf file can't interpret ${SS} environment variables or '~', so you would have to specify the full path</li>
<li>We just used the command-line option for this tutorial since this path will vary by user.</li>
</ul>
</div>
</div>
</div>
</div>
We need to add these BAMs to our index
We need to add these BAMs to our index
* Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
* Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
cat ${IN}/bams/bam.index >> ${OUT}/bam.index
** '''Be sure to do this command just once'''
cat ${SS}/bams/bam.index >> ${OUT}/bam.index
* ">>" will append to the file that follows it
* ">>" will append to the file that follows it
** Check that your BAM index is the correct size
**:<pre>wc -l ${OUT}/bam.index</pre>
*** <code>wc -l</code> counts the number of lines in the file
*** Should be 62
Verify your BAM index contains the additional BAMs
Verify your BAM index contains the additional BAMs
less ${IN}/bams/bam.index
less ${OUT}/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
Remember, use <code>'q'</code> to exit out of <code>less</code>
</ul>
</ul>
</ul>
</ul>
=== GotCloud Configuration File ===
We will use the same configuration file as we used yesterday in GotCloud Align.
See [[SeqShop:_Sequence Mapping and Assembly Practical, June 2014#GotCloud Configuration File|SeqShop: Alignment: GotCloud Configuration File]] for more details
* Note we want to limit snpcall to just chr22 so the configuration already has <code>CHRS = 22</code> (default was 1-22 & X).
For more information on configuration, see: [[GotCloud:_Variant_Calling_Pipeline#Configuration_File|GotCloud snpcall: Configuration File]]
* Contains information on how to configure for exome/targeted sequencing
== Run GotCloud SnpCall ==
== Run GotCloud SnpCall ==
Now that we have all of our input files, we need just a simple command to run:
Now that we have all of our input files, we need just a simple command to run:
${GC}/gotcloud snpcall --conf ${IN}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000
${GC}/gotcloud snpcall --conf ${SS}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000 --base_prefix ${SS} --outdir ${OUT}
* <code>${GC}/gotcloud</code> runs GotCloud
* <code>align</code> tells GotCloud you want to run the alignment pipeline.
* <code>--conf</code> tells GotCloud the name of the configuration file to use.
** The configuration for this test was downloaded with the seqshop input files.
* --numjobs tells GotCloud how many jobs to run in parallel
* --numjobs tells GotCloud how many jobs to run in parallel
** Depends on your system
** Depends on your system
* --region 22:36000000-37000000
* --region 22:36000000-37000000
** The sample files are just a small region of chromosome 22, so to save time, we tell GotCloud to ignore the other regions
** The sample files are just a small region of chromosome 22, so to save time, we tell GotCloud to ignore the other regions
* <code>--base_prefix</code> tells GotCloud the prefix to append to relative paths.
** The Configuration file cannot read environment variables, so we need to tell GotCloud the path to the input files, ${SS}
** Alternatively, gotcloud.conf could be updated to specify the full paths
* <code>--out_dir</code> tells GotCloud where to write the output.
** This could be specified in gotcloud.conf, but to allow you to use the ${OUT} to change the output location, it is specified on the command-line
<div class="mw-collapsible mw-collapsed" style="width:500px">
<div class="mw-collapsible mw-collapsed" style="width:500px">
If you cancelled GotCloud part way through, just rerun your GotCloud command and it will pick up where it left off.
If you cancelled GotCloud part way through, just rerun your GotCloud command and it will pick up where it left off.
If you want to understand more detailed step of GotCloud SNP calling, here is a schematic picture with a little bit more details
[[File:Gotcloudoverview.png|600px]]
== Examining GotCloud SnpCall Output ==
== Examining GotCloud SnpCall Output ==
<div class="mw-collapsible mw-collapsed" style="width:250px">
<div class="mw-collapsible mw-collapsed" style="width:250px">
View Annotated Screenshot
View Screenshot
<div class="mw-collapsible-content">
<div class="mw-collapsible-content">
[[File:filterSum.png]]
[[File:filterSum.png]]
</div>
</div>
'''To understand how to interpret the filtering summary statistics, please refer to [[Understanding vcf-summary output]]'''
=== Filtered VCF ===
=== Filtered VCF ===
Now, let's look at the filtered file with genotypes.
Now, let's look at the filtered file with genotypes.
zless -S ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz
* Scroll down until you find some variants.
* Scroll down until you find some variants.
** Use down arrow to move down one line
** Use down arrow to move down one line
* Scroll right until you should see per sample genotype information
* Scroll right until you should see per sample genotype information
Remember, use <code>'q'</code> to exit out of <code>less</code>
q
<div class="mw-collapsible mw-collapsed" style="width:250px">
<div class="mw-collapsible mw-collapsed" style="width:250px">
<div class="mw-collapsible-content">
<div class="mw-collapsible-content">
[[File:SvmFiltGL.png|550px]]
[[File:SvmFiltGL.png|550px]]
</div>
</div>
=== Passing SNPs ===
Let's look at the file of just the pass sites:
zless -S ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz
* Scroll down: they all look like they <code>PASS</code>
Remember, use 'q' to exit out of less
q
Let's check if they are all PASS.
zcat ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
It will return nothing since there are no non-passing variants in this file.
<div class="mw-collapsible mw-collapsed" style="width:450px">
;Want an explanation of this command?
<div class="mw-collapsible-content">
* zcat ...: uncompress the zipped VCF
* '|' : this takes the output of one command and sends it as input to the next
* grep -v "^#" : exclude any lines that start with "#" - headers
* cut -f 7 : extract the FILTER column (the 7th column)
* grep -v "PASS" : exclude any rows that have a "PASS" in the FILTER column
</div>
</div>
Compare that to the filtered file we looked at before:
zcat ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
;Do you see any filters?
<div class="mw-collapsible mw-collapsed" style="width:450px">
*Answer
<div class="mw-collapsible-content">
* Yes
** It should have scrolled and you should see filters like:
*** INDEL5;SVM
*** INDEL5
*** SVM
</div>
</div>
</div>
</div>
== GotCloud Genotype Refinement ==
== GotCloud Genotype Refinement ==
To improve the quality of the genotypes, we run a genotype refinement pipeline.
This pipeline runs [http://faculty.washington.edu/browning/beagle/beagle.html Beagle] & thunder.
=== Genotype Refinement Input ===
The GotCloud genotype refinement pipeline takes as input ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz (the VCF file of PASS'ing SNPs from snpcall).
The bam index and the configuration file we used for GotCloud snpcall will tell GotCloud genotype refinement everything it needs to know, so no new input files need to be prepared.
Note: the configuration file overrides the THUNDER command to make it go faster than the default settings so the tutorial will run faster:
[[File:thunderConf.png|600px]]
=== Running GotCloud Genotype Refinement ===
Since everything is setup, just run the following command (very similar to snpcall).
${GC}/gotcloud ldrefine --conf ${SS}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000 --base_prefix ${SS} --outdir ${OUT}
* Beagle will take about 2-3 minutes to complete
* Thunder will automatically run and will take another 3-4 minutes
<div class="mw-collapsible mw-collapsed" style="width:350px">
When completed, it should look like this:
<div class="mw-collapsible-content">
[[File:GcldrefineOut.png]]
</div>
</div>
=== Genotype Refinement Output ===
; What's new in the output directory?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer</li>
<div class="mw-collapsible-content">
:<pre>ls ${OUT}</pre>
<ul>
<li><code>beagle</code> directory : Beagle output</li>
<li><code>thunder</code> directory : Thunder output</li>
<li><code>umake.beagle.*</code> : Contain the configuration & steps used in GotCloud beagle</li>
<li><code>umake.thunder.*</code> files : Contain the configuration & steps used in GotCloud thunder</li>
</ul>
</div>
</div>
</ul>
Let's take a look at that interesting location we found in the [[SeqShop:_Sequence_Mapping_and_Assembly_Practical, June 2014#Accessing_BAMs_by_Position|alignment tutorial]] : chromosome 22, positions 36907000-36907100
Use tabix to extract that from the VCFs:
${GC}/bin/tabix ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
Remember, type 'q' to quit less.
q
;Are there any variants in this region?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>Yes!</li>
<li>Positions:</li>
<ul>
<li><code>36907001</code>; Ref: T, Alt: C - that's what we saw before</li>
<li><code>36907098</code>; Ref: T, Alt: C - that's what we saw before</li>
</ul>
</ul>
</div>
</div>
</ul>
;What is HG00551's genotype at these positions?
#First check which sample number HG00551 is:
zcat ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz |grep "#CHROM"
* That will help you figure out it's genotype.
* Rerun the tabix command and scroll to find HG00551's genotype:
${GC}/bin/tabix ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>It is the first sample</li>
<li><code>0|1</code>: Heterozygous</li>
<li><code>1|1</code>; Homozygous Alt (C)</li>
</ul>
</div>
</div>
</ul>
Remember, type 'q' to quit less.
q
=== Did I find interesting variants? ===
The region we selected contains ''APOL1'' gene, which is known to play an important role in kidney diseases such as nephrotic syndrome. One of the non-synonymous risk allele, <code>rs73885139</code> located at position <code>22:36661906</code> increases the risk of nephrotic syndrome by >2-folds. Let's see if we found the interesting variant by looking at the VCF file by position.
${GC}/bin/tabix ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1
Did you see a variant at the position?
${GC}/bin/tabix ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;
FIC=-0.0092;SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;
OBS=145,191,0,0,3,2,0,0;STR=-0.040;STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;
AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;MQ30=0.000;SVM=1.51214
GT:DP:GQ:PL 0/0:4:28:0,12,65
Let's check the sequence data to confirm that the variant really exists
${GC}/bin/samtools tview ${SS}/bams/HG01242.recal.bam ${SS}/ref22/human.g1k.v37.chr22.fa
* Type 'g' to go to a specific position
* Type 22:36661906 to move to the position
* Press arrows to move between positions
* Press 'b' if you want to color by base quality
* Press '?' for more help
<div class="mw-collapsible mw-collapsed" style="width:250px">
View Screenshot
<div class="mw-collapsible-content">
[[File:Samtoolstviewsnp.png|600px]]
</div>
</div>
=== Improvements ===
Let's get some information on the BEAGLE VCF:
perl ${SS}/ext/bed-diff.pl --vcf1 ${SS}/ref22/1kg.omni.chr22.36Mb.vcf.gz --vcf2 ${OUT}/beagle/chr22/chr22.filtered.PASS.beagled.ALL.vcf.gz --gcRoot ${GC} --out ${OUT}/bedDiff.beagle
Look at the results:
more ${OUT}/bedDiff.beagle.summary
<div class="mw-collapsible mw-collapsed" style="width:400px">
*Results
<div class="mw-collapsible-content">
OVERALL: 43601 44293 0.9844
NREF-EITHER: 19667 20359 0.9660
NMAJ-EITHER: 14585 15277 0.9547
HOMREF: 23934 100 1 0.9958
HET: 329 11959 175 0.9596
HOMALT: 4 83 7708 0.9888
HOMMAJ: 29016 126 2 0.9956
HET: 364 11959 140 0.9596
HOMMIN: 3 57 2626 0.9777
</div>
</div>
Now, let's see if it improved after running Thunder VCF:
perl ${SS}/ext/bed-diff.pl --vcf1 ${SS}/ref22/1kg.omni.chr22.36Mb.vcf.gz --vcf2 ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --gcRoot ${GC} --out ${OUT}/bedDiff.thunder
Look at the results:
more ${OUT}/bedDiff.thunder.summary
<div class="mw-collapsible mw-collapsed" style="width:400px">
*Results
<div class="mw-collapsible-content">
OVERALL: 43685 44293 0.9863
NREF-EITHER: 19758 20366 0.9701
NMAJ-EITHER: 14688 15296 0.9603
HOMREF: 23927 106 2 0.9955
HET: 286 12057 120 0.9674
HOMALT: 6 88 7701 0.9879
HOMMAJ: 28997 144 3 0.9950
HET: 286 12057 120 0.9674
HOMMIN: 5 50 2631 0.9795
</div>
</div>
There is an improvement.
== What is GotCloud snpcall doing? ==
To run GotCloud, you really just needed a single command.
Well, that one command runs many steps. Here is a diagram of all the steps.
[[File:SnpcallPipeline.jpg|1200px]]
Aren't you glad you didn't have to configure & run each one yourself?
