Ylaine/2009-8-17

From ZhangLabWiki
Jump to navigation Jump to search

New 1K Genomes Reference[edit]

  • BowtieMaq: Get 5 more matches and 5 fewer misses than when trio file used for comparison
    • Should be more. Too many mismatches? Debugging code.
    • Found bug: using reference instead of 1K call for comparison
../Scripts/hetOnly.pl no_dbSNP.new1K >no_dbSNP.new1K.het
  • Quality and coverage alone are the best combination of parameters
  • For n=50, p=0.3, f = 0.75:
' Mean Std
FP 25.26 7.32
MD 15.20 6.17
Overall 18.32 4.02
  • Try with SOAP/SAM
../Scripts/compare1K.pl hiqual.no_dbSNP.het>hiqual.no_dbSNP.het.1K
    • Place in separate folder

New Data[edit]

  • SNPs thresholded on 8x coverage and 30x quality

Part 1[edit]

"What is the total size of exonic region that were covered by >=8x in this data set?"

WorkSpace/Exome/Solexa/NA12878/NA12878_061009_061109_080509_081309_40bp_sequence.bowtie.pileup
  • Re-named 'orig.pileup'
genome-tech:Aug17Data ygerardin$ ../Scripts/threshold.pl orig.pileup 8 >cov8.temp&
wc -l cov8.temp
  • ANSWER: 23,830,845

Part 2[edit]

"How many variants found are within these regions and how many located outside of our target regions (non-specific capture)?"

  • Use 'calculateSNPeff.pl', which outputs the number of probes covering the region as well as median efficiency of the probes
../Scripts/calculateSNPeff.pl orig.snp >snp.eff &
  • ANSWER: 12412/20761 or 59.79% (for coverage >= 8 there are 10656/14708 or 72.45%)
    • EDIT 8/18: Ran different script (snpTargeted.pl) and got much higher percentage.

Part 5[edit]

"Using 1KG data as reference, what is our false positive rate as the function of the filter we used? This could be plotted as a curve. Also what is the total variants we can call as the function of false-positive rate?"

  • Compare to 1KG data
../Scripts/compare1K.pl orig.snp>snp.1k
  • Separate into matches and misses. There are 74 heterozygous and 522 homozygous mismatches
grep 'match' snp.1k |grep 'het'>1k.match.het.txt
grep 'match' snp.1k |grep 'homo'>1k.match.homo.txt
grep 'miss' snp.1k |grep 'homo'>1k.miss.homo.txt
grep 'miss' snp.1k |grep 'het'>1k.miss.het.txt
    • EDIT 8/18: mismatches put into match file
  • Use [pPresent, numPresent] output of Matlab script 'snpCurve.m'
  • False positive rate as function of filter:

File:Aug17 percent1k.jpg

    • Lower left hand corner:

File:Aug17 percent1k closeup.jpg

  • ROC (each "stripe" is either a single quality or coverage threshold)

File:Aug17 roc.jpg File:Aug17 roc2.jpg

Other[edit]

File:Aug17 match1k.jpg

Thoughts about classification[edit]

  • We need a universal classification scheme that doesn't rely on knowing the classes of the training set.
  • A more analogous problem is one of clustering