Matt:LabNotes/2013-8-9: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Mzcai
(Created page with "==Quantifying Errors in CA12k and Agi26k Oligo Pools (Ignoring low base quality substitution errors)== ===Plan outline=== *Going to use combination of CIGAR string, baseQuali...")
 
>Mzcai
No edit summary
 
(5 intermediate revisions by the same user not shown)
Line 1: Line 1:
==Quantifying Errors in CA12k and Agi26k Oligo Pools (Ignoring low base quality substitution errors)==
==Quantifying Errors in CA12k and Agi26k Oligo Pools (Ignoring low base quality substitution errors)==
 
*For previous quantifying errors without considering base quality, see: http://genome-tech.ucsd.edu/LabNotes/index.php/Matt:LabNotes/2013-8-3
===Plan outline===
===Plan outline===
*Going to use combination of CIGAR string, baseQuality, and readBases with matching bases replaced with '=' (using samtools calmd -e) because easier to parse
*Going to use combination of CIGAR string, baseQuality, and readBases with matching bases replaced with '=' (using samtools calmd -e) because easier to parse
Line 32: Line 32:
   perl [[Media:Count_mismatches_from_SAMwEquals_quantifyOligov2.txt | Count_mismatches_from_SAMwEquals_quantifyOligov2.pl]]
   perl [[Media:Count_mismatches_from_SAMwEquals_quantifyOligov2.txt | Count_mismatches_from_SAMwEquals_quantifyOligov2.pl]]


Error Rate: 0.785889045410645%<br>
Error Rate: 0.00785889045410645 = 0.79%<br>
Error Rate of Insertions: 0.000640143838273694 = 0.064%<br>
Error Rate of Deletions: 2.23512184987101e-05 = 0.0022%<br>
Error Rate of Substitutions: 0.00719639539733404 = 0.72%<br>
Result stats:<br>
Result stats:<br>
{| {{table}}
{| {{table}}
Line 49: Line 52:
|}
|}
*11 Probes had 0 reads pass through filter
*11 Probes had 0 reads pass through filter
*22 Probes that passed through filter had 0 perfect reads
*27 Probes had 0 perfect reads
*12,322 Probes had at least 1 perfect read (99.73%)
*12,328 Probes had at least 1 perfect read (99.78%)


===Count mismatches of Agi26k0gap===
===Count mismatches of Agi26k0gap===
Line 64: Line 67:


   perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)
   perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)
Error Rate: 0.662328568034046%<br>
Error Rate: 0.00662328568034046 = 0.66%<br>
Error Rate of Insertions: 0.00079017808669708 = 0.079%<br>
Error Rate of Deletions: 1.83753980980316e-05 = 0.0018%<br>
Error Rate of Substitutions: 0.00581473219554535 = 0.58%<br>
Result stats:<br>
Result stats:<br>
{| {{table}}
{| {{table}}
Line 81: Line 87:
|}
|}
*All probes had at least 5 reads pass through filter
*All probes had at least 5 reads pass through filter
*7 Probes had 0 perfect reads (almost all of which were due to substitution errors)
*1 probe had 0 perfect reads (all single substitution error)
*12,957 Probes had at least 1 perfect read (99.95%)
*12,963 Probes had at least 1 perfect read (99.99%)
 
===Count mismatches of Agi26k20gap===
===Count mismatches of Agi26k20gap===
   samtools view -bS Readsalign2split_revcomp_20gap_default.txt > Agi26k20gap_default.bam
   samtools view -bS Readsalign2split_revcomp_20gap_default.txt > Agi26k20gap_default.bam
Line 96: Line 103:
   perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)
   perl Count_mismatches_from_SAMwEquals_quantifyOligov2.pl (same as for CA12k above but different file I/O)


Error Rate: 0.636130751158818%<br>
Error Rate: 0.00636130751158818 = 0.64%<br>
Error Rate of Insertions: 0.000959521086196716 = 0.096%<br>
Error Rate of Deletions: 1.98355089143817e-05 = 0.0020%<br>
Error Rate of Substitutions: 0.00538195091647708 = 0.54%<br>
Result stats:<br>
Result stats:<br>
{| {{table}}
{| {{table}}
Line 112: Line 122:
| 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%)
| 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
*4 probes had 0 reads pass through filter
*17 Probes that passed through filter had 0 perfect reads
*13 probes had 0 perfect reads
*13,158 Probes had at least 1 perfect read (99.84%)
*13,166 Probes had at least 1 perfect read (99.9%)
 
==Error Rate Comparison Summary==
{| {{table}}
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''CA12k'''
| align="center" style="background:#f0f0f0;"|'''Agi26k0gap'''
| align="center" style="background:#f0f0f0;"|'''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
|}

Latest revision as of 21:38, 22 August 2013

Quantifying Errors in CA12k and Agi26k Oligo Pools (Ignoring low base quality substitution errors)[edit]

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)
  1. Turn sam -> bam
  2. Sort bam file
  3. Index sorted bam file
  4. Filter out unmapped reads and mapped reads with MAPQ < 8 (16% alignment is wrong - Understanding Mapping Quality)
  5. Remove ':' from probe names in sam file and reference fasta file (Because ':' is used by samtools to signify end of contig name)
  6. Replace any matched bases in readBases line to '='
  7. 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
  8. 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