Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

gtfファイルからrefFlat形式への変換

みなさんgtfファイルからrefFlatに変換する時ってどうされてるんですかね?Rを使っている?
自分は自作のツール "gtf2refFlat" を使っているので、ここではそれを紹介します。

gtf形式については↓をご覧ください。

https://bi.biopapyrus.jp/rnaseq/mapping/gtf.html

gtfファイルは最も詳細な遺伝子アノテーションと思いますが、1行1遺伝子にはなっていないので使いづらいことがあります。
私が作ったツール群はだいたい遺伝子アノテーションファイルをrefFlat形式で受け付けているのですが、EnsemblはrefFlat形式を提供していないので、自作しました。

gtfファイルのダウンロード

Ensemblサイトから、主要な染色体の遺伝子のみ格納した Homo_sapiens.GRCh38.97.chr.gtf.gz をダウンロードします。

 # ダウンロード
 $ wget ftp://ftp.ensembl.org/pub/release-97/gtf/homo_sapiens/Homo_sapiens.GRCh38.97.chr.gtf.gz
 # 解凍
 $ gunzip Homo_sapiens.GRCh38.97.chr.gtf.gz
 # 確認
 $ head Homo_sapiens.GRCh38.97.chr.gtf
#!genome-build GRCh38.p12
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.27
#!genebuild-last-updated 2019-03
1       havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
1       havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; tag "basic"; transcript_support_level "1";
1       havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";

ChIPseqToolsのインストール

gtf2refFlatはChIPseqToolsの中に含まれています。 以下はgitがインストールされている前提です。gitについては以下の記事を参照ください。

GitHubからプログラムをダウンロード・インストール - Palmsonntagmorgen

 $ git clone --recursive https://github.com/rnakato/ChIPseqTools.git
 $ cd ChIPseqTools
 $ make

ChIPseqTools/bin の中に gtf2refFlat が生成されていればOKです。 ヘルプを表示してみましょう。

 $ ChIPseqTools/bin/gtf2refFlat

Options:
  -g [ --gtf ] arg      Gene file
  -n [ --name ]         Output name instead of id
  -u [ --unique ]       Only output one transcript per one gene (default: all
                        transcripts)
  -h [ --help ]         Print this message

transcript単位のrefFlatを出力

gtf2refFlatはデフォルトではisoformを個別に格納したrefFlatを出力します。

 $ ChIPseqTools/bin/gtf2refFlat -g Homo_sapiens.GRCh38.97.chr.gtf \
   > Homo_sapiens.GRCh38.97.transcript.refFlat
 # 確認
 $  head -n5 Homo_sapiens.GRCh38.97.transcript.refFlat
ENSG00000210196 ENST00000387461 chrMT   -       15956   16023   15956   16023   1       15956,  16023,  Mt_tRNA Mt_tRNA
ENSG00000198804 ENST00000361624 chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding protein_coding
ENSG00000209082 ENST00000386347 chrMT   +       3230    3304    3230    3304    1       3230,   3304,   Mt_tRNA Mt_tRNA
ENSG00000210077 ENST00000387342 chrMT   +       1602    1670    1602    1670    1       1602,   1670,   Mt_tRNA Mt_tRNA
ENSG00000210127 ENST00000387392 chrMT   -       5587    5655    5587    5655    1       5587,   5655,   Mt_tRNA Mt_tRNA

1列目はgene ID, 2列目がtranscript IDです。一番右の列には遺伝子のタイプが2列表示されていますが、右から2列目がgene type, 最右列がtranscript typeです。coding geneでもtranscriptがnoncodingであることはあり得ますので、このようになっています。

Ensembl IDでなく transcript name で出力したい場合は -n オプションをつけます。

$ ChIPseqTools/bin/gtf2refFlat -n -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat
$ head -n5 Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat
MT-TP   MT-TP-201       chrMT   -       15956   16023   15956   16023   1       15956,  16023,  Mt_tRNA Mt_tRNA
MT-CO1  MT-CO1-201      chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding  protein_coding
MT-TL1  MT-TL1-201      chrMT   +       3230    3304    3230    3304    1       3230,   3304,   Mt_tRNA Mt_tRNA
MT-TV   MT-TV-201       chrMT   +       1602    1670    1602    1670    1       1602,   1670,   Mt_tRNA Mt_tRNA
MT-TA   MT-TA-201       chrMT   -       5587    5655    5587    5655    1       5587,   5655,   Mt_tRNA Mt_tRNA

トップにミトコンドリアが来ているのがアレですが、とにかく出力できています。

Gene単位のrefFlatを出力

transcript単位ではなくgene単位でrefFlatを出力したい場合は -u オプションをつけます。遺伝子名にしたい場合は、transcriptの時と同様に -n オプションをつけます。

 $ ChIPseqTools/bin/gtf2refFlat -u -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.gene.chr.refFlat
 $ head -n5  Homo_sapiens.GRCh38.97.gene.chr.refFlat
ENSG00000270016 ENST00000602595 chr15   -       30607695        30608193        30607695        30608193        1      30607695,        30608193,       lncRNA  lncRNA
ENSG00000271983 ENST00000607458 chr15   -       80693216        80693707        80693216        80693707        1      80693216,        80693707,       lncRNA  lncRNA
ENSG00000222139 ENST00000410207 chr15   +       80663772        80663878        80663772        80663878        1      80663772,        80663878,       snRNA   snRNA
ENSG00000259175 ENST00000558208 chr15   -       80554609        80562944        80554609        80562944        3      80562862,80560056,80554609,      80562944,80560167,80556933,     lncRNA  lncRNA
ENSG00000184140 ENST00000328882 chr15   +       101803509       101806887       101805720       101806658       2      101803509,101805687,     101803545,101806887,    protein_coding  protein_coding
 # 遺伝子名を出力する場合
 $ ChIPseqTools/bin/gtf2refFlat -u -n -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.gene.chr.name.refFlat
 $ head -n5 Homo_sapiens.GRCh38.97.gene.chr.name.refFlat
AC026150.3      AC026150.3-201  chr15   -       30607695        30608193        30607695        30608193        1      30607695,        30608193,       lncRNA  lncRNA
AC023302.1      AC023302.1-201  chr15   -       80693216        80693707        80693216        80693707        1      80693216,        80693707,       lncRNA  lncRNA
RNU6-380P       RNU6-380P-201   chr15   +       80663772        80663878        80663772        80663878        1      80663772,        80663878,       snRNA   snRNA
AC108451.1      AC108451.1-201  chr15   -       80554609        80562944        80554609        80562944        3      80562862,80560056,80554609,      80562944,80560167,80556933,     lncRNA  lncRNA
OR4F6   OR4F6-201       chr15   +       101803509       101806887       101805720       101806658       2       101803509,101805687,    101803545,101806887,    protein_coding  protein_coding

遺伝子ごとに出力する場合、gtf2refFlatはgtfに含まれるtranscriptの中で最もcanonicalなもの (Consensus CDS (CCDS) タグを持つtranscriptの中で最も遺伝子長が長いもの)を出力します。左から2列目がそのtranscriptのID(または名前)になります。

おまけ:coding geneだけを出力したい場合

gawkを用いてcoding gene だけを出力してみます。一番右のアノテーションが"protein_coding"となっている行だけを出力するというワンライナーです。

$ cat Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat \
  | gawk '{ if ($13=="protein_coding") print }' \
  > Homo_sapiens.GRCh38.97.transcript.coding.name.chr.refFlat 
$ head -n5 Homo_sapiens.GRCh38.97.transcript.coding.name.chr.refFlat
MT-CO1  MT-CO1-201      chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding  protein_coding
MT-ND2  MT-ND2-201      chrMT   +       4470    5511    4470    5511    1       4470,   5511,   protein_coding  protein_coding
MT-ND6  MT-ND6-201      chrMT   -       14149   14673   14149   14673   1       14149,  14673,  protein_coding  protein_coding
MT-CO3  MT-CO3-201      chrMT   +       9207    9990    9207    9990    1       9207,   9990,   protein_coding  protein_coding
MT-ND4  MT-ND4-201      chrMT   +       10760   12137   10760   12137   1       10760,  12137,  protein_coding  protein_coding