Difference between revisions of "SeqShop: Analysis of Structural Variation Practical, June 2014"
(39 intermediate revisions by 2 users not shown) | |||
Line 1: | Line 1: | ||
− | '' | + | '''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 == | == Goals of This Session == | ||
Line 7: | 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 | ||
− | * GotCloud program | + | This will setup some environment variables to point you to |
+ | * [[GotCloud]] program | ||
* Tutorial input files | * Tutorial input files | ||
* Setup an output directory | * Setup an output directory | ||
** It will leave your output directory from the previous tutorial in tact. | ** It will leave your output directory from the previous tutorial in tact. | ||
− | source /home/ | + | source /home/mktrost/seqshop/setup.txt |
* 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 | + | ** If you want to view the detail of the setup, type |
− | less /home/ | + | 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]] | |
− | + | </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]] | ||
+ | |||
+ | 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 == | == Examining GotCloud/GenomeSTRiP Input files == | ||
Line 44: | Line 75: | ||
=== Sequnce Alignment Files: BAM Files and Index Files=== | === Sequnce Alignment Files: BAM Files and Index Files=== | ||
− | The GotCloud | + | The GotCloud GenomeSTRiP structural variant caller takes the same inputs as GotCloud snpcall. |
* BAMs->SVs rather than BAMs->SNPs | * 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#Examining_GotCloud_SnpCall_Input_files|GotCloud SnpCall Input Files]] | + | 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 | If you want to check if you still have the bam index file, run | ||
− | head $OUT/bam.index | + | head ${OUT}/bam.index |
<ul> | <ul> | ||
Line 73: | Line 104: | ||
Also, make sure that you have only 62 samples (you did not append new files twice) | Also, make sure that you have only 62 samples (you did not append new files twice) | ||
− | wc -l $OUT/bam.index | + | wc -l ${OUT}/bam.index |
+ | |||
+ | Your expected output is similar to this. | ||
62 /net/seqshop-server/hmkang/out/bam.index | 62 /net/seqshop-server/hmkang/out/bam.index | ||
− | |||
=== Reference Files === | === Reference Files === | ||
Line 90: | Line 122: | ||
# 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 | + | We looked at them in previous tutorials, but you can take another look at the chromosome 22 reference files included for this tutorial: |
− | ls $ | + | ls ${SS}/ref22 |
<ul> | <ul> | ||
Line 98: | Line 130: | ||
<li>View Results</li> | <li>View Results</li> | ||
<div class="mw-collapsible-content" style="width:800px"> | <div class="mw-collapsible-content" style="width:800px"> | ||
− | + | [[File:RefDir.png|500px]] | |
− | |||
− | |||
− | |||
− | |||
− | |||
</div> | </div> | ||
</div> | </div> | ||
Line 109: | Line 136: | ||
− | Additional | + | Additional Reference files required just for Structural Variation: |
− | + | ls ${SS}/svtoolkit/ref | |
− | ls $ | ||
− | |||
<ul> | <ul> | ||
<div class="mw-collapsible mw-collapsed" style="width:200px"> | <div class="mw-collapsible mw-collapsed" style="width:200px"> | ||
Line 123: | Line 148: | ||
− | ls $ | + | Parameters files required just for Structural Variation: |
+ | ls ${SS}/svtoolkit/conf | ||
<ul> | <ul> | ||
Line 135: | Line 161: | ||
=== GotCloud Configuration File === | === GotCloud Configuration File === | ||
− | We will use | + | We will use the same configuration file we used for the GotCloud Align tutorial. |
− | See [[SeqShop: | + | 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). | * 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]] | 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> | <ul> | ||
Line 151: | Line 175: | ||
<li>View Results</li> | <li>View Results</li> | ||
<div class="mw-collapsible-content" style="width:800px"> | <div class="mw-collapsible-content" style="width:800px"> | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
############################## | ############################## | ||
## GenomeSTRIP | ## 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_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_PLOIDY_MAP = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/humgen_g1k_v37_ploidy.chr22.map | ||
Line 196: | Line 190: | ||
=== Why use GenomeSTRiP?=== | === Why use GenomeSTRiP?=== | ||
− | # GenomeSTRiP is a mature software for detecting and genotyping large deletions (and duplications soon to be implemented). In 1000 Genomes, 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. | # 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 | + | === 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. | + | # 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 | # 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 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 == | == Overview of GotCloud/GenomeSTRiP pipeline == | ||
Line 223: | Line 219: | ||
In principle, the metadata can be created from the input BAM files by running the following command | 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 $ | + | ls ${SS}/svtoolkit/metadata |
cpt depth depth.dat gcprofile gcprofiles.zip genome_sizes.txt isd isd.dist.bin spans spans.dat | cpt depth depth.dat gcprofile gcprofiles.zip genome_sizes.txt isd isd.dist.bin spans spans.dat | ||
Line 235: | Line 233: | ||
See [[http://gatkforums.broadinstitute.org/discussion/1514/svpreprocess-queue-script]] for the details of the Preprocess 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 | + | 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 == | == 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 | 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. | |
− | |||
− | This will take ~2 minutes to finish. | ||
Let's see the final outputs produced. | Let's see the final outputs produced. | ||
− | less $OUT/sv/discovery/discovery.vcf | + | less ${OUT}/sv/discovery/discovery.vcf |
You will see output file that looks like this | 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]] | [[File:Genomestrip discovery screenshot.png|800px]] | ||
+ | </div> | ||
+ | </div> | ||
=== How many variants are filtered out? === | === How many variants are filtered out? === | ||
Line 286: | Line 302: | ||
[[File:Genomestrip filter description.png|600px]] | [[File:Genomestrip filter description.png|600px]] | ||
− | |||
== Running GotCloud/GenomeSTRiP Genotyping Pipeline == | == 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. | 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 $ | + | 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. | This will take ~3 minutes to finish. | ||
Line 306: | Line 322: | ||
[[File:Genomestrip genotype screenshot.png|800px]] | [[File:Genomestrip genotype screenshot.png|800px]] | ||
</div> | </div> | ||
− | <div> | + | </div> |
== Running GotCloud/GenomeSTRiP 3rd-party Site Genotyping/Filtering Pipeline == | == 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. | 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. | |
− | |||
− | This will take ~ | ||
You can also check the output by running | You can also check the output by running | ||
Line 327: | Line 343: | ||
</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"> | <div class="mw-collapsible mw-collapsed" style="width:600px"> | ||
Line 338: | Line 355: | ||
</div> | </div> | ||
</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
- Login to the windows machine
- The username/password for the Windows machine should be written on the right-hand monitor
- 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.
- Open putty
- Start->Enter "putty" in the search and select "PuTTY" from the program list
- Configure PuTTY in the PuTTY Configuration window
- Host Name:
seqshop-server.sph.umich.edu
- 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
- Click
Open
- If it prompts about a key, click
OK
- 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.
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.
- Point to where you installed GotCloud
- Point to where you installed the seqshop files
- 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.
- 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
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?
- 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 [[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
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
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
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
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.