CRISPRで編集がかかったゲノムシークエンスを読むプログラム

CRISPRで編集がかかったと思われるT0サンプルのゲノムをシークエンスすると、変異は２本の染色体にランダムに起こるため、基本的にはヘテロとなってシークエンスデータの波形が２つ分重なって検出されることになるのは先日に書いたとおり。
これをパズルのように分解して再構築して２本の染色体に分離し、それぞれの変異を同定しているわけだが、なんとか自動できないものかと。
というわけでプログラムを作ってみた。

ref.txt:入力としてreference配列のファイル
motif.txt:とりあえず２つのシークエンスが混ざっているのを記述したファイル
例えば
ref.txt
ATGGTGAGCATCAGCCTCTGGGTGG

ベースコールでは
ATGGTGAGCATGACTCTCTGTGG
となっていて、ATGGTGAGCまでは山がずれてなくて、以後２番手の山がかぶっている
ATGGTGAGCCACCGCCTGGGGGT
のように読み取れたとするなら
motif.txtに
ずれ始めから
A/CT/AG/CA/CC/GT/CCTC/GT/GGT/GGG/T
というふうに２つが重なっている部分は/で２つを併記する。

これを
crisgt.py

import sys
import argparse
import re

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set reference seq file", type=str)
        parser.add_argument("motif", help="please set search motif file", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        w = list1.read().replace('a','A').replace('c','C').replace('g','G').replace('t','T')
        with open(motif) as list2:
            k = list2.read().strip()
            k = k.replace('a','A').replace('c','C').replace('g','G').replace('t','T').replace('N', '[ACGT]').replace('n', '[ACGT]').replace('A/T','[AT]').replace('A/C','[AC]').replace('A/G','[AG]').replace('C/G','[CG]').replace('C/T','[CT]').replace('G/T','[GT]')
            m = re.finditer(k,w)
            for match in m:
               print(match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()

$ python3 crisgt.py ref.txt motif.txt
ATCAGCCTCTGGGT 9 23

というふうに処理すると
あれれ。２つ目の位置が記録されないぞ？

qiita.com

お、これか。pythonのmatchは重複を認めない設定になっているのか。こりゃ以前に作ったプログラムも見直さないといかんな。

            m = re.finditer(k,w)

import regex
            m = regex.finditer(k,w,overlapped=True)

これで解決。

$ python3 crisgtype.py ref.txt motif.txt 
ATCAGCCTCTGGGT 9 23
CAGCCTCTGGGTGG 11 25

これにより、９（10塩基め）から正常な配列に対し、2塩基抜けて12塩基目がずれて来ていることがわかる。

あとは重なった波形から自動的に2つのピークを抽出してきたりしてくれるとすごくいいが、そこはちょっとむずかしいかも。

import sys
import argparse
import re
import regex

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set reference seq file", type=str)
        parser.add_argument("seq1", help="please set seq1 file", type=str)
        parser.add_argument("seq2", help="please set seq2 file", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'seq1'):
        seq1 = args.seq1
    if hasattr(args, 'seq2'):
        seq2 = args.seq2

    with open(file) as list1:
        w = list1.read().replace('a','A').replace('c','C').replace('g','G').replace('t','T')
        with open(seq1) as s1:
            sq1 = s1.read().strip()
        with open(seq2) as s2:
            sq2 = s2.read().strip()
        s = "".join('['+i+j+']' for i, j in zip(sq1, sq2))
        m = regex.finditer(s,w,overlapped=True)

        for match in m:
            print(match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()

お、ずいぶんスッキリしてきたぞ。
山が2つ重なりだすところからのベースコールデータと、2番手の山を読み取った配列をそれぞれseq1.txt、seq2.txtとしておいてこれを

$ python3 crisgt.py ref.txt seq1.txt seq2.txt
ATCAGCCTCTGGGT 9 23
CAGCCTCTGGGTGG 11 25

こりゃもう実用レベルじゃね？