# GGGenomeでRefSeq転写産物の3’UTRを検索する
### 目的
* siRNAやmiRNAがseed matchするような転写産物を検索したい。
* **seed配列に依存的なオフターゲット作用** の評価に。
* mRNAの全長ではなく **3’UTR** の配列だけを検索したい。
※ seedとは、siRNAのガイド鎖（アンチセンス鎖）またはmatureなmiRNAの、5'末端から数えて2〜8番目の7merを中心とする配列。文献により2〜7番目の6mer、2〜9番目の8merなどとされる場合もある。seedと相補的な配列がmRNAの3'UTRに存在すると、そのmRNAの分解や翻訳抑制が起こる場合がある。
### 方針を検討
* GGGenomeでは、RefSeqに収録された全生物種の転写産物 `NM_`, `NR_`, `XM_`, `XR_` を検索することができる。ただし最新のRefSeqリリースのみ。
* https://gggenome.dbcls.jp/refseq/
* RefSeqリリースのキリ番（200, 205, 210, 215, ...）については、ヒトの `NM_`, `NR_` に限定した検索が用意されている。さらに、RefSeqが更新されてもキリ番リリースは削除されずに残されている。
* 例： https://gggenome.dbcls.jp/hsnm_refseq215/
* 今回は、ヒトの `NM_`, `NR_` をGGGenomeで検索し、そこから3'UTRにヒットしたものを抽出することにする。
* GGGenomeの検索結果をBED形式で取得し、3'UTRの範囲を記載したBEDファイルと `bedtools intersect` すればよい。
### 3’UTRの範囲を記載したBEDファイルを準備
* [RefSeq転写産物の3’UTRの範囲をBEDに書き出す](https://hackmd.io/@meso/ByQ3Q_gan) を参照
* 完成版： [hs_refseq215_3utr.bed](https://gist.github.com/meso-cacase/97301382511db5e5dd61366991dda8d7) （66,326行）
### GGGenomeによる検索例
たとえば、下記のsiRNAのseed `UGAACUC` と相補的な `GAGUUCA` という配列を3′UTRにもつ（seedマッチする）mRNAを検索する。
* `U` を `T` に置き換え、`GAGTTCA` をGGGenomeで検索する。
* https://gggenome.dbcls.jp/hsnm_refseq215/+/GAGTTCA
* 下記のように検索しても同じ結果が得られる。
* https://gggenome.dbcls.jp/hsnm_refseq215/-/UGAACUC
* `U` は `T` に変換される。`/-/` で相補的な配列を検索。
* URLの最後に `.bed` をつければ検索結果をBED形式で取得できる。
* https://gggenome.dbcls.jp/hsnm_refseq215/-/UGAACUC.bed
* 不要な部分を削除して使いやすくしておく。
* `sed 1d` で1行目（ヘッダ）を削除
* `awk` で RefSeq ID, start, end のみ抜き出して出力
```
curl -s https://gggenome.dbcls.jp/hsnm_refseq215/-/UGAACUC.bed \
| sed 1d \
| awk -F "\t" '{ OFS="\t" ; gsub(/ .*/, "", $1) ; print $1,$2,$3 }' \
> GGGenome.bed
```
（結果 = 28,288行）
```
NR_028389.1 1998 2005
NR_003928.2 2210 2217
NR_028028.2 3776 3783
NR_110015.1 846 853
NR_110014.1 914 921
NR_002714.1 533 540
NM_001256867.1 1344 1351
NM_001242328.1 1344 1351
NM_001256873.1 1344 1351
NM_001242326.1 1344 1351
[...]
```
このファイルから3’UTRにヒットしたものを抽出する。
```
bedtools intersect -wa -f 1 -a GGGenome.bed -b hs_refseq215_3utr.bed \
> GGGenome.UTR3.bed
```
（メモ）
* `-f 1` なし：1塩基でも重なれば出力する。
* `-f 1` あり：全長が3'UTRに含まれる場合のみ出力する。
* 今回はこちらを採用。3'UTRの配列DBを作成して検索した場合と同じ結果が得られるようにする。
（結果 = 11,059行）
```
NM_032333.5 3691 3698
NM_014506.3 1682 1689
NM_173555.4 3111 3118
NM_019596.6 1496 1503
NM_001040273.3 2794 2801
NM_001162496.3 1415 1422
NM_020398.4 863 870
NM_001305225.4 7177 7184
NM_015090.4 9372 9379
NM_001160332.2 9387 9394
[...]
```
#### 同じ転写産物に対するヒット数
同じ転写産物の3'UTRに複数のseedマッチがある場合は強く抑制される傾向があるので、転写産物ごとのヒット数を調べられるようにしておく。
* `awk` でカウントしてみる。
```
cat GGGenome.UTR3.bed \
| cut -f 1 \
| awk -v 'OFS=\t' 'a[$0]++ {} END {
for( key in a ){ print key, a[key] }
}'
```
（結果 = 9,039行）
```
NM_001375706.1 1
NM_001010989.3 1
NM_001350659.2 3
NM_016089.3 1
NM_006386.5 1
NM_017970.4 2
NM_198395.2 1
NM_003012.5 1
NM_000643.2 1
NM_001348186.2 5
[...]
```
（確認）NM_001348186.2 の 3'UTR には5箇所seedマッチしているか？ → OK
* オレンジ部分がCDS
Sign in with Wallet
Connect another wallet