Changes

From Genome Analysis Wiki
Jump to: navigation, search

SeqShop: Calling Your Own Genome, May 2015

3,714 bytes removed, 12:10, 22 May 2015
Annotation / Lookup against dbSNP
export SAMPLE=NA12878
Point to your GotCloud & your output directory:
export GC=~/seqshop/gotcloud
export OUT=~/$SAMPLE/output
* Relative path, so assumes running from your home directory (I prefer absolute paths, but for simplicity of the workshop, we just use relative path).
=== Configuring Indel SNPCALL ===No special Configuration settings for SNP calling
cat ~/$SAMPLE/gotcloud.conf
** 2-3 of you on the machine - 3*8 = 24 jobs will be running in parallel on that machine
${GC}/gotcloud snpcall --conf $SAMPLE/gotcloud.conf --numjobs 8 --outdir $OUT
* Only need the configuration & , number of threads, and the output directory, rest is specified within the configuration.
This will run overnight. We will check if it completed at the practical in the morning.
exit PuTTY
 
 
=== FEEDBACK! ===
Please provide feedback on today:
 
https://docs.google.com/forms/d/1ADTkBjzT-QNj2lrejyqGqDaahTponrw20kSgDNwqwH4/viewform
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
== Wednesday Thursday ==
<div class="mw-collapsible-content">
<div class="mw-collapsible" style="width:500px">=== Checking if INDEL snpcall Completed ===<div class="mw-collapsible-content">==== Logging Back in Resume screen to Check Jobs ====
;How do you log back into screen?
This will resume an already running screen.
Verify Your screen session still has your environment variables set, so you got a "completed successfully" messagedo not need to reset them.
How long did INDEL calling take? Look at the log message - time in seconds.
=== Detach From screen===Detach from screenVerify you got a "completed successfully" message. We will resume it again later when we start SNPCall. Ctrl-a d
</div></div>How long did snpcall calling take? Look at the log message - time in seconds.
=== Structural Variation Tutorial ===
Now we are going to run the Structural Variation Practical
 
Please go to: [[SeqShop: Analysis of Structural Variation Practical, May 2015]]
 
We will Start SNP Calling after the practical.
 
 
<div class="mw-collapsible" style="width:500px">
 
=== Start SNP Calling ===
<div class="mw-collapsible-content">
==== Resume screen ====
 
;How do you log back into screen?
screen -r
This will resume an already running screen.
 
Your screen session still has your environment variables set, so you do not need to reset them.
==== List of BAMs ====
* Relative path, so assumes running from your home directory (I prefer absolute paths, but for simplicity of the workshop, we just use relative path).
 ==== GotCloud SnpCall INDEL Configuration ====
cat ~/$SAMPLE/gotcloud.conf
You will see something like this:<pre># Cluster SettingsBATCH_TYPE = BATCH_OPTS =  OUT_DIR = Sample13/output # Align SettingsMAP_TYPE = BWA_MEMBWA_THREADS = -t 24FASTQ_LIST = fastq.list # SNP Call SettingsUNIT_CHUNK = 20000000 # Chunk size of SNP Same as it looked the other day with no special Configuration settings for INDEL calling : 20MbVCF_EXTRACT = /net/seqshop-server/home/mktrost/seqshop/singleSample/snpOnly.vcf.gzMODEL_GLFSINGLE = TRUEMODEL_SKIP_DISCOVER = FALSEMODEL_AF_PRIOR = TRUE EXT_DIR = /net/seqshop-server/home/mktrost/seqshop/singleSample/extEXT = $(EXT_DIR)/ALL.chrCHR.phase3.combined.sites.unfiltered.vcf.gz $(EXT_DIR)/chrCHR.filtered.sites.vcf.gz</pre>
===== Configuration Updates ====='''No configuration updates from the original settings are necessary.'''* Originally you were going to update the configuration to do exome calling only, but we have decided to do whole genome** If it doesn't finish tonight, we will kill it at the tutorial in the morning, and take advantage of the GotCloud restart capability after the tutorial tomorrow. If you started making updates to your configuration yesterday, you can try to make it match above (mktrost is actually the path you want in the paths above).* The only difference from above should be your OUT_DIR.Check it with (leave as ~mktrost - you are comparing to my NA12878 conf file): diff $SAMPLE/gotcloud.conf ~mktrost/NA12878/gotcloud.conf '''If you see differences other than OUT_DIR, you can modify your file using nedit (or your favorite editor as you used yesterday).'''* '''Or you can copy my file and just change OUT_DIR:''' cp ~mktrost/NA12878/gotcloud.conf $SAMPLE/gotcloud.conf nedit $SAMPLE/gotcloud.conf&* Edit OUT_DIR OUT_DIR = Sample##/output* Save, and close ==== Running SnpCall INDEL ====Run GotCloud snpcall indel with 6 jobs running in parallel* Why 6? ** You want to run as many as you can.** 5 of you on the machine - 5*6 = 30 jobs will be running in parallel on that machine ${GC}/gotcloud snpcall indel --conf $SAMPLE/gotcloud.conf --numjobs 6--outdir $OUT* Only need the configuration & , number of threads, and the output directory, rest is specified within the configuration.
This will run overnight. We will check if it completed at the practical in the morning.
exit PuTTY
=== FEEDBACK!===
Please provide feedback for today.
https://docs.google.com/a/umich.edu/forms/d/1iES6usHxLB7Ec9hRxtqYgH7v05lU3Ume4VJcksx8Ogg/viewform
</div>
</div>
</div>
</div>
 
<div class="mw-collapsible mw-collapsed" style="width:500px">
 
== Thursday ==
<div class="mw-collapsible-content">
=== Checking if SnpCall Completed ===
==== Logging Back in to Check Jobs ====
 
;How do you log back into screen?
screen -r
This will resume an already running screen.
 
==== Checking Completion ====
 
Did you get a "completed successfully" message?
 
If yes, how long did SNP calling take? Look at the log message - time in seconds.
 
If no, are you running on seqshop-server? If so, KILL it. ssh -X to one of the seqshop machines and run on there.
Ctrl-c
 
Detach from screen. We will resume it again later when we restart SNPCall (if necessary).
Ctrl-a d
 
=== Ancestry Tutorial ===
Now we are going to run the Structural Variation Practical
 
Please go to: [[SeqShop: Estimates of Genetic Ancestry Practical, May 2015]]
 
We will Resume SNP Calling (if necessary) after the practical.
 
=== Restart SnpCall ===
==== Resume screen ====
 
;How do you log back into screen?
screen -r
This will resume an already running screen.
 
Your screen session still has your environment variables set, so you do not need to reset them.
 
==== Running SnpCall ====
Run GotCloud snpcall with 6 jobs running in parallel
${GC}/gotcloud snpcall --conf $SAMPLE/gotcloud.conf --numjobs 6
* GotCloud will pick up after the last completed step from before.
 
This will run overnight. We will check if it completed at the practical in the morning.
 
==== Log Out ====
;Want to log out and leave your job running?
In the screen window, type:
Ctrl-a d
(Hold down Ctrl and type 'a', let go of both and type 'd')
* This will "detach" from your screen session while your alignment continues to run.
 
exit PuTTY
 
 
</div>
</div>
[[SeqShop: Ancestry On Your Own Genome, May 2015]]
=== Association Analysis Tutorial Setup Variables ===Now we are going Set these values. Also, be sure to run the Association Analysis Practicalspecify your sample name instead of SampleXX export SAMPLE=SampleXXor export SAMPLE=NA12878
Please go Point toyour GotCloud & your output directory: [[SeqShop: Association Analysis, May 2015]] export GC=~/seqshop/gotcloudWe will look at our own genomes again after the practical. === Return to SeqShop: Ancestry On Your Own Genome, May 2015 = export OUT==Return to [[SeqShop:_Ancestry_On_Your_Own_Genome,_May_2015#Checking_if_Pileup_finished]]~/$SAMPLE/output
=== Reviewing Indel Results ===
Set these values. Also, be sure to specify your sample name (or NA12878) instead of SampleXX
export SAMPLE=SampleXX
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
 
Look in the output directory
ls ~/$SAMPLE/output
The insertion deletion ratio increases from 0.91 to 0.92.
=== Return to SeqShop: Ancestry On Your Own Genome, May 2015 ===
Return to [[SeqShop:_Ancestry_On_Your_Own_Genome,_May_2015#Checking_if_Pileup_finished]]
=== Friday: Reviewing SNPCALL Results ===
=== Friday : More SNP Analysis ===
In addition, set another environmental variable for locating the binaries for custom analysis
export HK=/net/seqshop-server/home/hmkang/apigenome/bin export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts/ export REF=/net/seqshop-server/home/mktrost/seqshop/singleSample/ref/gotcloud.ref export GC= Environmental Variables ====~/seqshop/gotcloud
If you didn't set the environmental variable, you can set it again
source /net/seqshop-server/home/mktrost/seqshop/setup.txt export SAMPLE=SampleXX (MAKE SURE TO CHANGE XX to your number or use NA12878 instead) source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt In addition, set another environmental variable for locating the binaries for custom analysis  export HKOUT=~/net/seqshop-server/home/hmkang/seqshop$SAMPLE/binoutput
==== Annotation / Lookup against dbSNP ====
If you want to add rsIDs to your variant files, you can do this by running the following command
$HK/vcf-add-rsid -vcf $OUT/vcfs/chr1/chr1.filtered.vcf.gz --db $HK/../data/dbSNP.b138/dbsnp_138dbsnp_142.b37.vcf.gz --out $OUT/vcfs/chr1/chr1.filtered.rsid.vcf.gz
If you want to run this command across all chromosomes in parallel, you can use the special script run-command-wgs
$HK/run-commandmake -wgs -repeat-chr --cmd "$HK/vcf-add-rsid -vcf $OUT/vcfs/chr1/chr1.filtered.vcf.gz --db $HK/../data/dbSNP.b138/dbsnp_138dbsnp_142.b37.vcf.gz --out $OUT/vcfs/chr1/chr1.filtered.rsid.vcf.gz" --numjobs 6--out runmake.rsid
Looking up SNPs by rsID is possible by (for example, rs17766217)-- How can we find its position? $HK/vcf-lookup-rsid --vcf tabix $OUT/vcfs/chr1chr8/chr1chr8.filtered.rsid.vcf.gz --sepchr --rs rs177662178:128504497 | less
* Be sure to look at the QUAL & your sample's PL, and not just the GL field. Check if QUAL is 0 or PL is 0,0,0 - NS is also probably 0; DP is probably 0. That means you probably didn't have any copies, so your GT may not be correct/is unknown.
If you want to browse the rsIDs of known GWAS SNPs, you can do this by
cut -f 1,8,12,13,22 $HK/../data/gwascatalog/gwascatalog.txt | lessgrep -w rs17766217
==== Annotating your genome ====
Or you can run multiple chromosomes in parallel in one command
$HK/run-commandmake -wgs -repeat-chr --cmd "$EPACTS/bin/epacts anno --in $OUT/vcfs/chr1/chr1.filtered.rsid.vcf.gz --out $OUT/vcfs/chr1/chr1.filtered.rsid.anno.vcf.gz" --numjobs 6 --out runmake.anno
==== Extracting only exonic SNPs ====
If you want to look at the exonic SNPs, you can extract using the following command
$HK/run-commandmake -wgs -repeat-chr --cmd "($HK/tabix -H $OUT/vcfs/chr1/chr1.filtered.rsid.anno.vcf.gz; zcat $OUT/vcfs/chr1/chr1.filtered.rsid.anno.vcf.gz | grep Exon;)| $HK/bgzip -c > $OUT/vcfs/chr1/chr1.filtered.rsid.anno.exon.vcf.gz" --numjobs 6--out runmake.exome
And they can be combined as follows
(zcat $OUT/vcfs/chr1/chr1.filtered.rsid.anno.exon.vcf.gz; zcat $OUT/vcfs/chr[2-9]/chr*.filtered.rsid.anno.exon.vcf.gz $OUT/vcfs/chr??/chr*.filtered.rsid.anno.exon.vcf.gz $OUT/vcfs/chrX/chrX.filtered.rsid.anno.exon.vcf.gz | grep -v ^#) | $HK/bgzip -c > $OUT/wgs.filtered.rsid.anno.exon.vcf.gz
$HK/tabix -pvcf $OUT/wgs.filtered.rsid.anno.exon.vcf.gz
==== Exonic Variants NOT found by 1000G ====
Want to see this from the BAM file? Use samtools tview:
$GC/bin/samtools tview $SAMPLE/output/bams/$SAMPLE.recal.bam $GCREF/gotcloud.ref/human.g1k.v37hs37d5.fa
Use 'g' & enter the Chr:Pos
* Some patterns may indicate not real variants.
The phred score at the last column quantifies the degree of functional significance
 
</div>
</div>
 
== OVERALL COURSE FEEDBACK! ==
Please provide feedback:
https://docs.google.com/forms/d/1pxfPXKwWfA71ZJM99Sevs3MwAUz2UbHAR8dnRI-kRNM/viewform

Navigation menu