View source for SeqShop: Analysis of Structural Variation Practical, June 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 [[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}}

=== 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 <code>source</code>
** 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">
View setup.txt
<div class="mw-collapsible-content">
[[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>

== 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 : [[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)

 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.
* For this practical, I already downloaded them for you.
* See [[GotCloud: Genetic Reference and Resource Files]] for more information on downloading/generating reference files

Similar to SNP and Indel calling, you need
# Reference genome FASTA file

For running GenomeSTRiP, you additionally need:
# Masked FASTA file to exclude hard-to-align regions
# 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

<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Results</li>
<div class="mw-collapsible-content" style="width:800px">
[[File:RefDir.png|500px]]
</div>
</div>
</ul>


Additional Reference files required just for Structural Variation:
 ls ${SS}/svtoolkit/ref
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Results</li>
<div class="mw-collapsible-content" style="width:800px">
 human_g1k_v37.chr22.mask.100.fasta  human_g1k_v37.chr22.mask.100.fasta.dict  human_g1k_v37.chr22.mask.100.fasta.fai
</div>
</div>
</ul>


Parameters files required just for Structural Variation:
 ls ${SS}/svtoolkit/conf

<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Results</li>
<div class="mw-collapsible-content" style="width:800px">
 genstrip_parameters.txt  humgen_g1k_v37_ploidy.chr22.map  humgen_g1k_v37_ploidy.map
</div>
</div>
</ul>

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

See [[SeqShop:_Sequence_Mapping_and_Assembly_Practica, June 2014l#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]]

Check out the GenomeStrip specific settings at the end of the configuration file
 tail -n 8 ${SS}/gotcloud.conf

<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Results</li>
<div class="mw-collapsible-content" style="width:800px">
 ##############################
 ## 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
</div>
</div>
</ul>

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

=== Why use GenomeSTRiP?===
# 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.

=== Why do we use GotCloud/GenomeSTRiP pipeline? ===
# 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.

== 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 [[http://gatkforums.broadinstitute.org/discussion/1514/svpreprocess-queue-script]] 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}
* <code>${SS}/svtoolkit/bin/genomestrip.pl -run-discovery</code> runs the GenomeSTRiP Discovery Pipeline
* <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>

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

<div class="mw-collapsible mw-collapsed" style="width:400px">
* Show Example
<div class="mw-collapsible-content" style="width:800px">
[[File:Genomestrip discovery screenshot.png|800px]]
</div>
</div>

=== 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

<div class="mw-collapsible mw-collapsed" style="width:400px">
You will see the following output
<div class="mw-collapsible-content" style="width:800px">
      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
</div>
</div>

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 

 [[File:Genomestrip filter description.png|600px]]

== 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 <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.

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

<div class="mw-collapsible mw-collapsed" style="width:600px">
You will see the output with genotype information
<div class="mw-collapsible-content" style="width:800px">
 [[File:Genomestrip genotype screenshot.png|800px]]
</div>
</div>

== 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 <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

This will take ~1 minute to finish.

You can also check the output by running

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

<div class="mw-collapsible mw-collapsed" style="width:600px">
You will see the output with genotype information
<div class="mw-collapsible-content" style="width:800px">
[[File:Genomestrip thirdparty screenshot.png|800px]]
</div>
</div>

== 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.

<div class="mw-collapsible mw-collapsed" style="width:600px">
Do you understand why this is a likely SV?
<div class="mw-collapsible-content" style="width:800px">
[[File:Igvsvexample.png|800px]]
</div>
</div>