Ylaine/2009-7-29

From ZhangLabWiki
Jump to navigation Jump to search

False Negatives[edit]

  • Combine output files
cat fn.1Ka* > fn.1K
  • Compare to MAQ/SAM SNPs
../Scripts/fnCompare.pl snp.chr.txt fn.1K>fn.1Kcompare.txt
  • Of true SNPs (not equal to reference):
Total: 10151
Miss 1719
Mismatch 1179
Match 7253
  • Get statistics for higher coverage (8x, 20x)
  • Split fn file based on matching into 3 files containing only depth
  • To calculate # of matches with at least 8x coverage, add these two outputs:
grep '^.$' fn.temp.match | grep -c '[89]'
grep -c -v '^.$' fn.temp.match

esult

' Match Miss Mismatch Total
2x 7253 1719 1179 10151
8x 4272 179 57 4508
10x 3552 113 37 3702
20x 1580 27 3 1610

File:Fn maqsam.jpg

False Positives[edit]

  • Gold standard data set: 1000 Genomes data with high coverage and quality (20 and 40)
  • Downloaded output of MAQ SNP filter:
ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/working/20080812_ceu_trio/NA12878/
    • These should all already be different from the reference!!
  • From MAQ manual:

"Each line consists of chromosome, position, reference base, consensus base, Phred-like consensus quality, read depth, the average number of hits of reads covering this position, the highest mapping quality of the reads covering the position, the minimum consensus quality in the 3bp flanking regions at each side of the site (6bp in total), the second best call, log likelihood ratio of the second best and the third best call, and the third best call."

  • Parse out relevant information
awk '{print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6}' NA12878.snp.flt >allsnp.qd
  • Write Perl script to filter out coverage <20 and quality <40, convert SNPs with "oneToThree" hash
../Scripts/qrFilter.pl allsnp.qd >snp.gold
    • 1,431,676 in this set
  • Filter by target region
../Scripts/snpInRange.pl snp.gold>snp.gold.inRange &
    • Takes forever to run!

MAQ/SAM specific[edit]

  • Filter by coverage (compare to pileup file, edit 'snpCovered.pl')
    • Need to also compare to consensus: get pileup with call instead of reference
awk '{print $2 "\t" $8 "\t" $4}' all.pileup>all.pileup.temp &
    • Remove coverage = 1 regions
grep -v '   1       ' all.pileup.temp>pileup.noone.temp
    • Convert to chromosome locations
../Scripts/convertCnsGenome.pl pileup.noone.temp>pileup.call.chr &