目的
- 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 になっている。