3,045
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchSeqShop: Calling Your Own Genome, May 2015 (view source)
Revision as of 12:10, 22 May 2015
, 12:10, 22 May 2015→Annotation / Lookup against dbSNP
export HK=/net/seqshop-server/home/hmkang/apigenome/bin
export HK=/net/seqshop-server/home/hmkang/apigenome/bin
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts/
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 ====
==== Annotation / Lookup against dbSNP ====
Looking up SNPs by rsID is possible by (for example, rs17766217) -- How can we find its position?
Looking up SNPs by rsID is possible by (for example, rs17766217) -- How can we find its position?
$HK/tabix ~/NA12878/output/vcfs/chr8/chr8.filtered.rsid.vcf.gz 8:128504497 | less
$HK/tabix $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.
* 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.
And they can be combined as follows
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
(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 ====
==== Exonic Variants NOT found by 1000G ====
Want to see this from the BAM file? Use samtools tview:
Want to see this from the BAM file? Use samtools tview:
$GC/bin/samtools tview $SAMPLE/output/bams/$SAMPLE.recal.bam $GC/gotcloud.ref/human.g1k.v37.fa
$GC/bin/samtools tview $SAMPLE/output/bams/$SAMPLE.recal.bam $REF/hs37d5.fa
Use 'g' & enter the Chr:Pos
Use 'g' & enter the Chr:Pos
* Some patterns may indicate not real variants.
* Some patterns may indicate not real variants.
The phred score at the last column quantifies the degree of functional significance
The phred score at the last column quantifies the degree of functional significance
