1,120
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchSeqShop: Analysis of Structural Variation Practical, June 2014 (view source)
Revision as of 01:41, 17 June 2014
, 01:41, 17 June 2014Created page with "==Introduction== See the introductory slides for an intro to this tutorial. == Goals of This Session == * What we want to learn is ca..."
==Introduction==
See the [[Media:SeqShop - GotCloud svcall.pdf|introductory slides]] for an intro to this tutorial.
== 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.
{{SeqShopLogin}}
== Setup your run environment==
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in. 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 <code>source</code>
** 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.
<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content">
[[File:setup.png|500px]]
</div>
</div>
== Examining GotCloud/GenomeSTRiP Input files ==
=== Sequnce Alignment Files: BAM Files ===
Per sample BAM files contain sequence reads that are mapped to positions in the genome.
For a reminder on how to look at/read BAM files, see: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#BAM_Files|SeqShop Aligment: BAM Files]]
For this tutorial, we will use the 4 BAMs produced in the [[SeqShop: Sequence Mapping and Assembly Practical]] as well as with 58 BAMs that were pre-aligned to that 1MB region of chromosome 22.
=== 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
For GotCloud/GenomeSTRiP, you need:
# Reference genome FASTA file
# 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}
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:RefDir.png|500px]]
</div>
</div>
</ul>
=== GotCloud BAM Index File ===
The BAM index file points GotCloud to the BAM files
* generated by the alignment pipeline
Look at the BAM index file the alignment pipeline generated
cat ${OUT}/bam.index
;What is the path to the BAM file for sample HG00640?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>/home/YourUserName/out/bams/HG00640.recal.bam</li>
[[File:Bamindex.png|500px]]
</div>
</div>
</ul>
</ul>
The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.
* The other 58 samples were already aligned:
ls ${IN}/bams
Look at the BAM index for those BAMs:
less ${IN}/bams/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
q
;Do you notice a difference between this index and yours?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>It doesn't have a full path to the BAM file, while your index has /home/...</li>
[[File:Bamindex1.png|300px]]
<li>That's ok, <code>gotcloud.conf</code> contains the path to those BAMs</li>
[[File:BamindexConf.png|300px]]
</div>
</div>
</ul>
</ul>
We need to add these BAMs to our index
* Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
cat ${IN}/bams/bam.index >> ${OUT}/bam.index
* ">>" will append to the file that follows it
Verify your BAM index contains the additional BAMs
less ${OUT}/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
q
;Do you see both sets of BAMs?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Annotated Screenshot:</li>
<div class="mw-collapsible-content">
<ul>
<li>If not, let me know</li>
[[File:Bamindex2.png|400px]]
</div>
</div>
</ul>
</ul>
=== GotCloud Configuration File ===
We will use the same 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
== Run GotCloud SnpCall ==
[[File:SnpcallDiagram.png|500px]]
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
<div class="mw-collapsible mw-collapsed" style="width:500px">
Curious if it started running properly? Check out this screenshot:
<div class="mw-collapsible-content">
[[File:SnpcallStart.png|750px]]
</div>
</div>
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: <code>Commands finished in 289 secs with no errors reported</code>
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?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* View Annotated Screenshot:
<div class="mw-collapsible-content">
[[File:gcsnpcallOut.png|500px]]
</div>
</div>
Let's look at the vcfs directory:
ls ${OUT}/vcfs
Just a <code>chr22</code> directory, so look inside of there:
ls ${OUT}/vcfs/chr22
;Can you identify the final filtered VCF and the associated summary file?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* Answer & annotated directory listing:
<div class="mw-collapsible-content">
<ul>
<li>Filtered VCF (SVM & hard filters): chr22.filtered.vcf.gz</li>
<li>Summary file: chr22.filtered.sites.vcf.summary</li>
</ul>
[[File:vcfsout.png|600px]]
</div>
</div>
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?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* Answer:
<div class="mw-collapsible-content">
<ul>
<li>chr22.filtered.PASS.vcf.gz</li>
</ul>
[[File:splitOut.png|600px]]
</div>
</div>
=== Filtering Summary Statistics ===
cat ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf.summary
<div class="mw-collapsible mw-collapsed" style="width:250px">
View Screenshot
<div class="mw-collapsible-content">
[[File:filterSum.png]]
</div>
</div>
'''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?
<div class="mw-collapsible mw-collapsed" style="width:250px">
* Answer:
<div class="mw-collapsible-content">
<ul>
<li>It failed SVM filter</li>
</ul>
[[File:SvmFilt.png|550px]]
</div>
</div>
Remember, use <code>'q'</code> to exit out of <code>less</code>
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
<div class="mw-collapsible mw-collapsed" style="width:250px">
* View annotated screenshot:
<div class="mw-collapsible-content">
[[File:SvmFiltGL.png|550px]]
</div>
</div>
=== 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 <code>PASS</code>
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.
<div class="mw-collapsible mw-collapsed" style="width:450px">
;Want an explanation of this command?
<div class="mw-collapsible-content">
* 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
</div>
</div>
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?
<div class="mw-collapsible mw-collapsed" style="width:450px">
*Answer
<div class="mw-collapsible-content">
* Yes
** It should have scrolled and you should see filters like:
*** INDEL5;SVM
*** INDEL5
*** SVM
</div>
</div>
== GotCloud Genotype Refinement ==
To improve the quality of the genotypes, we run a genotype refinement pipeline.
This pipeline runs [http://faculty.washington.edu/browning/beagle/beagle.html 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:
[[File:thunderConf.png|600px]]
=== 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?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer</li>
<div class="mw-collapsible-content">
<ul>
<li><code>beagle</code> directory : Beagle output</li>
<li><code>thunder</code> directory : Thunder output</li>
<li><code>umake.beagle.conf</code> : Configuration values used for GotCloud beagle</li>
<li><code>umake.beagle.Makefile</code> : GNU makefile for commands run as part of GotCloud beagle</li>
<li><code>umake.beagle.Makefile.log</code> : Log of the GotCloud beagle run</li>
<li><code>umake.thunder.*</code> files : Contain the configuration & steps used in GotCloud thunder</li>
</ul>
</div>
</div>
</ul>
Let's take a look at that interesting location we found in the [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Accessing_BAMs_by_Position|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?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>Yes!</li>
<li>Positions:</li>
<ul>
<li><code>36907001</code>; Ref: T, Alt: C - that's what we saw before</li>
<li><code>36907098</code>; Ref: T, Alt: C - that's what we saw before</li>
</ul>
</ul>
</div>
</div>
</ul>
;What is HG00551's genotype at these positions?
#First check which sample number HG00551 is:
zcat ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz |grep "#CHROM"
* 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
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>It is the first sample</li>
<li><code>0|1</code>: Heterozygous</li>
<li><code>1|1</code>; Homozygous Alt (C)</li>
</ul>
</div>
</div>
</ul>
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, <code>rs73885139</code> located at position <code>22:36661906</code> 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?
See the [[Media:SeqShop - GotCloud svcall.pdf|introductory slides]] for an intro to this tutorial.
== 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.
{{SeqShopLogin}}
== Setup your run environment==
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in. 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 <code>source</code>
** 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.
<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content">
[[File:setup.png|500px]]
</div>
</div>
== Examining GotCloud/GenomeSTRiP Input files ==
=== Sequnce Alignment Files: BAM Files ===
Per sample BAM files contain sequence reads that are mapped to positions in the genome.
For a reminder on how to look at/read BAM files, see: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#BAM_Files|SeqShop Aligment: BAM Files]]
For this tutorial, we will use the 4 BAMs produced in the [[SeqShop: Sequence Mapping and Assembly Practical]] as well as with 58 BAMs that were pre-aligned to that 1MB region of chromosome 22.
=== 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
For GotCloud/GenomeSTRiP, you need:
# Reference genome FASTA file
# 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}
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:RefDir.png|500px]]
</div>
</div>
</ul>
=== GotCloud BAM Index File ===
The BAM index file points GotCloud to the BAM files
* generated by the alignment pipeline
Look at the BAM index file the alignment pipeline generated
cat ${OUT}/bam.index
;What is the path to the BAM file for sample HG00640?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>/home/YourUserName/out/bams/HG00640.recal.bam</li>
[[File:Bamindex.png|500px]]
</div>
</div>
</ul>
</ul>
The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.
* The other 58 samples were already aligned:
ls ${IN}/bams
Look at the BAM index for those BAMs:
less ${IN}/bams/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
q
;Do you notice a difference between this index and yours?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>It doesn't have a full path to the BAM file, while your index has /home/...</li>
[[File:Bamindex1.png|300px]]
<li>That's ok, <code>gotcloud.conf</code> contains the path to those BAMs</li>
[[File:BamindexConf.png|300px]]
</div>
</div>
</ul>
</ul>
We need to add these BAMs to our index
* Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
cat ${IN}/bams/bam.index >> ${OUT}/bam.index
* ">>" will append to the file that follows it
Verify your BAM index contains the additional BAMs
less ${OUT}/bam.index
Remember, use <code>'q'</code> to exit out of <code>less</code>
q
;Do you see both sets of BAMs?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Annotated Screenshot:</li>
<div class="mw-collapsible-content">
<ul>
<li>If not, let me know</li>
[[File:Bamindex2.png|400px]]
</div>
</div>
</ul>
</ul>
=== GotCloud Configuration File ===
We will use the same 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
== Run GotCloud SnpCall ==
[[File:SnpcallDiagram.png|500px]]
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
<div class="mw-collapsible mw-collapsed" style="width:500px">
Curious if it started running properly? Check out this screenshot:
<div class="mw-collapsible-content">
[[File:SnpcallStart.png|750px]]
</div>
</div>
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: <code>Commands finished in 289 secs with no errors reported</code>
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?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* View Annotated Screenshot:
<div class="mw-collapsible-content">
[[File:gcsnpcallOut.png|500px]]
</div>
</div>
Let's look at the vcfs directory:
ls ${OUT}/vcfs
Just a <code>chr22</code> directory, so look inside of there:
ls ${OUT}/vcfs/chr22
;Can you identify the final filtered VCF and the associated summary file?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* Answer & annotated directory listing:
<div class="mw-collapsible-content">
<ul>
<li>Filtered VCF (SVM & hard filters): chr22.filtered.vcf.gz</li>
<li>Summary file: chr22.filtered.sites.vcf.summary</li>
</ul>
[[File:vcfsout.png|600px]]
</div>
</div>
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?
<div class="mw-collapsible mw-collapsed" style="width:350px">
* Answer:
<div class="mw-collapsible-content">
<ul>
<li>chr22.filtered.PASS.vcf.gz</li>
</ul>
[[File:splitOut.png|600px]]
</div>
</div>
=== Filtering Summary Statistics ===
cat ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf.summary
<div class="mw-collapsible mw-collapsed" style="width:250px">
View Screenshot
<div class="mw-collapsible-content">
[[File:filterSum.png]]
</div>
</div>
'''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?
<div class="mw-collapsible mw-collapsed" style="width:250px">
* Answer:
<div class="mw-collapsible-content">
<ul>
<li>It failed SVM filter</li>
</ul>
[[File:SvmFilt.png|550px]]
</div>
</div>
Remember, use <code>'q'</code> to exit out of <code>less</code>
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
<div class="mw-collapsible mw-collapsed" style="width:250px">
* View annotated screenshot:
<div class="mw-collapsible-content">
[[File:SvmFiltGL.png|550px]]
</div>
</div>
=== 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 <code>PASS</code>
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.
<div class="mw-collapsible mw-collapsed" style="width:450px">
;Want an explanation of this command?
<div class="mw-collapsible-content">
* 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
</div>
</div>
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?
<div class="mw-collapsible mw-collapsed" style="width:450px">
*Answer
<div class="mw-collapsible-content">
* Yes
** It should have scrolled and you should see filters like:
*** INDEL5;SVM
*** INDEL5
*** SVM
</div>
</div>
== GotCloud Genotype Refinement ==
To improve the quality of the genotypes, we run a genotype refinement pipeline.
This pipeline runs [http://faculty.washington.edu/browning/beagle/beagle.html 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:
[[File:thunderConf.png|600px]]
=== 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?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer</li>
<div class="mw-collapsible-content">
<ul>
<li><code>beagle</code> directory : Beagle output</li>
<li><code>thunder</code> directory : Thunder output</li>
<li><code>umake.beagle.conf</code> : Configuration values used for GotCloud beagle</li>
<li><code>umake.beagle.Makefile</code> : GNU makefile for commands run as part of GotCloud beagle</li>
<li><code>umake.beagle.Makefile.log</code> : Log of the GotCloud beagle run</li>
<li><code>umake.thunder.*</code> files : Contain the configuration & steps used in GotCloud thunder</li>
</ul>
</div>
</div>
</ul>
Let's take a look at that interesting location we found in the [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Accessing_BAMs_by_Position|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?
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>Yes!</li>
<li>Positions:</li>
<ul>
<li><code>36907001</code>; Ref: T, Alt: C - that's what we saw before</li>
<li><code>36907098</code>; Ref: T, Alt: C - that's what we saw before</li>
</ul>
</ul>
</div>
</div>
</ul>
;What is HG00551's genotype at these positions?
#First check which sample number HG00551 is:
zcat ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz |grep "#CHROM"
* 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
<ul>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Answer:</li>
<div class="mw-collapsible-content">
<ul>
<li>It is the first sample</li>
<li><code>0|1</code>: Heterozygous</li>
<li><code>1|1</code>; Homozygous Alt (C)</li>
</ul>
</div>
</div>
</ul>
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, <code>rs73885139</code> located at position <code>22:36661906</code> 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?
