Difference between revisions of "SeqShop: Analysis of Structural Variation Practical, June 2014"

From Genome Analysis Wiki
Jump to navigationJump to search
(Created page with "==Introduction== See the introductory slides for an intro to this tutorial. == Goals of This Session == * What we want to learn is ca...")
 
 
(59 intermediate revisions by 2 users not shown)
Line 1: Line 1:
−
==Introduction==
+
'''Note:''' the latest version of this practical is available at: [[SeqShop: Analysis of Structural Variation Practical]]
−
See the [[Media:SeqShop - GotCloud svcall.pdf|introductory slides]] for an intro to this tutorial.
+
* The ones here is the original one from the June workshop (updated to be run from elsewhere)
  
  
Line 9: Line 9:
 
** How to perform genotyping for large deletions
 
** How to perform genotyping for large deletions
 
** How to perform variant discovery and filtering from third party sites.
 
** How to perform variant discovery and filtering from third party sites.
 +
 +
Please refer to [[Media:Seqshop cnv partb 2014 06.pdf|Lecture slides]] for more general background.
 +
 +
== GenomeSTRiP ==
 +
GenomeSTRiP was developed at the Broad Institute and at the McCarroll Lab at the Harvard Medical School Department of Genetics: http://www.broadinstitute.org/software/genomestrip/
 +
 +
If you use GenomeSTRiP for your research, please cite it:
 +
Handsaker RE, Korn JM, Nemesh J, McCarroll SA
 +
Discovery and genotyping of genome structural polymorphism by sequencing on a population scale.
 +
Nature genetics 43, 269-276 (2011)
 +
PMID: 21317889
 +
 +
GenomeStrip is currently included in with the seqshop example data under the svtoolkit directory.  We have added the bin/ sub-directory to add a high level pipeline that will run genomestrip in the same framework as GotCloud.
 +
 +
== 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
Line 22: Line 42:
 
* 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 set up, type
+
** If you want to view the detail of the setup, type
 
  less /home/mktrost/seqshop/setup.txt
 
  less /home/mktrost/seqshop/setup.txt
 
and press 'q' to finish.
 
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 class="mw-collapsible-content">
 
<div class="mw-collapsible-content">
 
[[File:setup.png|500px]]
 
[[File:setup.png|500px]]
 +
</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]]
 +
 +
It also uses the bam.index file created in the SnpCall Tutorial.  If you have not yet run that tutorial, please follow the directions at: [[SeqShop:_Variant_Calling_and_Filtering_for_SNPs_Practical, June 2014#GotCloud_BAM_Index_File|GotCloud BAM Index File]]
 +
 +
 +
{{SeqShopRemoteEnv}}
 
</div>
 
</div>
 
</div>
 
</div>
Line 34: Line 73:
 
== Examining GotCloud/GenomeSTRiP Input files ==
 
== Examining GotCloud/GenomeSTRiP Input files ==
  
−
=== Sequnce Alignment Files: BAM Files ===
+
=== Sequnce Alignment Files: BAM Files and Index Files===
 +
 
 +
The GotCloud GenomeSTRiP structural variant caller takes the same inputs as GotCloud snpcall.
 +
* BAMs->SVs rather than BAMs->SNPs
 +
 
 +
If you want a reminder, of what they look like, here is a link to the previous tutorial : [[SeqShop:_Variant_Calling_and_Filtering_for_SNPs_Practical, June 2014#Examining_GotCloud_SnpCall_Input_files|GotCloud SnpCall Input Files]]
 +
 
 +
If you want to check if you still have the bam index file, run
 +
 
 +
head ${OUT}/bam.index
 +
 
 +
<ul>
 +
<div class="mw-collapsible mw-collapsed" style="width:200px">
 +
<li>View Results</li>
 +
<div class="mw-collapsible-content" style="width:800px">
 +
HG00641 ALL /net/seqshop-server/hmkang/out/bams/HG00641.recal.bam
 +
HG00640 ALL /net/seqshop-server/hmkang/out/bams/HG00640.recal.bam
 +
HG00551 ALL /net/seqshop-server/hmkang/out/bams/HG00551.recal.bam
 +
HG00553 ALL /net/seqshop-server/hmkang/out/bams/HG00553.recal.bam
 +
HG00554 ALL bams/HG00554.recal.bam
 +
HG00637 ALL bams/HG00637.recal.bam
 +
HG00638 ALL bams/HG00638.recal.bam
 +
HG00734 ALL bams/HG00734.recal.bam
 +
HG00736 ALL bams/HG00736.recal.bam
 +
HG00737 ALL bams/HG00737.recal.bam
 +
</div>
 +
</div>
 +
</ul>
 +
 
 +
Also, make sure that you have only 62 samples (you did not append new files twice)
  
−
Per sample BAM files contain sequence reads that are mapped to positions in the genome.
+
wc -l ${OUT}/bam.index
  
−
For a reminder on how to look at/read BAM files, see: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#BAM_Files|SeqShop Aligment: BAM Files]]
+
Your expected output is similar to this.
  
−
For this tutorial, we will use the 4 BAMs produced in the [[SeqShop: Sequence Mapping and Assembly Practical]] as well as with 58 BAMs that were pre-aligned to that 1MB region of chromosome 22.
+
62 /net/seqshop-server/hmkang/out/bam.index
  
 
=== Reference Files ===
 
=== Reference Files ===
Line 47: Line 115:
 
* 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
  
−
For GotCloud/GenomeSTRiP, you need:
+
Similar to SNP and Indel calling, you need
 
# Reference genome FASTA file
 
# Reference genome FASTA file
 +
 +
For running GenomeSTRiP, you additionally need:
 
# Masked FASTA file to exclude hard-to-align regions
 
# Masked FASTA file to exclude hard-to-align regions
 
# PloidyMap file indicating the regions of genomes with unusual ploidy (e.g. chrX, chrY)
 
# PloidyMap file indicating the regions of genomes with unusual ploidy (e.g. chrX, chrY)
  
−
We looked at them yesterday, but you can take another look at the chromosome 22 reference files included for this tutorial:
+
We looked at them in previous tutorials, but you can take another look at the chromosome 22 reference files included for this tutorial:
−
  ls ${REF}
+
 
 +
  ls ${SS}/ref22
  
 
<ul>
 
<ul>
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
−
<li>View Screenshot</li>
+
<li>View Results</li>
−
<div class="mw-collapsible-content">
+
<div class="mw-collapsible-content" style="width:800px">
 
[[File:RefDir.png|500px]]
 
[[File:RefDir.png|500px]]
 
</div>
 
</div>
Line 64: Line 135:
 
</ul>
 
</ul>
  
−
=== GotCloud BAM Index File ===
 
−
The BAM index file points GotCloud to the BAM files
 
−
* generated by the alignment pipeline
 
  
−
Look at the BAM index file the alignment pipeline generated
+
Additional Reference files required just for Structural Variation:
−
  cat ${OUT}/bam.index
+
  ls ${SS}/svtoolkit/ref
−
 
 
−
;What is the path to the BAM file for sample HG00640?
 
 
<ul>
 
<ul>
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
−
<li>Answer:</li>
+
<li>View Results</li>
−
<div class="mw-collapsible-content">
+
<div class="mw-collapsible-content" style="width:800px">
−
<ul>
+
human_g1k_v37.chr22.mask.100.fasta  human_g1k_v37.chr22.mask.100.fasta.dict  human_g1k_v37.chr22.mask.100.fasta.fai
−
<li>/home/YourUserName/out/bams/HG00640.recal.bam</li>
 
−
[[File:Bamindex.png|500px]]
 
 
</div>
 
</div>
 
</div>
 
</div>
 
</ul>
 
</ul>
−
</ul>
 
−
 
−
The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.
 
−
* The other 58 samples were already aligned:
 
−
ls ${IN}/bams
 
  
−
Look at the BAM index for those BAMs:
 
−
less ${IN}/bams/bam.index
 
  
−
Remember, use <code>'q'</code> to exit out of <code>less</code>
+
Parameters files required just for Structural Variation:
−
q
+
ls ${SS}/svtoolkit/conf
  
−
;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:200px">
−
<li>Answer:</li>
+
<li>View Results</li>
−
<div class="mw-collapsible-content">
+
<div class="mw-collapsible-content" style="width:800px">
−
<ul>
+
genstrip_parameters.txt  humgen_g1k_v37_ploidy.chr22.map  humgen_g1k_v37_ploidy.map
−
<li>It doesn't have a full path to the BAM file, while your index has /home/...</li>
 
−
[[File:Bamindex1.png|300px]]
 
−
<li>That's ok, <code>gotcloud.conf</code> contains the path to those BAMs</li>
 
−
[[File:BamindexConf.png|300px]]
 
 
</div>
 
</div>
 
</div>
 
</div>
−
</ul>
 
 
</ul>
 
</ul>
  
 +
=== GotCloud Configuration File ===
 +
We will use the same configuration file we used for the GotCloud Align tutorial.
  
−
We need to add these BAMs to our index
+
See [[SeqShop:_Sequence_Mapping_and_Assembly_Practica, June 2014l#GotCloud Configuration File|SeqShop: Alignment: GotCloud Configuration File]] for more details
−
* Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
+
* Note we want to limit snpcall to just chr22 so the configuration already has <code>CHRS = 22</code> (default was 1-22 & X).
−
cat ${IN}/bams/bam.index >> ${OUT}/bam.index
 
−
* ">>" will append to the file that follows it
 
  
−
Verify your BAM index contains the additional BAMs
+
For more information on configuration, see: [[GotCloud:_Variant_Calling_Pipeline#Configuration_File|GotCloud snpcall: Configuration File]]
−
  less ${OUT}/bam.index
 
  
−
Remember, use <code>'q'</code> to exit out of <code>less</code>
+
Check out the GenomeStrip specific settings at the end of the configuration file
−
  q
+
  tail -n 8 ${SS}/gotcloud.conf
  
−
;Do you see both sets of BAMs?
 
 
<ul>
 
<ul>
−
<div class="mw-collapsible mw-collapsed" style="width:500px">
+
<div class="mw-collapsible mw-collapsed" style="width:200px">
−
<li>Annotated Screenshot:</li>
+
<li>View Results</li>
−
<div class="mw-collapsible-content">
+
<div class="mw-collapsible-content" style="width:800px">
−
<ul>
+
##############################
−
<li>If not, let me know</li>
+
## GenomeSTRIP
−
[[File:Bamindex2.png|400px]]
+
#############################
 +
GENOMESTRIP_OUT = $(OUT_DIR)/sv
 +
GENOMESTRIP_SVTOOLKIT_DIR = svtoolkit
 +
GENOMESTRIP_MASK_FASTA = $(GENOMESTRIP_SVTOOLKIT_DIR)/ref/human_g1k_v37.chr22.mask.100.fasta
 +
GENOMESTRIP_PLOIDY_MAP = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/humgen_g1k_v37_ploidy.chr22.map
 +
GENOMESTRIP_PARAM = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/genstrip_parameters.txt
 
</div>
 
</div>
 
</div>
 
</div>
−
</ul>
 
 
</ul>
 
</ul>
  
−
=== GotCloud Configuration File ===
+
== Before starting... a few 'why' questions.. ==
−
We will use the same configuration file as we used yesterday in GotCloud Align.
 
  
−
See [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#GotCloud Configuration File|SeqShop: Alignment: GotCloud Configuration File]] for more details
+
=== Why use GenomeSTRiP?===
−
* Note we want to limit snpcall to just chr22 so the configuration already has <code>CHRS = 22</code> (default was 1-22 & X).
+
# GenomeSTRiP is a mature software for detecting and genotyping large deletions (and duplications soon to be implemented). In 1000 Genomes, GenomeSTRiP was demonstrated as one of the top-performing SV caller in most evaluation metrics.
 +
# GenomeSTRiP is a great tool to integrate across multiple structural variant calls. When multiple structural variant calls exists, all the other variants can be genotyped and filtered with GenomeSTRiP, and that is how 1000 Genomes structural variant call sets were made.
 +
# Currently, GenomeSTRiP only allows calling large deletions, but duplicate calling pipeline is under way.
  
−
For more information on configuration, see: [[GotCloud:_Variant_Calling_Pipeline#Configuration_File|GotCloud snpcall: Configuration File]]
+
=== Why do we use GotCloud/GenomeSTRiP pipeline? ===
−
* Contains information on how to configure for exome/targeted sequencing
+
# The main purpose of GotCloud pipelines is to provide a pipeline for users with limited knowledge and experience with high performance computing environment.
 +
#* GotCloud/GenomeSTRiP provide a simple interface consistent to alignment, SNP, and indel calling.
 +
#* GenomeSTRiP itself also provides a straightforward pipeline to use as standalone software
 +
# GotCloud supports a variety of cluster environment that is not currently supported by GenomeSTRiP
 +
#* GenomeSTRiP is designed based on a framework called Qscript, which provide a nice support for LSF cluster system
 +
#* GotCloud support many additional cluster environments such as MOSIX or SLURM we use locally at Michigan.
 +
# GotCloud also provide a fault-tolerant solution for large-scale jobs.
 +
#* GotCloud automatically picks up jobs from the point where it failed. This allows easier and simpler run against potential technical glitches in the system.
  
−
== Run GotCloud SnpCall ==
+
== Overview of GotCloud/GenomeSTRiP pipeline ==
−
[[File:SnpcallDiagram.png|500px]]
+
GotCloud/GenomeSTRiP pipeline consists of three separate steps.
  
−
Now that we have all of our input files, we need just a simple command to run:
+
* Preprocess step : Create metadata summarizing the GC profiles, depth distribution, insert size distribution for accurate discovery and genotyping of structural variants.
−
${GC}/gotcloud snpcall --conf ${IN}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000
+
* Discovery step : Perform variant discovery split by region, across all samples. Also, perform variant filtering based on expert knowledge.
−
* --numjobs tells GotCloud how many jobs to run in parallel
+
* Genotyping step : Iterate discovered variants across the samples and calculate the genotype likelihood of for each possible genotype.
−
** Depends on your system
 
−
* --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
 
  
−
<div class="mw-collapsible mw-collapsed" style="width:500px">
+
In addition, if one wants to genotype structural variants from other structural variant caller, there is a step available.
−
Curious if it started running properly?  Check out this screenshot:
+
* Third-party Genotyping and Filtering step : Perform genotyping on the variant sites specified by an input VCF, and also perform variant filtering.
−
<div class="mw-collapsible-content">
 
−
[[File:SnpcallStart.png|750px]]
 
−
</div>
 
−
</div>
 
−
This should take about 5 minutes to run.
 
−
* After about 4 minutes of running, GotCloud snpcall will output some text to the screen.  Don't worry, that is expected and is just output from some of the intermediate tools.
 
−
* It should end with a line like: <code>Commands finished in 289 secs with no errors reported</code>
 
  
−
If you cancelled GotCloud part way through, just rerun your GotCloud command and it will pick up where it left off.
+
== Running GotCloud/GenomeSTRiP Metadata Pipeline ==
  
−
== Examining GotCloud SnpCall Output ==
+
We first need to create metadata summarizing genome-wide statistics such as GC profiles, depth distribution, insert size distributions.
−
Let's look at the output directory:
 
−
ls ${OUT}
 
  
−
;Do you see any new files or directories?
+
In principle, the metadata can be created from the input BAM files by running the following command
−
<div class="mw-collapsible mw-collapsed" style="width:350px">
+
perl ${SS}/svtoolkit/bin/genomestrip.pl -run-metadata --conf ${SS}/gotcloud.conf --numjobs 2 --base-prefix ${SS} --outdir ${OUT}
−
* View Annotated Screenshot:
 
−
<div class="mw-collapsible-content">
 
−
[[File:gcsnpcallOut.png|500px]]
 
−
</div>
 
−
</div>
 
  
−
Let's look at the vcfs directory:
+
'''WAIT!!!!! DO NOT RUN THIS COMMAND, because it will take ~50 minutes to finish'''.
−
ls ${OUT}/vcfs
 
−
Just a <code>chr22</code> directory, so look inside of there:
 
−
ls ${OUT}/vcfs/chr22
 
−
;Can you identify the final filtered VCF and the associated summary file?
 
−
<div class="mw-collapsible mw-collapsed" style="width:350px">
 
−
* Answer & annotated directory listing:
 
−
<div class="mw-collapsible-content">
 
−
<ul>
 
−
<li>Filtered VCF (SVM & hard filters): chr22.filtered.vcf.gz</li>
 
−
<li>Summary file: chr22.filtered.sites.vcf.summary</li>
 
−
</ul>
 
−
[[File:vcfsout.png|600px]]
 
−
</div>
 
−
</div>
 
  
−
Now, let's look in the split directory for the VCF with just the passing variants:
+
Instead, let's look what the output would have looked like.
−
ls ${OUT}/split/chr22
 
−
;Which file do you think is the one you want?
 
−
<div class="mw-collapsible mw-collapsed" style="width:350px">
 
−
* Answer:
 
−
<div class="mw-collapsible-content">
 
−
<ul>
 
−
<li>chr22.filtered.PASS.vcf.gz</li>
 
−
</ul>
 
−
[[File:splitOut.png|600px]]
 
−
</div>
 
−
</div>
 
  
−
=== Filtering Summary Statistics ===
+
ls ${SS}/svtoolkit/metadata
  
−
  cat ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf.summary
+
  cpt  depth  depth.dat  gcprofile  gcprofiles.zip  genome_sizes.txt  isd  isd.dist.bin  spans  spans.dat
  
−
<div class="mw-collapsible mw-collapsed" style="width:250px">
+
The directory contains metadata output and other intermediate files produced by "GenomeSTRiP SVProcess" step.  
−
View Screenshot
 
−
<div class="mw-collapsible-content">
 
−
[[File:filterSum.png]]
 
−
</div>
 
−
</div>
 
  
−
'''To understand how to interpret the filtering summary statistics, please refer to [[Understanding vcf-summary output]]'''
+
See [[http://gatkforums.broadinstitute.org/discussion/1514/svpreprocess-queue-script]] for the details of the Preprocess step.
  
−
=== Filtered VCF ===
+
NOTE: You don't always have to create the metadata on your own. You can in principle use the public metadata generated for 1000G samples, under the assumption that the metadata share similar characteristics to your samples. But if you have enough computing resources, the best practice is to create metadata specifically for your sequence data.
  
−
Let's look at the filtered sites file.
+
== Running GotCloud/GenomeSTRiP Discovery Pipeline ==
−
less -S ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf
 
  
−
* Scroll down until you find some variants.
+
To discover large deletions from the 62 BAMs we are using for this workshop, you can run the following command
−
** Use space bar to jump a full page
+
time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-discovery --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT}
−
** Use down arrow to move down one line
+
* <code>${SS}/svtoolkit/bin/genomestrip.pl -run-discovery</code> runs the GenomeSTRiP Discovery Pipeline
−
* Scroll right: lots of info fields, but no per sample genotype information
+
* <code>--metadata ${SS}/svtoolkit/metadata</code> points to the pre-made metadata file as explained in the previous section, [[#Running GotCloud/GenomeSTRiP Metadata Pipeline|Running GotCloud/GenomeSTRiP Metadata Pipeline]].
 +
* <code>--conf ${SS}/gotcloud.conf</code> points to the configuration file to use.
 +
** The configuration for this test was downloaded with the seqshop input files (same as other tutorials).
 +
* <code>--numjobs</code> tells how many jobs to run in parallel
 +
** Depends on your system
 +
* <code>--region 22:36000000-37000000</code>
 +
** The sample files are just a small region of chromosome 22, so to save time, we tell the pipeline to ignore the other regions
 +
* <code>--base_prefix</code> tells the pipeline the prefix to append to relative paths.
 +
** The Configuration file cannot read environment variables, so we need to tell it 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
 +
** Based on <code>gotcloud.conf</code>, the GenomeSTRiP output will go in <code>$(OUT_DIR)/sv</code>
  
−
;What is the first filtered out variant that you find & what filter did it fail?
+
This will take ~2-3 minutes to finish.
−
<div class="mw-collapsible mw-collapsed" style="width:250px">
 
−
* Answer:
 
−
<div class="mw-collapsible-content">
 
−
<ul>
 
−
<li>It failed SVM filter</li>
 
−
</ul>
 
−
[[File:SvmFilt.png|550px]]
 
−
</div>
 
−
</div>
 
−
 
 
−
Remember, use <code>'q'</code> to exit out of <code>less</code>
 
−
q
 
  
 +
Let's see the final outputs produced.
  
−
Now, let's look at the filtered file with genotypes.
+
  less ${OUT}/sv/discovery/discovery.vcf
−
  zless -S ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz
 
  
−
* Scroll down until you find some variants.
+
You will see output file that looks like this
−
** Use space bar to jump a full page
 
−
** Use down arrow to move down one line
 
−
* Scroll right until you should see per sample genotype information
 
  
−
<div class="mw-collapsible mw-collapsed" style="width:250px">
+
<div class="mw-collapsible mw-collapsed" style="width:400px">
−
* View annotated screenshot:
+
* Show Example
−
<div class="mw-collapsible-content">
+
<div class="mw-collapsible-content" style="width:800px">
−
[[File:SvmFiltGL.png|550px]]
+
[[File:Genomestrip discovery screenshot.png|800px]]
 
</div>
 
</div>
 
</div>
 
</div>
  
−
=== Passing SNPs ===
+
=== How many variants are filtered out? ===
  
−
Let's look at the file of just the pass sites:
+
Run the following command to see filtering statistics.
−
zless -S ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz
+
  grep -v ^# $OUT/sv/discovery/discovery.vcf | cut -f 7 | sort | uniq -c
  
−
* Scroll down: they all look like they <code>PASS</code>
+
<div class="mw-collapsible mw-collapsed" style="width:400px">
−
 
+
You will see the following output
−
Let's check if they are all PASS.
+
<div class="mw-collapsible-content" style="width:800px">
−
zcat ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
+
      7 COHERENCE;COVERAGE;DEPTH;DEPTHPVAL
−
It will return nothing since there are no non-passing variants in this file.
+
    17 COHERENCE;COVERAGE;DEPTH;DEPTHPVAL;PAIRSPERSAMPLE
−
<div class="mw-collapsible mw-collapsed" style="width:450px">
+
      3 COHERENCE;COVERAGE;DEPTH;PAIRSPERSAMPLE
−
;Want an explanation of this command?
+
      2 COHERENCE;COVERAGE;DEPTHPVAL;PAIRSPERSAMPLE
−
<div class="mw-collapsible-content">
+
      1 COHERENCE;COVERAGE;PAIRSPERSAMPLE
−
* zcat ...: uncompress the zipped VCF
+
      3 COVERAGE
−
* '|' : this takes the output of one command and sends it as input to the next
+
      1 COVERAGE;DEPTH
−
* grep -v "^#" : exclude any lines that start with "#" - headers
+
    67 COVERAGE;DEPTH;DEPTHPVAL
−
* cut -f 7 : extract the FILTER column (the 7th column)
+
    270 COVERAGE;DEPTH;DEPTHPVAL;PAIRSPERSAMPLE
−
* grep -v "PASS" : exclude any rows that have a "PASS" in the FILTER column
+
      2 COVERAGE;DEPTH;PAIRSPERSAMPLE
 +
      4 COVERAGE;DEPTHPVAL
 +
      5 COVERAGE;DEPTHPVAL;PAIRSPERSAMPLE
 +
      5 COVERAGE;PAIRSPERSAMPLE
 
</div>
 
</div>
 
</div>
 
</div>
  
−
Compare that to the filtered file we looked at before:
+
What does it mean? There is no "PASS filter" variants! This is because the metadata was created from only a small fraction of genome (with very unusual distribution of depth across chr22!). If whole-genome metadata was used, the results will look more reasonable, and you will have some "PASS" variants. Trust me!
−
zcat ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
+
 
−
;Do you see any filters?
+
=== What does each filter mean? ===
−
<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>
 
  
−
== GotCloud Genotype Refinement ==
+
Probably the most useful documentation of GenomeSTRiP is the powerpoint presentation available at http://www.broadinstitute.org/software/genomestrip/sites/default/files/materials/GATKWorkshop_GenomeSTRiP_tutorial_Dec2012.pdf
−
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.
+
In slide 27, you will see the following description of the filters
  
−
=== Genotype Refinement Input ===
+
[[File:Genomestrip filter description.png|600px]]
−
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.
+
== Running GotCloud/GenomeSTRiP Genotyping Pipeline ==
  
−
Note: the configuration file overrides the THUNDER command to make it go faster than the default settings so the tutorial will run faster:
+
The discovery pipeline only performs discovery of variant sites with filtering. You will need to iterate BAMs again to perform genotyping.
−
[[File:thunderConf.png|600px]]
+
* If running on a small machine, you may want to reduce <code>--numjobs</code> from 4 to 1.
 +
time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-genotype --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT} --gcroot ${GC}
 +
* The added <code>--gcroot ${GC}</code> option directs the pipeline to tabix/bgzip programs found within gotcloud.
  
−
=== Running GotCloud Genotype Refinement ===
+
This will take ~3 minutes to finish.
−
Since everything is setup, just run the following command (very similar to snpcall).
 
−
${GC}/gotcloud ldrefine --conf ${IN}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000
 
  
−
* Beagle will take about 2-3 minutes to complete
+
You can check the output by running
−
* Thunder will automatically run and will take another 3-4 minutes
+
zless $OUT/sv/genotype/genotype.vcf.gz
  
−
=== Genotype Refinement Output ===
+
You will see output similar to this
  
−
; What's new in the output directory?
+
<div class="mw-collapsible mw-collapsed" style="width:600px">
−
<ul>
+
You will see the output with genotype information
−
<div class="mw-collapsible mw-collapsed" style="width:500px">
+
<div class="mw-collapsible-content" style="width:800px">
−
<li>Answer</li>
+
[[File:Genomestrip genotype screenshot.png|800px]]
−
<div class="mw-collapsible-content">
 
−
<ul>
 
−
<li><code>beagle</code> directory : Beagle output</li>
 
−
<li><code>thunder</code> directory : Thunder output</li>
 
−
<li><code>umake.beagle.conf</code> : Configuration values used for GotCloud beagle</li>
 
−
<li><code>umake.beagle.Makefile</code> : GNU makefile for commands run as part of GotCloud beagle</li>
 
−
<li><code>umake.beagle.Makefile.log</code> : Log of the GotCloud beagle run</li>
 
−
<li><code>umake.thunder.*</code> files : Contain the configuration & steps used in GotCloud thunder</li>
 
−
</ul>
 
 
</div>
 
</div>
 
</div>
 
</div>
−
</ul>
 
  
−
Let's take a look at that interesting location we found in the [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Accessing_BAMs_by_Position|alignment tutorial]] : chromosome 22, positions 36907000-36907100
+
== Running GotCloud/GenomeSTRiP 3rd-party Site Genotyping/Filtering Pipeline ==
  
−
Use tabix to extract that from the VCFs:
+
You can take a 3rd-party site and genotype with GenomeSTRiP. Here we take a 1000 Genomes phase 1 sites and genotype them.
−
  ${GC}/bin/tabix ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
+
* If running on a small machine, you may want to reduce <code>--numjobs</code> from 4 to 1.
 +
  time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-thirdparty --in-vcf ${SS}/ext/1kg.phase1.chr22.36Mb.sites.vcf --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT} --gcroot ${GC} --numjobs 4
  
−
Remember, type 'q' to quit less.
+
This will take ~1 minute to finish.
−
q
 
  
−
;Are there any variants in this region?
+
You can also check the output by running
−
<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?
+
  zless $OUT/sv/thirdparty/genotype.vcf.gz
−
#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:600px">
−
<div class="mw-collapsible mw-collapsed" style="width:500px">
+
You will see the output with genotype information
−
<li>Answer:</li>
+
<div class="mw-collapsible-content" style="width:800px">
−
<div class="mw-collapsible-content">
+
[[File:Genomestrip thirdparty screenshot.png|800px]]
−
<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>
 
</div>
 
</div>
−
</ul>
 
  
−
Remember, type 'q' to quit less.
+
== What does a real SV look like? ==
−
q
 
  
−
=== Did I find interesting variants? ===
+
samtools tview does not provide a good way to visualize structural variants due to limited resolution to show large-scale 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.
+
IGV provides a good alternative way to visualize structural variants as shown in the xample below.
  
−
$GC/bin/tabix $OUT/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1
+
<div class="mw-collapsible mw-collapsed" style="width:600px">
−
 
+
Do you understand why this is a likely SV?
−
Did you see a variant at the position?
+
<div class="mw-collapsible-content" style="width:800px">
 +
[[File:Igvsvexample.png|800px]]
 +
</div>
 +
</div>

Latest revision as of 14:37, 13 November 2014

Note: the latest version of this practical is available at: SeqShop: Analysis of Structural Variation Practical

  • The ones here is the original one from the June workshop (updated to be run from elsewhere)


Goals of This Session

  • What we want to learn is calling large deletions using GenomeSTRiP implemented in GotCloud pipeline
    • How to prepare metadata for running GenomeSTRiP.
    • How to perform variant discovery and filtering for large deletions
    • How to perform genotyping for large deletions
    • How to perform variant discovery and filtering from third party sites.

Please refer to Lecture slides for more general background.

GenomeSTRiP

GenomeSTRiP was developed at the Broad Institute and at the McCarroll Lab at the Harvard Medical School Department of Genetics: http://www.broadinstitute.org/software/genomestrip/

If you use GenomeSTRiP for your research, please cite it: 

Handsaker RE, Korn JM, Nemesh J, McCarroll SA
Discovery and genotyping of genome structural polymorphism by sequencing on a population scale.
Nature genetics 43, 269-276 (2011)
PMID: 21317889

GenomeStrip is currently included in with the seqshop example data under the svtoolkit directory. We have added the bin/ sub-directory to add a high level pipeline that will run genomestrip in the same framework as GotCloud.

Setup in person at the SeqShop Workshop

This section is specifically for the SeqShop Workshop computers.

If you are not running during the SeqShop Workshop, please skip this section.

Login to the seqshop-server Linux Machine

This section will appear redundantly in each session. If you are already logged in or know how to log in to the server, please skip this section

  1. Login to the windows machine
    • The username/password for the Windows machine should be written on the right-hand monitor
  2. Start xming so you can open external windows on our Linux machine
    • Start->Enter "Xming" in the search and select "Xming" from the program list
    • Nothing will happen, but Xming was started.
    • View Screenshot

      • Xming.png

  3. Open putty
    • Start->Enter "putty" in the search and select "PuTTY" from the program list
    • View Screenshot

      • PuttyS.png

  4. Configure PuTTY in the PuTTY Configuration window
    • Host Name: seqshop-server.sph.umich.edu
    • View Screenshot

      • Seqshop.png

    • Setup to allow you to open external windows:
      • In the left pannel: Connection->SSH->X11
        • Add a check mark in the box next to Enable X11 forwarding
        • View Screenshot

          • SeqshopX11.png

    • Click Open
    • If it prompts about a key, click OK
  5. Enter your provided username & password as provided


You should now be logged into a terminal on the seqshop-server and be able to access the test files.

  • If you need another terminal, repeat from step 3.

Login to the seqshop Machine

So you can each run multiple jobs at once, we will have you run on 4 different machines within our seqshop setup.

  • You can only access these machines after logging onto seqshop-server

3 users logon to: 

ssh -X seqshop1

3 users logon to: 

ssh -X seqshop2

2 users logon to: 

ssh -X seqshop3

2 users logon to: 

ssh -X seqshop4

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 will setup some environment variables to point you to

  • GotCloud program
  • Tutorial input files
  • Setup an output directory
    • It will leave your output directory from the previous tutorial in tact.
source /home/mktrost/seqshop/setup.txt
  • You won't see any output after running source
    • 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.

View setup.txt

Setup.png


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.

If you are running during the SeqShop Workshop, please skip this section.

This tutorial builds on the alignment tutorial, if you have not already, please first run that tutorial: Alignment Tutorial

It also uses the bam.index file created in the SnpCall Tutorial. If you have not yet run that tutorial, please follow the directions at: GotCloud BAM Index File


Setup your run environment

Environment variables will be used throughout the tutorial.

We recommend that you setup these variables so you won't have to modify every command in the tutorial.

  1. Point to where you installed GotCloud
  2. Point to where you installed the seqshop files
  3. Point to where you want the output to go
Using bash (replace the paths below with the appropriate paths):
export GC=~/seqshop/gotcloud
export SS=~/seqshop/example
export OUT=~/seqshop/output
Using tcsh (replace the paths below with the appropriate paths):
setenv GC ~/seqshop/gotcloud
setenv SS ~/seqshop/example
setenv OUT ~/seqshop/output

Examining GotCloud/GenomeSTRiP Input files

Sequnce Alignment Files: BAM Files and Index Files

The GotCloud GenomeSTRiP structural variant caller takes the same inputs as GotCloud snpcall.

  • BAMs->SVs rather than BAMs->SNPs

If you want a reminder, of what they look like, here is a link to the previous tutorial : GotCloud SnpCall Input Files

If you want to check if you still have the bam index file, run 

head ${OUT}/bam.index
  • View Results
    • HG00641	ALL	/net/seqshop-server/hmkang/out/bams/HG00641.recal.bam
HG00640	ALL	/net/seqshop-server/hmkang/out/bams/HG00640.recal.bam
HG00551	ALL	/net/seqshop-server/hmkang/out/bams/HG00551.recal.bam
HG00553	ALL	/net/seqshop-server/hmkang/out/bams/HG00553.recal.bam
HG00554	ALL	bams/HG00554.recal.bam
HG00637	ALL	bams/HG00637.recal.bam
HG00638	ALL	bams/HG00638.recal.bam
HG00734	ALL	bams/HG00734.recal.bam
HG00736	ALL	bams/HG00736.recal.bam
HG00737	ALL	bams/HG00737.recal.bam

Also, make sure that you have only 62 samples (you did not append new files twice) 

wc -l ${OUT}/bam.index

Your expected output is similar to this. 

62 /net/seqshop-server/hmkang/out/bam.index

Reference Files

Reference files can be downloaded with GotCloud or from other sources.

Similar to SNP and Indel calling, you need

  1. Reference genome FASTA file

For running GenomeSTRiP, you additionally need:

  1. Masked FASTA file to exclude hard-to-align regions
  2. PloidyMap file indicating the regions of genomes with unusual ploidy (e.g. chrX, chrY)

We looked at them in previous tutorials, but you can take another look at the chromosome 22 reference files included for this tutorial: 

ls ${SS}/ref22
  • View Results

    • RefDir.png


Additional Reference files required just for Structural Variation: 

ls ${SS}/svtoolkit/ref
  • View Results
    • human_g1k_v37.chr22.mask.100.fasta  human_g1k_v37.chr22.mask.100.fasta.dict  human_g1k_v37.chr22.mask.100.fasta.fai


Parameters files required just for Structural Variation: 

ls ${SS}/svtoolkit/conf
  • View Results
    • genstrip_parameters.txt  humgen_g1k_v37_ploidy.chr22.map  humgen_g1k_v37_ploidy.map

GotCloud Configuration File

We will use the same configuration file we used for the GotCloud Align tutorial.

See SeqShop: Alignment: GotCloud Configuration File for more details

  • Note we want to limit snpcall to just chr22 so the configuration already has CHRS = 22 (default was 1-22 & X).

For more information on configuration, see: GotCloud snpcall: Configuration File

Check out the GenomeStrip specific settings at the end of the configuration file 

tail -n 8 ${SS}/gotcloud.conf
  • View Results
    • ##############################
## GenomeSTRIP
#############################
GENOMESTRIP_OUT = $(OUT_DIR)/sv
GENOMESTRIP_SVTOOLKIT_DIR = svtoolkit
GENOMESTRIP_MASK_FASTA = $(GENOMESTRIP_SVTOOLKIT_DIR)/ref/human_g1k_v37.chr22.mask.100.fasta
GENOMESTRIP_PLOIDY_MAP = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/humgen_g1k_v37_ploidy.chr22.map
GENOMESTRIP_PARAM = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/genstrip_parameters.txt

Before starting... a few 'why' questions..

Why use GenomeSTRiP?

  1. GenomeSTRiP is a mature software for detecting and genotyping large deletions (and duplications soon to be implemented). In 1000 Genomes, GenomeSTRiP was demonstrated as one of the top-performing SV caller in most evaluation metrics.
  2. GenomeSTRiP is a great tool to integrate across multiple structural variant calls. When multiple structural variant calls exists, all the other variants can be genotyped and filtered with GenomeSTRiP, and that is how 1000 Genomes structural variant call sets were made.
  3. Currently, GenomeSTRiP only allows calling large deletions, but duplicate calling pipeline is under way.

Why do we use GotCloud/GenomeSTRiP pipeline?

  1. The main purpose of GotCloud pipelines is to provide a pipeline for users with limited knowledge and experience with high performance computing environment.
    • GotCloud/GenomeSTRiP provide a simple interface consistent to alignment, SNP, and indel calling.
    • GenomeSTRiP itself also provides a straightforward pipeline to use as standalone software
  2. GotCloud supports a variety of cluster environment that is not currently supported by GenomeSTRiP
    • GenomeSTRiP is designed based on a framework called Qscript, which provide a nice support for LSF cluster system
    • GotCloud support many additional cluster environments such as MOSIX or SLURM we use locally at Michigan.
  3. GotCloud also provide a fault-tolerant solution for large-scale jobs.
    • GotCloud automatically picks up jobs from the point where it failed. This allows easier and simpler run against potential technical glitches in the system.

Overview of GotCloud/GenomeSTRiP pipeline

GotCloud/GenomeSTRiP pipeline consists of three separate steps.

  • Preprocess step : Create metadata summarizing the GC profiles, depth distribution, insert size distribution for accurate discovery and genotyping of structural variants.
  • Discovery step : Perform variant discovery split by region, across all samples. Also, perform variant filtering based on expert knowledge.
  • Genotyping step : Iterate discovered variants across the samples and calculate the genotype likelihood of for each possible genotype.

In addition, if one wants to genotype structural variants from other structural variant caller, there is a step available.

  • Third-party Genotyping and Filtering step : Perform genotyping on the variant sites specified by an input VCF, and also perform variant filtering.

Running GotCloud/GenomeSTRiP Metadata Pipeline

We first need to create metadata summarizing genome-wide statistics such as GC profiles, depth distribution, insert size distributions.

In principle, the metadata can be created from the input BAM files by running the following command 

perl ${SS}/svtoolkit/bin/genomestrip.pl -run-metadata --conf ${SS}/gotcloud.conf --numjobs 2 --base-prefix ${SS} --outdir ${OUT}

WAIT!!!!! DO NOT RUN THIS COMMAND, because it will take ~50 minutes to finish.

Instead, let's look what the output would have looked like. 

ls ${SS}/svtoolkit/metadata

cpt  depth  depth.dat  gcprofile  gcprofiles.zip  genome_sizes.txt  isd  isd.dist.bin  spans  spans.dat

The directory contains metadata output and other intermediate files produced by "GenomeSTRiP SVProcess" step.

See [[1]] for the details of the Preprocess step.

NOTE: You don't always have to create the metadata on your own. You can in principle use the public metadata generated for 1000G samples, under the assumption that the metadata share similar characteristics to your samples. But if you have enough computing resources, the best practice is to create metadata specifically for your sequence data.

Running GotCloud/GenomeSTRiP Discovery Pipeline

To discover large deletions from the 62 BAMs we are using for this workshop, you can run the following command 

time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-discovery --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT}
  • ${SS}/svtoolkit/bin/genomestrip.pl -run-discovery runs the GenomeSTRiP Discovery Pipeline
  • --metadata ${SS}/svtoolkit/metadata points to the pre-made metadata file as explained in the previous section, Running GotCloud/GenomeSTRiP Metadata Pipeline.
  • --conf ${SS}/gotcloud.conf points to the configuration file to use.
    • The configuration for this test was downloaded with the seqshop input files (same as other tutorials).
  • --numjobs tells how many jobs to run in parallel
    • Depends on your system
  • --region 22:36000000-37000000
    • The sample files are just a small region of chromosome 22, so to save time, we tell the pipeline to ignore the other regions
  • --base_prefix tells the pipeline the prefix to append to relative paths.
    • The Configuration file cannot read environment variables, so we need to tell it the path to the input files, ${SS}
    • Alternatively, gotcloud.conf could be updated to specify the full paths
  • --out_dir 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
    • Based on gotcloud.conf, the GenomeSTRiP output will go in $(OUT_DIR)/sv

This will take ~2-3 minutes to finish.

Let's see the final outputs produced. 

less ${OUT}/sv/discovery/discovery.vcf

You will see output file that looks like this

  • Show Example

Genomestrip discovery screenshot.png

How many variants are filtered out?

Run the following command to see filtering statistics. 

 grep -v ^# $OUT/sv/discovery/discovery.vcf | cut -f 7 | sort | uniq -c

You will see the following output 

     7 COHERENCE;COVERAGE;DEPTH;DEPTHPVAL
    17 COHERENCE;COVERAGE;DEPTH;DEPTHPVAL;PAIRSPERSAMPLE
     3 COHERENCE;COVERAGE;DEPTH;PAIRSPERSAMPLE
     2 COHERENCE;COVERAGE;DEPTHPVAL;PAIRSPERSAMPLE
     1 COHERENCE;COVERAGE;PAIRSPERSAMPLE
     3 COVERAGE
     1 COVERAGE;DEPTH
    67 COVERAGE;DEPTH;DEPTHPVAL
   270 COVERAGE;DEPTH;DEPTHPVAL;PAIRSPERSAMPLE
     2 COVERAGE;DEPTH;PAIRSPERSAMPLE
     4 COVERAGE;DEPTHPVAL
     5 COVERAGE;DEPTHPVAL;PAIRSPERSAMPLE
     5 COVERAGE;PAIRSPERSAMPLE

What does it mean? There is no "PASS filter" variants! This is because the metadata was created from only a small fraction of genome (with very unusual distribution of depth across chr22!). If whole-genome metadata was used, the results will look more reasonable, and you will have some "PASS" variants. Trust me!

What does each filter mean?

Probably the most useful documentation of GenomeSTRiP is the powerpoint presentation available at http://www.broadinstitute.org/software/genomestrip/sites/default/files/materials/GATKWorkshop_GenomeSTRiP_tutorial_Dec2012.pdf

In slide 27, you will see the following description of the filters 

Genomestrip filter description.png

Running GotCloud/GenomeSTRiP Genotyping Pipeline

The discovery pipeline only performs discovery of variant sites with filtering. You will need to iterate BAMs again to perform genotyping.

  • If running on a small machine, you may want to reduce --numjobs from 4 to 1.
time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-genotype --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT} --gcroot ${GC}
  • The added --gcroot ${GC} option directs the pipeline to tabix/bgzip programs found within gotcloud.

This will take ~3 minutes to finish.

You can check the output by running 

zless $OUT/sv/genotype/genotype.vcf.gz

You will see output similar to this

You will see the output with genotype information 

Genomestrip genotype screenshot.png

Running GotCloud/GenomeSTRiP 3rd-party Site Genotyping/Filtering Pipeline

You can take a 3rd-party site and genotype with GenomeSTRiP. Here we take a 1000 Genomes phase 1 sites and genotype them.

  • If running on a small machine, you may want to reduce --numjobs from 4 to 1.
time perl ${SS}/svtoolkit/bin/genomestrip.pl -run-thirdparty --in-vcf ${SS}/ext/1kg.phase1.chr22.36Mb.sites.vcf --metadata ${SS}/svtoolkit/metadata --conf ${SS}/gotcloud.conf --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT} --gcroot ${GC} --numjobs 4

This will take ~1 minute to finish.

You can also check the output by running 

zless $OUT/sv/thirdparty/genotype.vcf.gz

You will see the output with genotype information

Genomestrip thirdparty screenshot.png

What does a real SV look like?

samtools tview does not provide a good way to visualize structural variants due to limited resolution to show large-scale variants.

IGV provides a good alternative way to visualize structural variants as shown in the xample below.

Do you understand why this is a likely SV?

Igvsvexample.png

Retrieved from "http://genome.sph.umich.edu/w/index.php?title=SeqShop:_Analysis_of_Structural_Variation_Practical,_June_2014&oldid=11992"