SeqShop: Analysis of Structural Variation Practical, December 2014

Introduction

Main Workshop wiki page: SeqShop: December 2014

See the introductory slides for an intro to this tutorial.

Goals of This Session

• What we want to learn
  • 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.

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

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

Open putty

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

Configure PuTTY in the PuTTY Configuration window

• Host Name: seqshop-server.sph.umich.edu
View Screenshot



• 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



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

View setup.txt

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

View Results

HG00641	/net/seqshop-server/home/hmkang/out/bams/HG00641.recal.bam
HG00640	/net/seqshop-server/home/hmkang/out/bams/HG00640.recal.bam
HG00551	/net/seqshop-server/home/hmkang/out/bams/HG00551.recal.bam
HG00553	/net/seqshop-server/home/hmkang/out/bams/HG00553.recal.bam
HG00554	bams/HG00554.recal.bam
HG00637	bams/HG00637.recal.bam
HG00638	bams/HG00638.recal.bam
HG00734	bams/HG00734.recal.bam
HG00736	bams/HG00736.recal.bam
HG00737	bams/HG00737.recal.bam

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

wc -l ${OUT}/bam.list

Your expected output is similar to this.

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

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

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

1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz
1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz.tbi
1000G.snps_indels.22.sites.bcf
1000G.snps_indels.22.sites.bcf.csi
1kg.omni.chr22.36Mb.vcf.gz
1kg.pilot_release.merged.indels.sites.hg19.chr22.vcf
dbsnp.13147541variants.22.sites.bcf
dbsnp.13147541variants.22.sites.bcf.csi
dbsnp_135.b37.chr22.vcf.gz
dbsnp_135.b37.chr22.vcf.gz.tbi
dbsnp.bcf.vcf
dbsnp.vcf.gz.vcf
gencode.cds.22.bed.gz
hapmap_3.3.b37.sites.chr22.vcf.gz
hapmap_3.3.b37.sites.chr22.vcf.gz.tbi
human.g1k.v37.chr22-bs.umfa
human.g1k.v37.chr22.dict
human.g1k.v37.chr22.fa
human.g1k.v37.chr22.fa.amb
human.g1k.v37.chr22.fa.ann
human.g1k.v37.chr22.fa.bwt
human.g1k.v37.chr22.fa.fai
human.g1k.v37.chr22.fa.pac
human.g1k.v37.chr22.fa.sa
human_g1k_v37.chr22.mask.100.fasta
human_g1k_v37.chr22.mask.100.fasta.fai
human.g1k.v37.chr22.winsize100.gc
humgen_g1k_v37_ploidy.chr22.map
indel.reference.txt
mdust.22.bed.gz
mills.208620indels.22.sites.bcf
mills.208620indels.22.sites.bcf.csi
mills_indels_hg19.22.sites.bcf

Parameters files required just for Structural Variation:

ls ${GC}/src/svtoolkit/conf

View Results

genstrip_parameters.txt

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 5 ${SS}/gotcloud.conf

View Results

##############################
## GenomeSTRIP
#############################
GENOMESTRIP_MASK_FASTA = $(REF_DIR)/human_g1k_v37.chr22.mask.100.fasta
GENOMESTRIP_PLOIDY_MAP = $(REF_DIR)/humgen_g1k_v37_ploidy.chr22.map

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

computerc.args.list
cpt
depth
depth.args.list
depth.dat
gcprofile
gcprofiles.list
gcprofiles.zip
genome_sizes.txt
isd
isd.dist.args.list
isd.dist.bin
rccache
rccache.bin  
rccache.bin.idx  
rccache.list  
rccache.merge  
spans  
spans.args.list 
 spans.dat

WAIT!!!!! DO NOT RUN THIS COMMAND, because it will take >1 hour to finish.

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

ls ${OUT}/sv/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

perl ${GC}/bin/genomestrip.pl --run-discovery --metadata ${SS}/metadata --conf ${SS}/gotcloud.conf --numjobs 2 --conf ${SS}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000 --base-prefix ${SS} --outdir ${OUT}

• ${GC}/bin/genomestrip.pl -run-discovery runs the GenomeSTRiP Discovery Pipeline
• --metadata ${SS}/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

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 ${GC}/bin/genomestrip.pl --run-genotype --metadata ${SS}/metadata --conf ${SS}/gotcloud.conf --numjobs 2 --base-prefix ${SS} --outdir ${OUT}

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

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

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?

Starting SNP Call on your own Genome

Go to SeqShop: Calling Your Own Genome, December 2014 so we can run SNP calling overnight.