Changes

From Genome Analysis Wiki
Jump to navigationJump to search

SeqShop: Analysis of Structural Variation Practical, June 2014 (view source)

Revision as of 14:37, 13 November 2014

2,078 bytes added ,  14:37, 13 November 2014
no edit summary
Line 1: Line 1: −
''(WARNING: Under Construction)''
+
'''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 a '''slightly different''' setup from what 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
 
** 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/hmkang/seqshop/setup.txt
+
  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 set up, type
+
** If you want to view the detail of the setup, type
−
  less /home/hmkang/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" style="width:800px">
+
<div class="mw-collapsible-content">
−
export GC=/home/hmkang/seqshop/gotcloud
+
[[File:setup.png|500px]]
−
export IN=/home/hmkang/seqshop/inputs
+
</div>
−
export REF=/home/hmkang/seqshop/reference/chr22
+
</div>
−
export VTREF=/home/hmkang/seqshop/reference/vtRef
  −
export SV=/home/hmkang/seqshop/reference/svtoolkit
  −
export EXT=/home/hmkang/seqshop/external
  −
export OUT=~/out
  −
mkdir -p ${OUT}
   
</div>
 
</div>
 
</div>
 
</div>
   −
Do you notice what the differences were?
+
 
 +
== 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 Indel caller takes the same inputs as GotCloud snpcall.
+
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.
 
Your expected output is similar to this.
Line 91: 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 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>
Line 99: 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">
−
1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz            hapmap_3.3.b37.sites.chr22.vcf.gz.tbi  human.g1k.v37.chr22.fa.bwt
+
[[File:RefDir.png|500px]]
−
1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz.tbi        human.g1k.v37.chr22-bs.umfa            human.g1k.v37.chr22.fa.fai
  −
1kg.pilot_release.merged.indels.sites.hg19.chr22.vcf  human.g1k.v37.chr22.dict              human.g1k.v37.chr22.fa.pac
  −
dbsnp_135.b37.chr22.vcf.gz                            human.g1k.v37.chr22.fa                human.g1k.v37.chr22.fa.sa
  −
dbsnp_135.b37.chr22.vcf.gz.tbi                        human.g1k.v37.chr22.fa.amb            human.g1k.v37.chr22.winsize100.gc
  −
hapmap_3.3.b37.sites.chr22.vcf.gz                    human.g1k.v37.chr22.fa.ann
   
</div>
 
</div>
 
</div>
 
</div>
Line 110: Line 136:       −
Additional reference and parameters
+
Additional Reference files required just for Structural Variation:
−
 
+
  ls ${SS}/svtoolkit/ref
−
  ls $SV/ref
  −
 
   
<ul>
 
<ul>
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
 
<div class="mw-collapsible mw-collapsed" style="width:200px">
Line 124: Line 148:       −
  ls $SV/conf  
+
Parameters files required just for Structural Variation:
 +
  ls ${SS}/svtoolkit/conf
    
<ul>
 
<ul>
Line 136: Line 161:     
=== GotCloud Configuration File ===
 
=== GotCloud Configuration File ===
−
We will use a slightly modified version of configuration file as we used yesterday in GotCloud Align.
+
We will use the same configuration file we used for the GotCloud Align tutorial.
   −
See [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#GotCloud Configuration File|SeqShop: Alignment: GotCloud Configuration File]] for more details
+
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]]
−
* Contains information on how to configure for exome/targeted sequencing
  −
  −
Check out what was changed.
     −
  cat $IN/gotcloud.conf
+
Check out the GenomeStrip specific settings at the end of the configuration file
 +
  tail -n 8 ${SS}/gotcloud.conf
    
<ul>
 
<ul>
Line 152: 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">
−
IN_DIR = $(GOTCLOUD_ROOT)/../inputs
  −
  −
INDEX_FILE = $(IN_DIR)/align.index
  −
FASTQ_PREFIX = $(IN_DIR)/fastq
  −
BAM_PREFIX = $(IN_DIR)/
  −
  −
OUT_DIR = out
  −
BAM_INDEX =  $(OUT_DIR)/bam.index
  −
  −
############
  −
# References
  −
REF_DIR = $(GOTCLOUD_ROOT)/../reference/chr22
  −
AS = NCBI37  # Genome assembly identifier
  −
REF = $(REF_DIR)/human.g1k.v37.chr22.fa
  −
DBSNP_VCF =  $(REF_DIR)/dbsnp_135.b37.chr22.vcf.gz
  −
HM3_VCF = $(REF_DIR)/hapmap_3.3.b37.sites.chr22.vcf.gz
  −
INDEL_PREFIX = $(REF_DIR)/1kg.pilot_release.merged.indels.sites.hg19
  −
OMNI_VCF = $(REF_DIR)/1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz
  −
  −
MAP_TYPE = BWA_MEM
  −
  −
###############
  −
CHRS = 22
  −
  −
######### THUNDER ########
  −
# Update so it will run faster for the tutorial
  −
#  * 10 rounds instead of 30 (-r 10)
  −
#  * without --compact option
  −
#  Runs faster, but uses more memory, but not a lot for the small example
  −
THUNDER = $(BIN_DIR)/thunderVCF -r 10 --phase --dosage --inputPhased $(THUNDER_STATES)
  −
   
  ##############################
 
  ##############################
 
  ## GenomeSTRIP
 
  ## GenomeSTRIP
 
  #############################
 
  #############################
−
  GENOMESTRIP_SVTOOLKIT_DIR = $(GOTCLOUD_ROOT)/../reference/svtoolkit
+
  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 201: Line 194:  
# Currently, GenomeSTRiP only allows calling large deletions, but duplicate calling pipeline is under way.
 
# Currently, GenomeSTRiP only allows calling large deletions, but duplicate calling pipeline is under way.
   −
=== Why do we use GotCloud/GenomeSTRiP pipeline instead of directly using GenomeSTRiP itself? ===
+
=== 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.
−
#* Although GenomeSTRiP provides a reasonably straightforward pipeline, it still requires a detailed understanding of GATK framework and the details of parameter.
+
#* GotCloud/GenomeSTRiP provide a simple interface consistent to alignment, SNP, and indel calling.
−
#* GotCloud aims to provide more simpler way to run these procedure
+
#* 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, but it does not support many other cluster enviroments such as MOSIX or SLURM we use locally.
+
#* 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.
 
#* GotCloud automatically picks up jobs from the point where it failed. This allows easier and simpler run against potential technical glitches in the system.
Line 225: 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
−
  time perl $GC/bin/genomestrip.pl -run-metadata --out $OUT/sv --conf $IN/gotcloud.conf --numjobs 2
+
  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'''.
 
'''WAIT!!!!! DO NOT RUN THIS COMMAND, because it will take ~50 minutes to finish'''.
Line 231: Line 225:  
Instead, let's look what the output would have looked like.
 
Instead, let's look what the output would have looked like.
   −
  ls $IN/metadata
+
  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 239: 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 daat.
+
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>
   −
time perl $GC/bin/genomestrip.pl -run-discovery --metadata $IN/metadata --out $OUT/sv --conf $IN/gotcloud.conf --region 22:36000000-37000000 --numjobs 2
+
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
Line 299: Line 306:     
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 $GC/bin/genomestrip.pl -run-genotype --metadata $IN/metadata --out $OUT/sv --conf $IN/gotcloud.conf --region 22:36000000-37000000 --numjobs 4
+
  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 319: Line 327:     
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 $GC/bin/genomestrip.pl -run-thirdparty --in-vcf $EXT/1kg.phase1.chr22.36Mb.sites.vcf --metadata $IN/metadata --out $OUT/sv --conf $IN/gotcloud.conf --region 22:36000000-37000000 --numjobs 4
+
  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 ~1 minute to finish.
Retrieved from "http://genome.sph.umich.edu/wiki/Special:MobileDiff/10925...11992"

Navigation menu