遺伝子アノテーションファイルの処理

GTFファイル

GTF (gene feature format) ファイルは遺伝子アノテーションを記述しているファイルである。ファイルの 1 行は 1 つの feature を表す。ここでいう feature は、exon、intron、tRNA などのゲノム上の領域を指す。

各行には 9 項目のデータが記述されている。各項目は「タブ」によって区切られている。各項目は以下の順で記述される。9 項目をすべて記述することになっている。不明な場合はピリオド(.)を記入する。

  1. seqname - 染色体名
  2. source- データの種類またはプログラム名など
  3. feature - intron、exon、CDS、RNA、miRNA など
  4. start - feature の開始位置
  5. end - feature の終了位置
  6. score - ?
  7. strand - 鎖の方向( + または - )
  8. frame - コドンの開始位置( 0、1 または 2 )
  9. attribute - 属性など。Ensembl ID、遺伝子汎用名、エキソンの番号などが記載されている

実際の GTP ファイルをテキストエディタで開くと、以下のようになっている。

GL000213.1	protein_coding	CDS	138767	139287	.	-	0	 gene_id "ENSG00000237375"; transcript_id "ENST00000327822"; exon_number "1"; gene_name "BX072566.1"; gene_biotype "protein_coding"; transcript_name "BX072566.1-201"; protein_id "ENSP00000329990";
GL000213.1	protein_coding	start_codon	139285	139287	.	-	0	 gene_id "ENSG00000237375"; transcript_id "ENST00000327822"; exon_number "1"; gene_name "BX072566.1"; gene_biotype "protein_coding"; transcript_name "BX072566.1-201";
GL000213.1	protein_coding	exon	134276	134390	.	-	.	 gene_id "ENSG00000237375"; transcript_id "ENST00000327822"; exon_number "2"; gene_name "BX072566.1"; gene_biotype "protein_coding"; transcript_name "BX072566.1-201";
GL000213.1	protein_coding	CDS	134276	134390	.	-	1	 gene_id "ENSG00000237375"; transcript_id "ENST00000327822"; exon_number "2"; gene_name "BX072566.1"; gene_biotype "protein_coding"; transcript_name "BX072566.1-201"; protein_id "ENSP00000329990";
GL000213.1	protein_coding	exon	133943	134116	.	-	.	 gene_id "ENSG00000237375"; transcript_id "ENST00000327822"; exon_number "3"; gene_name "BX072566.1"; gene_biotype "protein_coding"; transcript_name "BX072566.1-201";

例えば 1 行目を見ると、これ ID が GL000213.1 染色体上のアノテーションであることがわかる。そして、この染色体上の 138767 番目の塩基から 139287 番目の塩基までがタンパク質をコーディングする領域(protein_coding)で、実際に翻訳される部分(CDS)である。また、この部分は遺伝子 ENSG00000237375 の一部である。

awk を利用した GTF ファイルの処理

GTF ファイルの記述規則はタブ区切りで非常に明快である。Perl や Python などを利用することで簡単に処理できる。ここでは、awk とよばれるスクリプトを利用して、GTP ファイルから必要な情報を取り出す方法を示す。利用する GTF ファイルは Ensembl FTP サイトからダウンロードできる。ダウンロードしたファイルを展開して利用する。

Ensembl ID と遺伝子名の対応リストの作成

awk を利用して行全体に対してマッチングを行い、後方参照により遺伝子 ID と名前を出力する。

gawk '{print gensub(/^.+gene_id "([^"]+)".+gene_name "([^"]+)".+$/,"\\1\t\\2", "")}' Homo_sapiens.GRCh37.68.gtf

awk スクリプトの実行結果。

ENSG00000265283	MIR3118-5
ENSG00000237375	BX072566.1
ENSG00000237375	BX072566.1
ENSG00000237375	BX072566.1

重複する部分があるため、awk の連想配列を利用して重複を削除する。

gawk '
   {print gensub(/^.+gene_id "([^"]+)".+gene_name "([^"]+)".+$/,"\\1\t\\2", "")}
' Homo_sapiens.GRCh37.68.gtf |
awk '
  BEGIN{FS="\t"}
  {if (a[$0]++ == 0) print $0}
'

awk の実行結果。

ENSG00000265283	MIR3118-5
ENSG00000237375	BX072566.1
ENSG00000238432	U6
ENSG00000262510	C21orf74

エキソンの情報を抽出

GTF ファイルから遺伝子名、エキソンの番号と位置を抽出する例。

gawk '
   BEGIN{FS="\t";OFS="\t"}
   {print gensub(/^.+exon_number "([^"]+)".+gene_name "([^"]+)".+$/,"\\2\t\\1", ""), $4, $5}
' Homo_sapiens.GRCh37.68.gtf

awk の実行結果。遺伝子名、エキソンの番号、開始位置、終了位置の順で出力される。

MIR3118-5	1	104742	104817
BX072566.1	1	138767	139339
BX072566.1	1	138767	139287
BX072566.1	1	139285	139287
BX072566.1	2	134276	134390
BX072566.1	2	134276	134390

RNA のアノンテーション情報を抽出

ヒトゲノムの GTF ファイルから RNA に関するアノンテーション情報だけを抽出する例。

awk 'BEGIN{FS="\t"}$2 ~ /RNA/{print $0}' Homo_sapiens.GRCh37.68.gtf

さらに対象を絞って miRNA のアノンテーション情報だけを取り出す場合は次のようにする。

awk 'BEGIN{FS="\t"}$2 ~ /^miRNA$/{print $0}' Homo_sapiens.GRCh37.68.gtf

結果。

GL000213.1	miRNA	exon	104742	104817	.	+	.	 gene_id "ENSG00000265283"; transcript_id "ENST00000578976"; exon_number "1"; gene_name "MIR3118-5"; gene_biotype "miRNA"; transcript_name "MIR3118-5-201";
HSCHR21_4_CTG1_1	miRNA	exon	34836774	34836866	.	+	.	 gene_id "ENSG00000266367"; transcript_id "ENST00000585000"; exon_number "1"; gene_name "AC207981.1"; gene_biotype "miRNA"; transcript_name "AC207981.1-201";
HSCHR1_2_CTG31	miRNA	exon	155221944	155222032	.	-	.	 gene_id "ENSG00000262403"; transcript_id "ENST00000572520"; exon_number "1"; gene_name "AL713999.2"; gene_biotype "miRNA"; transcript_name "AL713999.2-201";
HSCHR6_MHC_QBL	miRNA	exon	28743697	28743804	.	+	.	 gene_id "ENSG00000265877"; transcript_id "ENST00000580289"; exon_number "1"; gene_name "AL662890.13"; gene_biotype "miRNA"; transcript_name "AL662890.13-201";