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](https://gggenome.dbcls.jp/) では、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://qiita.com/meso_cacase/items/f83a4f2786a5188cd92b) を参照
  * 完成版： [hs_refseq215_3utr.bed](https://gist.github.com/meso-cacase/97301382511db5e5dd61366991dda8d7) （66,326行）

### GGGenomeによる検索例

たとえば、下記のようなsiRNAのseed `UGAACUC` と相補的な `GAGUUCA` という配列を3′UTRにもつ（seedマッチする）mRNAを検索する。

![siVIM-270](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69129/1fb9d837-3922-f20c-0179-0c05ee6adfef.png)

* `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 にはseedマッチが5箇所あるか？
* オレンジ部分がCDS
* 緑色のハイライトがseedマッチ → OK

![NM_001348186.2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69129/77c4bc04-9d19-4e2c-8539-7596bc51dad3.png)

#### 参考

* [RefSeq転写産物の3′UTRの範囲をBEDに書き出す](https://qiita.com/meso_cacase/items/f83a4f2786a5188cd92b)
* [マイクロアレイ学生実習2014 - mesoの実験ノート](https://data.dbcls.jp/~meso/meme/archives/2788/)