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
