Matt:LabNotes/2013-8-9

From ZhangLabWiki
Jump to navigation Jump to search

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