# RefSeq転写産物の3'UTRの範囲をBEDに書き出す

### 目的

* GGGenomeでRefSeq転写産物の3'UTRを検索したい。
* GGGenomeは検索結果をBED形式で出力できるので、3'UTRの範囲を記載したBEDファイルを用意しておけば、`bedtools intersect` で3'UTRに対するヒットのみを得ることができる。

### RefSeq転写産物のgbffを入手

ftp://ftp.ncbi.nlm.nih.gov/refseq/release/complete/

から `complete.*.rna.gbff.gz` をダウンロード。

* `complete` は全生物種を含む
* `rna` が転写産物（ほかに `genomic`, `protein` がある）
* `gbff` はGenBank形式（配列＋アノテーション）

RefSeq 215 では、圧縮された状態で合計 71GB ある。

### gbffのパースの準備

#### コマンド

zcat complete.*.rna.gbff.gz \
  | egrep '^(LOCUS|VERSION|SOURCE|\s+CDS|//)' \

（メモ）gbffは大きいため、最初に `egrep` で必要な行だけを抜き出すことにより処理の時間を短縮。行末の `\` に続けて下記コマンドを順次実行。

  | perl -ne 'BEGIN{ $/ = "\n//\n" }
              my $length  = /^LOCUS\ .*?(\d+)\ *bp/m ? $1 : "" ;
              my $version = /^VERSION\ +(\S*)/m ? $1 : "" ;
              my $source  = /^SOURCE\ +(.*?)\ *$/m ? $1 : "" ;
              my ($cds1, $cds2) = /^\ {5}CDS\ +(\d+)\.\.(\d+)/m ? ($1, $2) : () ;
              print $source  . "\t" .
                    $version . "\t" .
                    $cds1    . "\t" .
                    $cds2    . "\t" .
                    $length  . "\n"' \

（メモ）gbffでは各エントリが `//` で区切られているので、perlの行セパレータを `//` に変更して、1エントリをまるごと「1行」として読み込んで高速化。1エントリごとに、配列の長さ、生物種、version、CDS の start と end を抽出していく。

  | grep '^Homo sapiens (human)' \
  | cut -f 2- \

（メモ）ヒト `Homo sapiens (human)` のみ抽出。

  | grep -v '^X'

（メモ）predictedな転写産物 `XM_`, `XR_` を除外。

#### パースの途中の状態

gbff から必要な行だけ抜き出したもの。

LOCUS       NM_001387374            4051 bp    mRNA    linear   PRI 21-APR-2022
VERSION     NM_001387374.1
SOURCE      Homo sapiens (human)
     CDS             167..2026
//
LOCUS       NM_001387375            4027 bp    mRNA    linear   PRI 21-APR-2022
VERSION     NM_001387375.1
SOURCE      Homo sapiens (human)
     CDS             167..2002
//
LOCUS       NM_001387376            3973 bp    mRNA    linear   PRI 21-APR-2022
VERSION     NM_001387376.1
SOURCE      Homo sapiens (human)
     CDS             167..1948
//
[...]

最終的な出力。

NM_001387562.1  221     2017    4790
NM_001387563.1  221     1219    1358
NR_073515.3                     1315
NR_170667.1                     6439
NR_047672.2                     488
NR_038873.1                     2636
NM_001141920.2  236     781     2901
NM_001363651.3  171     1949    2081
NM_001365041.3  510     1526    3569
NM_175569.3     236     778     2898
NM_001367912.2  39      1844    1986
NM_001387759.1  243     2255    2923
NM_001387772.1  160     2079    2747
[...]

BED形式のように見えるが開始位置が 0-based start でない点には注意。

`NM_` で始まるものは mRNA なので CDS の記載あり。
`NR_` で始まるものは ncRNA なので CDS の記載なし。

#### （2023/08/22追記）

OAZ1, OAZ2, OAZ3, PEG10 遺伝子は ribosomal frameshift が起こるため、

     CDS             join(480..1436,1436..2605)

のように記載されている。そのような場合にも対応できるよう下記のように修正。

  | perl -ne 'BEGIN{ $/ = "\n//\n" }
              my $length  = /^LOCUS\ .*?(\d+)\ *bp/m ? $1 : "" ;
              my $version = /^VERSION\ +(\S*)/m ? $1 : "" ;
              my $source  = /^SOURCE\ +(.*?)\ *$/m ? $1 : "" ;
              my ($cds1, $cds2) = /^\ {5}CDS\ +(\d+)\.\.(\d+)/m ? ($1, $2) :
                                  /^\ {5}CDS\ +join\((\d+)\.\.\d+,\d+\.\.(\d+)\)/m ? ($1, $2) :
                                  () ;
              print $source  . "\t" .
                    $version . "\t" .
                    $cds1    . "\t" .
                    $cds2    . "\t" .
                    $length  . "\n"' \

### （並列化）

そのまま実行すると時間がかかるため、並列化する方法もメモ。

下記を `/tmp/parse.sh` に保存して `chmod 755` しておく。

egrep '^(LOCUS|VERSION|SOURCE|\s+CDS|//)' \
| perl -ne 'BEGIN{ $/ = "\n//\n" }
            my $length  = /^LOCUS\ .*?(\d+)\ *bp/m ? $1 : "" ;
            my $version = /^VERSION\ +(\S*)/m ? $1 : "" ;
            my $source  = /^SOURCE\ +(.*?)\ *$/m ? $1 : "" ;
            my ($cds1, $cds2) = /^\ {5}CDS\ +(\d+)\.\.(\d+)/m ? ($1, $2) :
                                /^\ {5}CDS\ +join\((\d+)\.\.\d+,\d+\.\.(\d+)\)/m ? ($1, $2) :
                                () ;
            print $source  . "\t" .
                  $version . "\t" .
                  $cds1    . "\t" .
                  $cds2    . "\t" .
                  $length  . "\n"' \
| grep '^Homo sapiens (human)' \
| cut -f 2- \
| grep -v '^X'

下記を実行。

/bin/ls *gbff* \
  | awk '{print "zcat "$1" | /tmp/parse.sh > /tmp/refseq215/"$1".txt" }' \
  | /backup/bin/parallel.pl -12

（メモ）[parallel.pl](https://github.com/meso-cacase/parallel/blob/master/parallel.pl) は並列化を支援する自作スクリプト。

終了したら、

setopt numeric_glob_sort
cat /tmp/refseq215/* > /tmp/hs_refseq215_cds.txt

でマージ完了。

### 3'UTRの範囲の抽出 → BED化

* CDS end の値をそのまま 3'UTR start とする。
* 0-based start のBEDにするため。
* length の値を 3'UTR end とする。
* 出力の対象を mRNA `MN_` のみとする。

`/tmp/parse.sh` を下記のように修正して同様に実行。

egrep '^(LOCUS|VERSION|SOURCE|\s+CDS|//)' \
| perl -ne 'BEGIN{ $/ = "\n//\n" }
            my $length  = /^LOCUS\ .*?(\d+)\ *bp/m ? $1 : "" ;
            my $version = /^VERSION\ +(\S*)/m ? $1 : "" ;
            my $source  = /^SOURCE\ +(.*?)\ *$/m ? $1 : "" ;
            my ($cds1, $cds2) = /^\ {5}CDS\ +(\d+)\.\.(\d+)/m ? ($1, $2) :
                                /^\ {5}CDS\ +join\((\d+)\.\.\d+,\d+\.\.(\d+)\)/m ? ($1, $2) :
                                () ;
            print $source  . "\t" .
                  $version . "\t" .
                  $cds2    . "\t" .
                  $length  . "\n"' \
| grep '^Homo sapiens (human)' \
| cut -f 2- \
| grep '^NM_'

/bin/ls *gbff* \
  | awk '{print "zcat "$1" | /tmp/parse.sh > /tmp/refseq215/"$1".txt" }' \
  | /backup/bin/parallel.pl -12

並列処理が終了したら、

cat /tmp/refseq215/* \
  | /backup/bin/sort_refseq_id.pl \
> /tmp/hs_refseq215_3utr.bed

でマージ完了。

（メモ）[sort_refseq_id.pl](https://gist.github.com/meso-cacase/7a287171a6b364b5b487bdb55c66b251) はテキストを RefSeq ID の順番にソートする自作スクリプト。

### 結果

[hs_refseq215_3utr.bed](https://gist.github.com/meso-cacase/97301382511db5e5dd61366991dda8d7) （66,326行）

NM_000014.6     4495    4610
NM_000015.3     943     1285
NM_000016.6     1345    2261
NM_000017.4     1299    1859
NM_000018.4     2015    2184
[...]

* BED形式なので 0-based start である点に注意。
* CDSしか登録されていないエントリがある。
* その場合は start = end になっている。

### 完了