以前、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文字分ずつ削っていって、合致する配列が現れるポイントを探る感じにしてみた。
下の方の数文字しかないところは全く意味がないので削る文字数の制限とかつけるともっと見やすい感じになるかも。
というわけでかなり作り込みました。
これはかなりいい感じじゃなかろうか