CRseekとは
- CRseekとはCRISPRのデザイン用に開発されたPythonのモジュールです。
- 4つの一般的な切断予測アルゴリズムが実装され、ドラフト品質のゲノムでも動作するように適合化されています。
- Biopythonのパッケージと完全に統合されているため、遺伝子編集の前後の配列の操作が簡単にできます。
CRseekを用いた配列予測
fasta形式用のgenomeとgeneのファイルパスを指定することで、テキスト形式で配列を5つ返してくれるスクリプトを作成しました。
インストール方法
- 動作環境
- 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
ともにfasta形式を想定しています。genome , gene
で動きます。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)
今後やってみたいこと
- ターゲット配列をゲノムにマッピングする。
- ターゲットがエクソン間を跨いでいる配列などは除外する。
- エクソンの端の方の配列は除外する。