Matt:LabNotes/2013-8-9
Jump to navigation
Jump to search
Quantifying Errors in CA12k and Agi26k Oligo Pools (Ignoring low base quality substitution errors)[edit]
- For previous quantifying errors without considering base quality, see: http://genome-tech.ucsd.edu/LabNotes/index.php/Matt:LabNotes/2013-8-3
Plan outline[edit]
- Going to use combination of CIGAR string, baseQuality, and readBases with matching bases replaced with '=' (using samtools calmd -e) because easier to parse
- Could parse MD tag in sam file to count ins/del/subs of each mapped read
- Could use NM (edit distance to the reference) to see how many changes necessary to make read equal to reference
- Could also use mpileup
- Going to convert CIGAR string into a string with length equal to readBases (example: 5M1I6M -> MMMMMIMMMMMM) to determine if readBases letter is a insertion or substitution
- Any base index of a read with substitution error can be checked with the base quality scores to see if quality score is >= 20 (Phred33 scale)
- Turn sam -> bam
- Sort bam file
- Index sorted bam file
- Filter out unmapped reads and mapped reads with MAPQ < 8 (16% alignment is wrong - Understanding Mapping Quality)
- Remove ':' from probe names in sam file and reference fasta file (Because ':' is used by samtools to signify end of contig name)
- Replace any matched bases in readBases line to '='
- Count ins/del/subs with Count_mismatches_from_SAMwEquals_quantifyOligov2.pl into 7 categories: perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ mismatches
- Count up every error and divide by total number of bases: Error rate is output to command line
Count mismatches of CA12k[edit]
samtools view -bS Readsalign2probes_CA_default.txt | samtools sort - CA_default_sorted samtools view -h -F 4 -q 8 CA_default_sorted.bam > CA_default_sorted_filtered.sam
perl -pi.back -e 's/Probe:/Probe/g;' CA_default_sorted_filtered.sam perl -pi.back -e 's/Probe:/Probe/g;' CAprobes_to_order_50bp.fa
samtools view -bS CA_default_sorted_filtered.sam > CA_default_sorted_filtered.bam
samtools calmd -eS CA_default_sorted_filtered.sam CAprobes_to_order_50bp.fa > CA_default_sf=.sam
perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl
Error Rate: 0.00785889045410645 = 0.79%
Error Rate of Insertions: 0.000640143838273694 = 0.064%
Error Rate of Deletions: 2.23512184987101e-05 = 0.0022%
Error Rate of Substitutions: 0.00719639539733404 = 0.72%
Result stats:
Probe ID | Total reads | Filtered reads | Perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ |
Sums: | 5195976 | 4533981 (87.26%) | 3176472 (70.1%) | 77874 (1.7%) | 962022 (21.2%) | 4934 (0.11%) | 178599 (3.9%) | 30019 (0.66%) | 104061 (2.3%) |
- 11 Probes had 0 reads pass through filter
- 27 Probes had 0 perfect reads
- 12,328 Probes had at least 1 perfect read (99.78%)
Count mismatches of Agi26k0gap[edit]
samtools view -bS Readsalign2split_revcomp_0gap_default.txt > Agi26k0gap_default.bam samtools sort Agi26k0gap_default.bam Agi26k0gap_default_sorted samtools index Agi26k0gap_default_sorted.bam Agi26k0gap_default_sorted.bai samtools view -h -F 4 -q 8 Agi26k0gap_default_sorted.bam > Agi26k0gap_default_sorted_filtered.sam
perl -pi.back -e 's/Probe:/Probe/g;' Agi26k0gap_default_sorted_filtered.sam perl -pi.back -e 's/Probe:/Probe/g;' Agi26k0gap_corrected.fa
samtools calmd -eS Agi26k0gap_default_sorted_filtered.sam Agi26k0gap_corrected.fa > Agi26k0gap_default_sf=.sam
perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)
Error Rate: 0.00662328568034046 = 0.66%
Error Rate of Insertions: 0.00079017808669708 = 0.079%
Error Rate of Deletions: 1.83753980980316e-05 = 0.0018%
Error Rate of Substitutions: 0.00581473219554535 = 0.58%
Result stats:
Probe ID | Total reads | Filtered reads | Perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ |
Sums: | 4953642 | 4448339 (89.8%) | 3294962 (74.1%) | 102309 (2.3%) | 808922 (18.2%) | 12240 (0.28%) | 129671 (2.9%) | 24455 (0.55%) | 75780 (1.7%) |
- All probes had at least 5 reads pass through filter
- 1 probe had 0 perfect reads (all single substitution error)
- 12,963 Probes had at least 1 perfect read (99.99%)
Count mismatches of Agi26k20gap[edit]
samtools view -bS Readsalign2split_revcomp_20gap_default.txt > Agi26k20gap_default.bam samtools sort Agi26k20gap_default.bam Agi26k20gap_default_sorted samtools index Agi26k20gap_default_sorted.bam Agi26k20gap_default_sorted.bai samtools view -h -F 4 -q 8 Agi26k20gap_default_sorted.bam > Agi26k20gap_default_sorted_filtered.sam
perl -pi.back -e 's/Probe:/Probe/g;' Agi26k20gap_default_sorted_filtered.sam perl -pi.back -e 's/Probe:/Probe/g;' Agi26k20gap_corrected.fa
samtools calmd -eS Agi26k20gap_default_sorted_filtered.sam Agi26k20gap_corrected.fa > Agi26k20gap_default_sf=.sam
perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)
Error Rate: 0.00636130751158818 = 0.64%
Error Rate of Insertions: 0.000959521086196716 = 0.096%
Error Rate of Deletions: 1.98355089143817e-05 = 0.0020%
Error Rate of Substitutions: 0.00538195091647708 = 0.54%
Result stats:
Probe ID | Total reads | Filtered reads | Perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ |
Sums: | 4052727 | 3578431 (88.3%) | 2676974 (74.8%) | 101104 (2.8%) | 619290 (17.3%) | 14445 (0.40%) | 92296 (2.6%) | 20184 (0.56%) | 54138 (1.5%) |
- 4 probes had 0 reads pass through filter
- 13 probes had 0 perfect reads
- 13,166 Probes had at least 1 perfect read (99.9%)
Error Rate Comparison Summary[edit]
' | CA12k | Agi26k0gap | Agi26k20gap |
Overall Error Rate | 0.0078 | 0.0066 | 0.0064 |
Rate of Insertions | 0.00064 | 0.00079 | 0.00096 |
Rate of Deletions | 0.000022 | 0.000018 | 0.00002 |
Rate of Substitutions | 0.0072 | 0.0058 | 0.0054 |