---
# System prepended metadata

title: CRISPRdirect自動化

---

# CRISPRdirect自動化

[CRISPRdirect](https://crispr.dbcls.jp) には結果をテキスト形式で書き出す機能があり、`curl`や`wget`などのコマンドラインツールと組み合わせれば、ガイドRNAの設計を自動化することができる。

### コマンドラインからCRISPRdirectを使ってみる

例）マウスPzp遺伝子（NM_007376.4）の塩基配列を`NM_007376.4.fa`という名前のFASTAファイルに保存しておく。CRISPRdirectでマウスゲノムmm39を選択し、NM_007376.4を標的とするガイドRNAを設計するには、次のコマンドを実行すればよい。

なお、`@NM_007376.4.fa`のところは`@path/to/NM_007376.4.fa`のようにパスを指定できる。

```
% curl -X POST -F db=mm39 -F format=txt -F upload=@NM_007376.4.fa https://crispr.dbcls.jp/
```

読みにくいので、バックスラッシュ`\`を使い改行を入れてコマンドを実行することにする。

```
% curl -X POST \
    -F db=mm39 \
    -F format=txt \
    -F upload=@NM_007376.4.fa \
    https://crispr.dbcls.jp/
```

結果 ==**(A)**==：

```
# [ CRISPRdirect | 2023-04-14 11:37:47 ]
# sequence_name:	NM_007376.4 Mus musculus PZP, alpha-2-macroglobulin like (Pzp), mRNA
# pam_sequence:	NGG
# specificity_check:	Mouse (Mus musculus) genome, GRCm39/mm39 (Jun, 2020)
# start	end	strand	sequence	GC	Tm	TTTT	RE_sites	hit_20mer	hit_12mer	hit_8mer
#
4	26	+	GGATCAGAGTTCGGGGGCTGAGG	65.00	78.83	0		1	11	14844
5	27	+	GATCAGAGTTCGGGGGCTGAGGG	60.00	76.96	0		1	9	9602
40	62	+	TCTCTGCCCTCTCCACCATGAGG	60.00	77.23	0		1	122	76542
46	68	-	CCCTCTCCACCATGAGGAGAAAC	50.00	71.73	0		1	108	11591
47	69	-	CCTCTCCACCATGAGGAGAAACC	55.00	73.58	0		1	30	9285
52	74	-	CCACCATGAGGAGAAACCAGCTG	55.00	73.68	0	MspA1I,PvuII	1	30	5968
55	77	-	CCATGAGGAGAAACCAGCTGCCC	60.00	78.15	0	MspA1I,PvuII	1	49	7749
68	90	-	CCAGCTGCCCACACCAGCTTTTC	60.00	77.65	0	MspA1I,PvuII	1	58	7374
75	97	-	CCCACACCAGCTTTTCTTTTACT	35.00	66.17	0		1	39	5821
[以下略]
```

### NCBIから塩基配列を取得する

例）マウスPzp遺伝子（NM_007376.4）の塩基配列をFASTA形式で取得するには、

https://www.ncbi.nlm.nih.gov/sviewer/?report=fasta&retmode=text&val=NM_007376.4

にアクセスすればよい。これを`NM_007376.4.fa`という名前のファイルに保存するには、

```
% curl -s 'https://www.ncbi.nlm.nih.gov/sviewer/?report=fasta&retmode=text&val=NM_007376.4' > NM_007376.4.fa
```

または

```
% curl -s \
    -F report=fasta \
    -F retmode=text \
    -F val=NM_007376.4 \
    https://www.ncbi.nlm.nih.gov/sviewer/ \
  > NM_007376.4.fa
```

を実行すればよい。

### 上記の2つのステップを合体

bashやzshのプロセス置換`<(command)`を使って合体させると、

```
% curl -s \
    -X POST \
    -F db=mm39 \
    -F format=txt \
    -F upload=@<( \
       curl -s \
         -F report=fasta \
         -F retmode=text \
         -F val=NM_007376.4 \
         https://www.ncbi.nlm.nih.gov/sviewer/ \
       ) \
    https://crispr.dbcls.jp/ \
  > NM_007376.4.crispr.txt
```

のように実行でき、 ==**(A)**== と同じ結果をワンステップで得ることができる。

### 結果の絞り込み

オフターゲットサイトの少ない「緑のgRNA」だけを抽出するには、第9列と第10列がどちらも`1`であるものを出力すればよい。

```
% cat NM_007376.4.crispr.txt | awk -F "\t" '$9 == 1 && $10 == 1 { print }'
```

結果 ==**(A1)**==：

```
1237	1259	-	CCATTTTTACTACGGATGAGCAT	35.00	65.49	0		1	1	20154
2129	2151	+	GTTTCAACACTATCCGGCAATGG	45.00	69.65	0		1	1	425
2130	2152	+	TTTCAACACTATCCGGCAATGGG	40.00	68.34	0		1	1	267
2133	2155	+	CAACACTATCCGGCAATGGGAGG	55.00	74.31	0		1	1	3868
2210	2232	+	AGCAATGGGCGTGCCGATGATGG	60.00	78.89	0		1	1	557
2211	2233	+	GCAATGGGCGTGCCGATGATGGG	60.00	77.89	0		1	1	397
3193	3215	+	GCAACACTCCGGGAAATACTTGG	50.00	73.03	0		1	1	5255
4391	4413	+	CCTCGGTTTCTCCTTCGCGGTGG	65.00	78.04	0		1	1	184
```

### CDSの判別とexon番号の判別(1)

設計されたガイドRNAの位置がCDS上かどうか、また何番目のexonなのかを調べたい。

例）マウスPzp遺伝子（NM_007376.4）の配列情報をGenBank形式で取得すると、CDSの範囲およびexonの範囲が記載されている。

https://www.ncbi.nlm.nih.gov/sviewer/?report=genbank&retmode=text&val=NM_007376.4


```
     exon            1..142
                     /gene="Pzp"
                     /gene_synonym="A1m; A2m; MAM"
                     /inference="alignment:Splign:2.1.0"
     CDS             57..4544
                     /gene="Pzp"
                     /gene_synonym="A1m; A2m; MAM"
                     /note="alpha 1 macroglobulin; alpha-2-macroglobulin;
                     alpha-2-M; pregnancy zone protein"
                     /codon_start=1
[以下略]
```

このような形で範囲が記載されているので、抜き出して利用する。

まず、`curl`コマンドを用いてNM_007376.4の情報をGenBank形式で取得。

```
% curl -s \
    -F report=genbank \
    -F retmode=text \
    -F val=NM_007376.4 \
    https://www.ncbi.nlm.nih.gov/sviewer/
```

CDSおよびexonの範囲をfeature tableから抜き出しBED形式に整形。

```
% curl -s \
    -F report=genbank \
    -F retmode=text \
    -F val=NM_007376.4 \
    https://www.ncbi.nlm.nih.gov/sviewer/ \
    | perl -ne '/^\s+(exon|CDS)\s+(\d+)\.\.(\d+)/ and \
        print "NM_007376.4\t$2\t$3\t$1\n"'
```

結果：

```
NM_007376.4	1	142	exon
NM_007376.4	57	4544	CDS
NM_007376.4	143	326	exon
NM_007376.4	327	480	exon
NM_007376.4	481	533	exon
NM_007376.4	534	554	exon
NM_007376.4	555	723	exon
NM_007376.4	724	808	exon
NM_007376.4	809	932	exon
NM_007376.4	933	1047	exon
[以下略]
```

exonの番号を付加。これを`annotation.bed`という名前のファイルに保存しておく。

```
% curl -s \
    -F report=genbank \
    -F retmode=text \
    -F val=NM_007376.4 \
    https://www.ncbi.nlm.nih.gov/sviewer/ \
    | perl -ne '/^\s+(exon|CDS)\s+(\d+)\.\.(\d+)/ and \
        print "NM_007376.4\t$2\t$3\t$1\n"' \
    | perl -ne 's/exon$/sprintf("exon_%03d",$count + 1)/e and \
        $count ++ ; print' \
  > annotation.bed
```

結果 ==**(B)**==：

```
NM_007376.4	1	142	exon_001
NM_007376.4	57	4544	CDS
NM_007376.4	143	326	exon_002
NM_007376.4	327	480	exon_003
NM_007376.4	481	533	exon_004
NM_007376.4	534	554	exon_005
NM_007376.4	555	723	exon_006
NM_007376.4	724	808	exon_007
NM_007376.4	809	932	exon_008
NM_007376.4	933	1047	exon_009
[以下略]
```

CDSとexonの情報を、別のファイルに分けて保存しておく。
CDSを含む行を`cds.bed` ==**(B1)**== 、exonを含む行を`exon.bed` ==**(B2)**== として保存した。

```
grep CDS  annotation.bed > cds.bed
grep exon annotation.bed > exon.bed
```

### CDSの判別とexon番号の判別(2)

CRISPRdirectの結果 ==**(A)**== もBED形式にそろえて`crisprdirect.bed`という名前のファイルに保存しておく。

```
% curl -s \
    -X POST \
    -F db=mm39 \
    -F format=txt \
    -F upload=@<( \
       curl -s \
         -F report=fasta \
         -F retmode=text \
         -F val=NM_007376.4 \
         https://www.ncbi.nlm.nih.gov/sviewer/ \
       ) \
    https://crispr.dbcls.jp/ \
  | grep -v '^#' \
  | grep -v '^$' \
  | sed 's/^/NM_007376.4\t/' \
  > crisprdirect.bed
```

結果 ==**(A')**==：

```
NM_007376.4	4	26	+	GGATCAGAGTTCGGGGGCTGAGG	65.00	78.83	0		1	11	14844
NM_007376.4	5	27	+	GATCAGAGTTCGGGGGCTGAGGG	60.00	76.96	0		1	9	9602
NM_007376.4	40	62	+	TCTCTGCCCTCTCCACCATGAGG	60.00	77.23	0		1	122	76542
NM_007376.4	46	68	-	CCCTCTCCACCATGAGGAGAAAC	50.00	71.73	0		1	108	11591
NM_007376.4	47	69	-	CCTCTCCACCATGAGGAGAAACC	55.00	73.58	0		1	30	9285
NM_007376.4	52	74	-	CCACCATGAGGAGAAACCAGCTG	55.00	73.68	0	MspA1I,PvuII	1	30	5968
NM_007376.4	55	77	-	CCATGAGGAGAAACCAGCTGCCC	60.00	78.15	0	MspA1I,PvuII	1	49	7749
NM_007376.4	68	90	-	CCAGCTGCCCACACCAGCTTTTC	60.00	77.65	0	MspA1I,PvuII	1	58	7374
NM_007376.4	75	97	-	CCCACACCAGCTTTTCTTTTACT	35.00	66.17	0		1	39	5821
NM_007376.4	76	98	-	CCACACCAGCTTTTCTTTTACTG	40.00	65.60	0		1	42	7370
[以下略]
```

同様に、緑のgRNAを抽出した結果 ==**(A1)**== もBED形式にして`crisprdirect_green.bed`に保存。

```
% curl -s \
    -X POST \
    -F db=mm39 \
    -F format=txt \
    -F upload=@<( \
       curl -s \
         -F report=fasta \
         -F retmode=text \
         -F val=NM_007376.4 \
         https://www.ncbi.nlm.nih.gov/sviewer/ \
       ) \
    https://crispr.dbcls.jp/ \
  | awk -F "\t" '$9 == 1 && $10 == 1 { print }' \
  | sed 's/^/NM_007376.4\t/' \
  > crisprdirect_green.bed
```

結果 ==**(A1')**==：

```
NM_007376.4	1237	1259	-	CCATTTTTACTACGGATGAGCAT	35.00	65.49	0		1	1	20154
NM_007376.4	2129	2151	+	GTTTCAACACTATCCGGCAATGG	45.00	69.65	0		1	1	425
NM_007376.4	2130	2152	+	TTTCAACACTATCCGGCAATGGG	40.00	68.34	0		1	1	267
NM_007376.4	2133	2155	+	CAACACTATCCGGCAATGGGAGG	55.00	74.31	0		1	1	3868
NM_007376.4	2210	2232	+	AGCAATGGGCGTGCCGATGATGG	60.00	78.89	0		1	1	557
NM_007376.4	2211	2233	+	GCAATGGGCGTGCCGATGATGGG	60.00	77.89	0		1	1	397
NM_007376.4	3193	3215	+	GCAACACTCCGGGAAATACTTGG	50.00	73.03	0		1	1	5255
NM_007376.4	4391	4413	+	CCTCGGTTTCTCCTTCGCGGTGG	65.00	78.04	0		1	1	184
```

### CDSの判別とexon番号の判別(3)

==**(A')**== と ==**(B)**== の重なりを`bedtools intersect`を用いて検出する。

```
% bedtools intersect -a crisprdirect.bed -b annotation.bed -wao
```

結果：

```
NM_007376.4	4	26	+	GGATCAGAGTTCGGGGGCTGAGG	65.00	78.83	0		1	11	14844	NM_007376.4	1	142	exon_001	22
NM_007376.4	5	27	+	GATCAGAGTTCGGGGGCTGAGGG	60.00	76.96	0		1	9	9602	NM_007376.4	1	142	exon_001	22
NM_007376.4	40	62	+	TCTCTGCCCTCTCCACCATGAGG	60.00	77.23	0		1	122	76542	NM_007376.4	1	142	exon_001	22
NM_007376.4	40	62	+	TCTCTGCCCTCTCCACCATGAGG	60.00	77.23	0		1	122	76542	NM_007376.4	57	4544	CDS	5
NM_007376.4	46	68	-	CCCTCTCCACCATGAGGAGAAAC	50.00	71.73	0		1	108	11591	NM_007376.4	1	142	exon_001	22
NM_007376.4	46	68	-	CCCTCTCCACCATGAGGAGAAAC	50.00	71.73	0		1	108	11591	NM_007376.4	57	4544	CDS	11
NM_007376.4	47	69	-	CCTCTCCACCATGAGGAGAAACC	55.00	73.58	0		1	30	9285	NM_007376.4	1	142	exon_001	22
NM_007376.4	47	69	-	CCTCTCCACCATGAGGAGAAACC	55.00	73.58	0		1	30	9285	NM_007376.4	57	4544	CDS	12
NM_007376.4	52	74	-	CCACCATGAGGAGAAACCAGCTG	55.00	73.68	0	MspA1I,PvuII	1	30	5968	NM_007376.4	1	142	exon_001	22
NM_007376.4	52	74	-	CCACCATGAGGAGAAACCAGCTG	55.00	73.68	0	MspA1I,PvuII	1	30	5968	NM_007376.4	57	4544	CDS	17
[以下略]
```

exonを判別している行とCDSを判別している行とが別々に出力されてしまう。これはわかりにくい。

これを回避するために、二段階に分けて`bedtools intersect`を実行してみる。
一段階目で ==**(A')**== と ==**(B2)**== の重なりを調べることによりexon番号を判別し、二段階目でさらに ==**(B1)**== との重なりを調べることによりCDSを判別する。最後に、不要な列を`cut -f`で削除。

```
% bedtools intersect -wao \
    -a <(bedtools intersect -wb \
           -a crisprdirect.bed \
           -b exon.bed) \
    -b cds.bed \
    | cut -f 1-12,16,20
```

結果：

```
NM_007376.4	4	26	+	GGATCAGAGTTCGGGGGCTGAGG	65.00	78.83	0		1	11	14844	exon_001	.
NM_007376.4	5	27	+	GATCAGAGTTCGGGGGCTGAGGG	60.00	76.96	0		1	9	9602	exon_001	.
NM_007376.4	40	62	+	TCTCTGCCCTCTCCACCATGAGG	60.00	77.23	0		1	122	76542	exon_001	CDS
NM_007376.4	46	68	-	CCCTCTCCACCATGAGGAGAAAC	50.00	71.73	0		1	108	11591	exon_001	CDS
NM_007376.4	47	69	-	CCTCTCCACCATGAGGAGAAACC	55.00	73.58	0		1	30	9285	exon_001	CDS
NM_007376.4	52	74	-	CCACCATGAGGAGAAACCAGCTG	55.00	73.68	0	MspA1I,PvuII	1	30	5968	exon_001	CDS
NM_007376.4	55	77	-	CCATGAGGAGAAACCAGCTGCCC	60.00	78.15	0	MspA1I,PvuII	1	49	7749	exon_001	CDS
NM_007376.4	68	90	-	CCAGCTGCCCACACCAGCTTTTC	60.00	77.65	0	MspA1I,PvuII	1	58	7374	exon_001	CDS
NM_007376.4	75	97	-	CCCACACCAGCTTTTCTTTTACT	35.00	66.17	0		1	39	5821	exon_001	CDS
NM_007376.4	76	98	-	CCACACCAGCTTTTCTTTTACTG	40.00	65.60	0		1	42	7370	exon_001	CDS
[以下略]
```

CRISPRdirectの結果の右に、exon番号の列とCDS判定の列が加わった！

### 多数の遺伝子に対してガイドRNAを設計する

準備中
