False Negatives[edit]
cat fn.1Ka* > fn.1K
../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
../Scripts/snpInRange.pl snp.gold>snp.gold.inRange &
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 &