Ylaine/2009-8-3

From ZhangLabWiki
Jump to navigation Jump to search
  • Look at Dr. Zhang's files in directory:
/Users/kunzhang/WorkSpace/Exome/Solexa/NA12878
  • Copy best quality one:
cp NA12878_061009_061109_40bp_sequence.qualTrimmed.soapsnp.pileup.variants.txt ~/NA12878/hiqual
  • Separate into hits and misses:
grep 'rs' hiqual>hiqual.in_dbSNP
grep -v 'rs' hiqual>hiqual.no_dbSNP
    • 1963/13394 not in dbSNP (14.7%)
  • Get low-quality SNPs:
cp NA12878_072109_lib_sequence.qualTrimmed.soapsnp.pileup.variants.txt ~/NA12878/loqual
  • 148846/169073 not in dbSNP (88.0%)
  • Similar distributions of quality and coverage of non-dbSNP SNPs

File:Zsnp qr.jpg

Other Parameters[edit]

  • Types of SNPs: edit transition.pl
awk '{print $1 "\t" $2 "\t" $3 "\t" $4}' hiqual.no_dbSNP>hiqual.no_dbSNP.refcall
Scripts/transition.pl hiqual.no_dbSNP.refcall >hiqual.no_dbSNP.trans
  • Results for high-quality SNPs:
' transition transversion
reference 283 315
transition 512 1
transversion 852

File:Hiqual trans.jpg

  • Low-quality SNPs:
' transition transversion
reference 5319 9867
transition 47545 10
transversion 86105

File:Loqual trans.jpg

  • There is a vast over-representation of homozygous SNPs of both types (transition-transition and transversion-transversion) and a slightly higher rate of single-transversion SNPs
  • Plot quality and coverage of SNP types
paste loqual.no_dbSNP.refcall loqual.no_dbSNP.qr.txt >loqual.no_dbSNP.refqr
Scripts/transition.pl loqual.no_dbSNP.refqr >loqual.no_dbSNP.trans 
grep '22$' loqual.no_dbSNP.trans|awk '{print $5"\t"$6}'>loqual.22.txt

File:Trans compare.jpg

    • EDIT 8/4: Legend should have n=852 for 'Hi 22'
  • Are the same trends observable for SNPs at dbSNP locations?

Error Characteristic[edit]

  • Make heatmap for high-quality data set (title is wrong...)

File:Hiqual dbSNP.jpg

  • What if we compare homozygous and heterozygous SNPs?
grep 'het' hiqual.in_dbSNP | awk '{print $6"\t"$8"\t1"}' >hiqual.in_dbSNP.het
mv hiqual.in_dbSNP.het hiqual.in_dbSNP.het.txt
grep 'het' hiqual.no_dbSNP | awk '{print $6"\t"$8"\t0"}' >hiqual.no_dbSNP.het.txt

File:Hiqual het dbSNP.jpg

  • Similarly extract data for homozygous SNPs

File:Hiqual homo dbSNP.jpg