# 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を設計する
準備中
Sign in with Wallet
Connect another wallet