heterogeneous サンプルのDNA sequencerデータから mutationを調べる

例えば、データベース上の配列として
AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT
こんな配列があったとして

シークエンス結果から1st peakと2nd peakを抽出すると

AACTCTGTTGTACTCTAATTCTTACCTCCGTCGATTTATT
AACTCTGTTGTACTCTgggaaaTcatTCgcctagcTaAac

こんなだったとする。２つのシークエンスが混ざり合っているのだが、どちらか一方または両方に何らかのdeletionやinsertionが入っているためにこのようなシークエンス結果となっている。
mutationの正解は

 ref: AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT
chr1: AACTCTGTTGTACTCTAGGTCATAATTCCCTCGGCTAATT
                      ^ A insertion
chr2: AACTCTGTTGTACTCTG   ATAATTCCCTCGGCTAATTTAAC
                       ^^^ 3 nt deletion

こういうふうになっているのだけれど、これを機械的に探し出す方法を考える。

$ python3
Python 3.6.5 (default, Apr 25 2018, 14:26:36) 
[GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import regex
>>> w = 'AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT'
>>> s1 = 'AACTCTGTTGTACTCTAATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'AACTCTGTTGTACTCTGGGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> m = regex.finditer(s,w,overlapped=True)
>>> for match in m:
...     print(match.start(), match.group(), match.end())
... 
>>> 

最もシンプルな方法はこのようにmatchを使って一致する位置を探してやるのだけれど、これはどちらか一方がWTであったときにしか成功しない。

>>> import regex
>>> w = 'AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT'
>>> s1 = 'ATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'GGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> m = regex.finditer(s,w,overlapped=True)
>>> for match in m:
...     print(match.start(), match.group(), match.end())
... 
16 GGTCATAATTCCCTCGGCTAATT 39
20 ATAATTCCCTCGGCTAATTTAAC 43
>>> 

このように余分な配列をうまく切り詰めるときっちりとmatchする部分が出てくるので、順次配列を切り詰めていってやればいいのだけれど、手動でやるのは面倒くさい。
そこで自動でスキャンすることにした。

>>> import regex
>>> w = 'AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT'
>>> s1 = 'AACTCTGTTGTACTCTAATTCTTACCTCCGTCGATTTATT'
>>> s2 = 'AACTCTGTTGTACTCTGGGAAATCATTCGCCTAGCTAAAC'
>>> s = "".join('['+i+j+']' for i, j in zip(s1, s2))
>>> n = 1
>>> while s:
...     m = regex.finditer(s,w,overlapped=True)
...     print(str(n)+': '+s)
...     for match in m:
...             print(match.start(), match.group(), match.end())
...     s = s[4:]
...     n += 1
... 
1: [AA][AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
2: [AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
3: [CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
4: [TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
5: [CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
6: [TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
7: [GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
8: [TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
9: [TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
10: [GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
11: [TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
12: [AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
13: [CC][TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
14: [TT][CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
15: [CC][TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
16: [TT][AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
17: [AG][AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
18: [AG][TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
16 GGTCATAATTCCCTCGGCTAATT 39
20 ATAATTCCCTCGGCTAATTTAAC 43
19: [TG][TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
17 GTCATAATTCCCTCGGCTAATT 39
21 TAATTCCCTCGGCTAATTTAAC 43
20: [TA][CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
18 TCATAATTCCCTCGGCTAATT 39
22 AATTCCCTCGGCTAATTTAAC 43
21: [CA][TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
19 CATAATTCCCTCGGCTAATT 39
23 ATTCCCTCGGCTAATTTAAC 43
22: [TA][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
20 ATAATTCCCTCGGCTAATT 39
24 TTCCCTCGGCTAATTTAAC 43
23: [TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
21 TAATTCCCTCGGCTAATT 39
25 TCCCTCGGCTAATTTAAC 43
24: [AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
22 AATTCCCTCGGCTAATT 39
26 CCCTCGGCTAATTTAAC 43
25: [CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
23 ATTCCCTCGGCTAATT 39
27 CCTCGGCTAATTTAAC 43
26: [CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
24 TTCCCTCGGCTAATT 39
28 CTCGGCTAATTTAAC 43
27: [TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
25 TCCCTCGGCTAATT 39
29 TCGGCTAATTTAAC 43
28: [CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
26 CCCTCGGCTAATT 39
30 CGGCTAATTTAAC 43
29: [CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
27 CCTCGGCTAATT 39
31 GGCTAATTTAAC 43
30: [GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
28 CTCGGCTAATT 39
32 GCTAATTTAAC 43
31: [TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]
29 TCGGCTAATT 39
33 CTAATTTAAC 43
32: [CT][GA][AG][TC][TT][TA][AA][TA][TC]
30 CGGCTAATT 39
34 TAATTTAAC 43
33: [GA][AG][TC][TT][TA][AA][TA][TC]
31 GGCTAATT 39
35 AATTTAAC 43
34: [AG][TC][TT][TA][AA][TA][TC]
32 GCTAATT 39
36 ATTTAAC 43
35: [TC][TT][TA][AA][TA][TC]
33 CTAATT 39
37 TTTAAC 43
36: [TT][TA][AA][TA][TC]
21 TAATT 26
34 TAATT 39
38 TTAAC 43
37: [TA][AA][TA][TC]
21 TAAT 25
22 AATT 26
34 TAAT 38
35 AATT 39
39 TAAC 43
38: [AA][TA][TC]
0 AAC 3
22 AAT 25
23 ATT 26
35 AAT 38
36 ATT 39
40 AAC 43
39: [TA][TC]
1 AC 3
3 TC 5
7 TT 9
11 AC 13
13 TC 15
18 TC 20
20 AT 22
23 AT 25
24 TT 26
25 TC 27
29 TC 31
36 AT 38
37 TT 39
38 TT 40
41 AC 43
40: [TC]
2 C 3
3 T 4
4 C 5
5 T 6
7 T 8
8 T 9
10 T 11
12 C 13
13 T 14
14 C 15
15 T 16
18 T 19
19 C 20
21 T 22
24 T 25
25 T 26
26 C 27
27 C 28
28 C 29
29 T 30
30 C 31
33 C 34
34 T 35
37 T 38
38 T 39
39 T 40
42 C 43
43 C 44
44 C 45
46 C 47
47 T 48

頭から17塩基分除いてやるとrefの16番目と20番目にずれてマッチするよ、という結果が読み取れる。
まあここからインターフェースを使いやすくいろいろ改善してやり、ついでに1st peakと2nd peakの抽出も自動で.ab1ファイルから行えるようにつなげてやると、使いやすくなるはず。現状前側から削っていってマッチする場所を探す事になっているが、後ろからも削ってやらないとうまくマッチしない場合もある。