以前、CRISPRによって編集された遺伝子配列を解析するプログラムを書いたが、100％マッチするときしか検出できないのは不便だな、ということで、アライメントを取るツールを利用できないだろうかと考え中。
pairwise2 | BioPython の pairwise2 ライブラリーを利用したペアワイズアラインメント
こちらの情報を参考に、biopythonのpairwise2を使ってみる。
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.
>>> from Bio import pairwise2
>>> from Bio.pairwise2 import format_alignment
>>> for a in pairwise2.align.localms("ACCGTCC", "ACG", 2, -1, -0.5, -0.1):
... print(format_alignment(*a))
...
ACCGTCC
| ||
A-CG---
Score=5.5
ACCGTCC
|| |
AC-G---
Score=5.5なるほど。これで比べる2つ目の配列をregexで混合した配列から使えるとよいのだけれど
例えば
AACTCTGTTGTACTCTGGTCATAATTCCCTCGGCTAATTTAACCCGCT
こんな配列があったとして
シークエンス結果が
AACTCTGTTGTACTCTaatTCtTAccTCCgTCGatTtATT AACTCTGTTGTACTCTGGGaaAtcATtcgCctaGCtAaac
こんなだったとする。
正解は
AACTCTGTTGTACTCTAGGTCATAATTCCCTCGGCTAATT
^
AACTCTGTTGTACTCTG ATAATTCCCTCGGCTAATTTAACこういうふうになっているのだけれど
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())
...
>>> こんなふうに以前作ったアプリではうまくヒットしない。一方はミスマッチのinsertionがあり、もう一方はそれとずれてdeletionが入っているからだ
とはいえpairwise2になんとなく組み込んでみても
>>> for a in pairwise2.align.localms(w, s, 2, -1, -0.5, -0.1):
... print(format_alignment(*a))
...
-----AA--C---T---C---T---G---TT------G---T---A---C---T---C---T----G---G--T-------CA--TAA--TT---C--C---CT----------CG--GC--T--------A--A-------TT--TA--A--------CCCGCT
|| | | | | | || | | | | | | | | | | || | | || | | || || || | | | || || | |
[AA][AA][CC][TT][CC][TT][GG][TT][TT][GG][TT][AA][CC][TT][CC][TT][AG][AG][TG][TA][CA][T-A][TT][AC][CA][CT][TT][CC][CG][GC][TC][CT][GA][AG][TC][TT][TA][AA][TA][TC]----
Score=59.9
................これじゃだめだね。
もうちょっと考えよう。
ちなみに
>>> 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
>>> こういうふうに投げる配列をうまく調整すると答えにだいたいたどり着けるのだけど、どれくらいの位置にどれくらいdeletionが入っているのかというのが運次第というのではちょっとね。
>>> 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というわけでちょっと改良して、前から1文字分ずつ削っていって、合致する配列が現れるポイントを探る感じにしてみた。
下の方の数文字しかないところは全く意味がないので削る文字数の制限とかつけるともっと見やすい感じになるかも。
というわけでかなり作り込みました。
これはかなりいい感じじゃなかろうか
