Matt:LabNotes/2013-8-3
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
- 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_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%)