From Genome Analysis Wiki
Jump to: navigation, search

SeqShop: Calling Your Own Genome, May 2015

239 bytes added, 12:10, 22 May 2015
Annotation / Lookup against dbSNP
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=~/seqshop/gotcloud
export SAMPLE=SampleXX
export OUT=~/$SAMPLE/output
==== Annotation / Lookup against dbSNP ====
Looking up SNPs by rsID is possible by (for example, rs17766217) -- How can we find its position?
$HK/tabix ~/NA12878/output$OUT/vcfs/chr8/chr8.filtered.rsid.vcf.gz 8: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.
Or you can run multiple chromosomes in parallel in one command
$HK/run-make --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

Navigation menu