目的

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_ を除外。

パースの途中の状態

egrep を実行した直後の状態。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
//
[...]

最終的な出力。RefSeq ID, CDS の start と end, length が出力されている。

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 は並列化を支援する自作スクリプト。ファイルまたは標準入力からコマンドリストを読み込み並列実行する。コマンドリストは1行に1コマンド。

終了したら、

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

でマージ完了。

3′UTRの範囲の抽出 → BED化

上記のコマンドは、RefSeq ID, CDS の start と end, length を出力するので、3′UTR start, end を出力するよう改変。

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 はテキストを RefSeq ID の順番にソートする自作スクリプト。

結果

hs_refseq215_3utr.bed （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そのものになっているものがある。
- その場合は 3′UTR start = 3′UTR end になっている。

参考

GGGenomeでRefSeq転写産物の3′UTRを検索する