CRseekを用いたCRISPRターゲットの自動予測

投稿日: 投稿者:

CRseekとは

  • CRseekとはCRISPRのデザイン用に開発されたPythonのモジュールです。
  • ４つの一般的な切断予測アルゴリズムが実装され、ドラフト品質のゲノムでも動作するように適合化されています。
  • Biopythonのパッケージと完全に統合されているため、遺伝子編集の前後の配列の操作が簡単にできます。

CRseekを用いた配列予測

fasta形式用のgenomeとgeneのファイルパスを指定することで、テキスト形式で配列を５つ返してくれるスクリプトを作成しました。

インストール方法

  • 動作環境
    • Python(anaconda3-5.3.0(python3.6))
    • 以下コマンドによってインストールを行います
conda config --add channels conda-forge
conda config --add channels defaults
conda config --add channels r
conda config --add channels bioconda
conda create -yq -n crseek-environment --file requirements.conda
source activate crseek-environment
pip install -q -r requirements.pip
python setup.py install

ターゲット予測

  • 以下、ターゲット予測用のスクリプトです。crseek_test.pyというファイル名で保存してください。
  • genome , geneともにfasta形式を想定しています。
  • python crseek_test.py gene_path genome_pathで動きます。
  • pyyamlのWarning(YAMLLoadWarning)を吐きますが、問題なく動きます。
  • 実行したディレクトリにoutput.txtという名前のファイルが作成されます。
import sys
import crseek
from crseek import utils
from crseek import estimators
from crseek import annotators

from Bio import SeqIO
from Bio import Seq
import gzip
import pandas as pd
import numpy as np

import warnings
warnings.simplefilter('ignore', category=Seq.BiopythonWarning)

gene_path = sys.argv[1]
genome_path = sys.argv[2]

# Load  genome Sequence
with open(genome_path) as handle:
    genome = list(SeqIO.parse(handle, 'fasta'))[0]

# Load gene Sequence
with open(gene_path) as handle:
    gene = list(SeqIO.parse(handle, 'fasta'))[0]

def main(gene, genome):
    possible_targets = utils.extract_possible_targets(gene, pams=('NGG'))
    print('%i possible targets' % len(possible_targets))

    estimator = estimators.CFDEstimator.build_pipeline()

    # cas-offinder is used to rapidly find all positions in the Thermofilum genome that are <=5 mismatches for each of the potential gRNAs.
    possible_binding = utils.cas_offinder(possible_targets, 5, locus=[genome])
    possible_binding["Score"] = estimator.predict_proba(
        possible_binding.values)

    results = possible_binding.groupby('spacer')['Score'].agg('max')
    results.sort_values(inplace=True)
    results = results.head().index

    with open('output.txt', 'w') as f:
        for i in range(5):
            DNA = results[i].back_transcribe()
            f.write(str(DNA) + '\n')

if __name__ == "__main__":
    main(gene, genome)

今後やってみたいこと

  • ターゲット配列をゲノムにマッピングする。
  • ターゲットがエクソン間を跨いでいる配列などは除外する。
  • エクソンの端の方の配列は除外する。