Difference between revisions of "SeqShop: Analysis of Structural Variation Practical, June 2014"
Line 38: | Line 38: | ||
* 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 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 | ||
+ | |||
+ | 62 /net/seqshop-server/hmkang/out/bam.index | ||
− | |||
=== Reference Files === | === Reference Files === | ||
Line 71: | Line 98: | ||
<div class="mw-collapsible mw-collapsed" style="width:200px"> | <div class="mw-collapsible mw-collapsed" style="width:200px"> | ||
<li>View Results</li> | <li>View Results</li> | ||
− | <div class="mw-collapsible-content" style="width: | + | <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 | 1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz hapmap_3.3.b37.sites.chr22.vcf.gz.tbi human.g1k.v37.chr22.fa.bwt | ||
1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz.tbi human.g1k.v37.chr22-bs.umfa human.g1k.v37.chr22.fa.fai | 1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz.tbi human.g1k.v37.chr22-bs.umfa human.g1k.v37.chr22.fa.fai | ||
Line 87: | Line 114: | ||
ls $SV/ref | ls $SV/ref | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
<ul> | <ul> | ||
<div class="mw-collapsible mw-collapsed" style="width:200px"> | <div class="mw-collapsible mw-collapsed" style="width:200px"> | ||
− | <li> | + | <li>View Results</li> |
− | <div class="mw-collapsible-content"> | + | <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> | ||
</div> | </div> | ||
− | |||
</ul> | </ul> | ||
− | |||
− | |||
− | |||
− | + | ls $SV/conf | |
− | |||
− | |||
− | |||
− | |||
− | |||
<ul> | <ul> | ||
− | <div class="mw-collapsible mw-collapsed" style="width: | + | <div class="mw-collapsible mw-collapsed" style="width:200px"> |
− | <li> | + | <li>View Results</li> |
− | <div class="mw-collapsible-content"> | + | <div class="mw-collapsible-content" style="width:800px"> |
− | + | genstrip_parameters.txt humgen_g1k_v37_ploidy.chr22.map humgen_g1k_v37_ploidy.map | |
− | |||
− | |||
− | |||
− | |||
</div> | </div> | ||
</div> | </div> | ||
</ul> | </ul> | ||
− | |||
+ | === GotCloud Configuration File === | ||
+ | We will use a slightly modified version of configuration file as we used yesterday in GotCloud Align. | ||
+ | |||
+ | See [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#GotCloud Configuration File|SeqShop: Alignment: GotCloud Configuration File]] for more details | ||
+ | * 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]] | |
− | * | + | * Contains information on how to configure for exome/targeted sequencing |
− | |||
− | |||
− | + | Check out what was changed. | |
− | |||
− | + | cat $IN/gotcloud.conf | |
− | |||
− | |||
<ul> | <ul> | ||
− | <div class="mw-collapsible mw-collapsed" style="width: | + | <div class="mw-collapsible mw-collapsed" style="width:200px"> |
− | <li> | + | <li>View Results</li> |
− | <div class="mw-collapsible-content"> | + | <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_SVTOOLKIT_DIR = $(GOTCLOUD_ROOT)/../reference/svtoolkit | ||
+ | GENOMESTRIP_MASK_FASTA = $(GENOMESTRIP_SVTOOLKIT_DIR)/ref/human_g1k_v37.chr22.mask.100.fasta | ||
+ | GENOMESTRIP_PLOIDY_MAP = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/humgen_g1k_v37_ploidy.chr22.map | ||
+ | GENOMESTRIP_PARAM = $(GENOMESTRIP_SVTOOLKIT_DIR)/conf/genstrip_parameters.txt | ||
</div> | </div> | ||
</div> | </div> | ||
</ul> | </ul> | ||
− | + | ||
− | + | == Run GotCloud SnpCall == | |
− | + | [[File:SnpcallDiagram.png|500px]] | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | == Run GotCloud SnpCall == | ||
− | [[File:SnpcallDiagram.png|500px]] | ||
− | Now that we have all of our input files, we need just a simple command to run: | + | Now that we have all of our input files, we need just a simple command to run: |
− | + | $GC/gotcloud snpcall --conf $IN/gotcloud.conf --numjobs 4 --region 22:36000000-37000000 | |
* --numjobs tells GotCloud how many jobs to run in parallel | * --numjobs tells GotCloud how many jobs to run in parallel | ||
** Depends on your system | ** Depends on your system |
Revision as of 14:56, 17 June 2014
(WARNING: Under Construction)
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.
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. It 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 set up, type
less /home/mktrost/seqshop/setup.txt
and press 'q' to finish.
Examining GotCloud/GenomeSTRiP Input files
Sequnce Alignment Files: BAM Files and Index Files
The GotCloud Indel 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
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 yesterday, but you can take another look at the chromosome 22 reference files included for this tutorial:
ls $REF
ls $REF
- View Results
1000G_omni2.5.b37.sites.PASS.chr22.vcf.gz hapmap_3.3.b37.sites.chr22.vcf.gz.tbi human.g1k.v37.chr22.fa.bwt 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
Additional reference and parameters
ls $SV/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
ls $SV/conf
- View Results
genstrip_parameters.txt humgen_g1k_v37_ploidy.chr22.map humgen_g1k_v37_ploidy.map
GotCloud Configuration File
We will use a slightly modified version of configuration file as we used yesterday in GotCloud Align.
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
- Contains information on how to configure for exome/targeted sequencing
Check out what was changed.
cat $IN/gotcloud.conf
- View Results
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_SVTOOLKIT_DIR = $(GOTCLOUD_ROOT)/../reference/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
== Run GotCloud SnpCall ==
Now that we have all of our input files, we need just a simple command to run: $GC/gotcloud snpcall --conf $IN/gotcloud.conf --numjobs 4 --region 22:36000000-37000000
- --numjobs tells GotCloud 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 GotCloud to ignore the other regions
This should take about 5 minutes to run.
- After about 4 minutes of running, GotCloud snpcall will output some text to the screen. Don't worry, that is expected and is just output from some of the intermediate tools.
- It should end with a line like:
Commands finished in 289 secs with no errors reported
If you cancelled GotCloud part way through, just rerun your GotCloud command and it will pick up where it left off.
Examining GotCloud SnpCall Output
Let's look at the output directory:
ls $OUT
- Do you see any new files or directories?
Let's look at the vcfs directory:
ls $OUT/vcfs
Just a chr22
directory, so look inside of there:
ls $OUT/vcfs/chr22
- Can you identify the final filtered VCF and the associated summary file?
- Answer & annotated directory listing:
Now, let's look in the split directory for the VCF with just the passing variants:
ls $OUT/split/chr22
- Which file do you think is the one you want?
Filtering Summary Statistics
cat $OUT/vcfs/chr22/chr22.filtered.sites.vcf.summary
To understand how to interpret the filtering summary statistics, please refer to Understanding vcf-summary output
Filtered VCF
Let's look at the filtered sites file.
less -S $OUT/vcfs/chr22/chr22.filtered.sites.vcf
- Scroll down until you find some variants.
- Use space bar to jump a full page
- Use down arrow to move down one line
- Scroll right: lots of info fields, but no per sample genotype information
- What is the first filtered out variant that you find & what filter did it fail?
Remember, use 'q'
to exit out of less
q
Now, let's look at the filtered file with genotypes.
zless -S $OUT/vcfs/chr22/chr22.filtered.vcf.gz
- Scroll down until you find some variants.
- Use space bar to jump a full page
- Use down arrow to move down one line
- Scroll right until you should see per sample genotype information
Passing SNPs
Let's look at the file of just the pass sites:
zless -S $OUT/split/chr22/chr22.filtered.PASS.vcf.gz
- Scroll down: they all look like they
PASS
Let's check if they are all PASS.
zcat $OUT/split/chr22/chr22.filtered.PASS.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
It will return nothing since there are no non-passing variants in this file.
- Want an explanation of this command?
- zcat ...: uncompress the zipped VCF
- '|' : this takes the output of one command and sends it as input to the next
- grep -v "^#" : exclude any lines that start with "#" - headers
- cut -f 7 : extract the FILTER column (the 7th column)
- grep -v "PASS" : exclude any rows that have a "PASS" in the FILTER column
Compare that to the filtered file we looked at before:
zcat $OUT/vcfs/chr22/chr22.filtered.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
- Do you see any filters?
- Answer
- Yes
- It should have scrolled and you should see filters like:
- INDEL5;SVM
- INDEL5
- SVM
- It should have scrolled and you should see filters like:
GotCloud Genotype Refinement
To improve the quality of the genotypes, we run a genotype refinement pipeline.
This pipeline runs Beagle & thunder.
Genotype Refinement Input
The GotCloud genotype refinement pipeline takes as input $OUT/split/chr22/chr22.filtered.PASS.vcf.gz (the VCF file of PASS'ing SNPs from snpcall.
The bam index and the configuration file we used for GotCloud snpcall will tell GotCloud genotype refinement everything it needs to know, so no new input files need to be prepared.
Note: the configuration file overrides the THUNDER command to make it go faster than the default settings so the tutorial will run faster:
Running GotCloud Genotype Refinement
Since everything is setup, just run the following command (very similar to snpcall).
$GC/gotcloud ldrefine --conf $IN/gotcloud.conf --numjobs 2 --region 22:36000000-37000000
- Beagle will take about 2-3 minutes to complete
- Thunder will automatically run and will take another 3-4 minutes
Genotype Refinement Output
- What's new in the output directory?
- Answer
beagle
directory : Beagle outputthunder
directory : Thunder outputumake.beagle.conf
: Configuration values used for GotCloud beagleumake.beagle.Makefile
: GNU makefile for commands run as part of GotCloud beagleumake.beagle.Makefile.log
: Log of the GotCloud beagle runumake.thunder.*
files : Contain the configuration & steps used in GotCloud thunder
Let's take a look at that interesting location we found in the alignment tutorial : chromosome 22, positions 36907000-36907100
Use tabix to extract that from the VCFs:
$GC/bin/tabix $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
Remember, type 'q' to quit less.
q
- Are there any variants in this region?
- Answer:
- Yes!
- Positions:
36907001
; Ref: T, Alt: C - that's what we saw before36907098
; Ref: T, Alt: C - that's what we saw before
- What is HG00551's genotype at these positions?
- First check which sample number HG00551 is:
$GC/bin/tabix -H $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz
- That will help you figure out it's genotype.
- Rerun the tabix command and scroll to find HG00551's genotype:
$GC/bin/tabix $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
- Answer:
- It is the first sample
0|1
: Heterozygous1|1
; Homozygous Alt (C)
Remember, type 'q' to quit less.
q
Did I find interesting variants?
The region we selected contains APOL1 gene, which is known to play an important role in kidney diseases such as nephrotic syndrome. One of the non-synonymous risk allele, rs73885139
located at position 22:36661906
increases the risk of nephrotic syndrome by >2-folds. Let's see if we found the interesting variant by looking at the VCF file by position.
$GC/bin/tabix $OUT/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1
Did you see a variant at the position?