Matt:LabNotes/2013-8-3

From ZhangLabWiki
Jump to navigation Jump to search

Quantifying Errors in CA12k and Agi26k Oligo Pools[edit]

Plan outline[edit]

  • Going to use combination of CIGAR string 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
  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_quantifyOligo.pl into 7 categories: perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ mismatches

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_quantifyOligo.pl

Result stats:

Total reads Filtered reads Perfect match 1 ins/del 1 sub 2 ins/del 2 sub 1 ins/del & 1 sub 3+
5195976 4533981 (87.26% of total) 2830659 (62.4%) 69425 (1.5%) 1170878 (25.8%) 6148 (0.14%) 250843 (5.5%) 31985 (0.7%) 174043 (3.8%)
  • 11 Probes had 0 reads pass through filter
  • 22 Probes that passed through filter had 0 perfect reads
  • 12,322 Probes had at least 1 perfect read (99.73%)

Validate counts with example[edit]

ILLUMINA-379DBF:2:25:10040:14820#0 0 Probe4075 1 42 50M * 0 0 =========T========================================  ;+?<.::><4B?C?BB<DBD?=;;=DDDDBDBBBB9>9C>D@@DDDBD@2 AS:i:253 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:9G40 YT:Z:UU
ILLUMINA-379DBF:2:61:4007:10295#0 0 Probe4075 1 40 6M1I43M * 0 0 ======T=========================================== IIIIEIIIIIIGIIIDIIIIIIIIIIBIHHEEIGIGIHIHHIFIIGFBFI AS:i:248 XN:i:0 XM:i:0 XO:i:1 XG:i:1 NM:i:1 MD:Z:49 YT:Z:UU
ILLUMINA-379DBF:2:81:4293:19366#0 0 Probe4075 1 42 50M * 0 0 =======A========================================== GGGGGGGGGGDGGG>GGGGGBDGGGHGGHHHHHGDDG>BGGBGGGA>)C# AS:i:251 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:7C42 YT:Z:UU

Probe ID Total reads Filtered reads Perfect match 1 ins/del 1 sub 2 ins/del 2 sub 1 ins/del & 1 sub 3+
Probe:4075 4 3 0 1 2 0 0 0 0

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_quantifyOligo.pl (same as for CA12k above but different file I/O)

Result stats:

Total reads Filtered reads Perfect match 1 ins/del 1 sub 2 ins/del 2 sub 1 ins/del & 1 sub 3+
4953642 4448339 (89.8%) 3037261 (68.3%) 88179 (2%) 966306 (21.7%) 11836 (0.27%) 181135 (4.1%) 29938 (0.67%) 133684 (3%)
  • 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)
  • 12,957 Probes had at least 1 perfect read (99.95%)

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_quantifyOligo.pl (same as for CA12k above but different file I/O)

Result stats:

Total reads Filtered reads Perfect match 1 ins/del 1 sub 2 ins/del 2 sub 1 ins/del & 1 sub 3+
4052727 3578431 (88.3% of total) 2453051 (68.6%) 84448 (2.4%) 756094 (21.1%) 12492 (0.35%) 135634 (3.8%) 27104 (0.76%) 109608 (3%)
  • 4 Probes had 0 reads pass through filter
  • 17 Probes that passed through filter had 0 perfect reads
  • 13,158 Probes had at least 1 perfect read (99.84%)