gtf から転写物配列取得 "gffread"

この記事では、gtf を使って、リファレンス配列から転写物の配列を切り出す方法を紹介します。

gffread は cufflinks の昨日のうちの一つで、
gtf <-> gff の変換もできます。

以下の記事を参考に実際に gffread を使って転写物の配列を切り出します。
gffread を使った transcripts fasta 転写物の配列取得 – バイオインフォ道場 [bioinfo-Dojo]

例として、CLOCK 遺伝子のトランスクリプトの配列の切り出しを行います。

取得したい転写物の位置情報 gtf ファイルの準備

まず、欲しい転写物の位置情報のみを gtf ファイルから切り出します。

gtf ファイルから CLOCK の遺伝子名で検索します。

$ gzcat FANTOM_CAT.lv4_stringent.gtf.gz | grep CLOCK  

chr4    FANTOM  gene    56292003    56413670    .   -   .   gene_id "ENSG00000134852.10"; geneSuperClass "all_mRNA";  geneClass "coding_mRNA";  geneSubClass "protein_coding"; gene_type "protein_coding"; gene_name "CLOCK"; coding_status "coding"; cumulative_support "ENCODE:FANTOM:GENCODE:HUBDMAP:MITRANS"; geneCategory "coding_mRNA"; DHS_type "DHS_promoter";
chr4    FANTOM  transcript  56300990    56412783    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000381322.1"; transcript_type "protein_coding"; transcript_name "CLOCK-002"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "65.83909";
chr4    FANTOM  transcript  56348940    56412783    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000509151.1"; transcript_type "processed_transcript"; transcript_name "CLOCK-009"; coding_status "nonCoding"; cumulative_support "GENCODE:FANTOM"; TIEScore "75.73369";
chr4    FANTOM  transcript  56345014    56413048    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000506923.1"; transcript_type "processed_transcript"; transcript_name "CLOCK-011"; coding_status "nonCoding"; cumulative_support "GENCODE:FANTOM"; TIEScore "80.83459";
chr4    FANTOM  transcript  56294070    56376198    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000309964.4"; transcript_type "protein_coding"; transcript_name "CLOCK-005"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "67.56291";
chr4    FANTOM  transcript  56352561    56412783    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513033.1"; transcript_type "processed_transcript"; transcript_name "CLOCK-008"; coding_status "nonCoding"; cumulative_support "GENCODE:FANTOM"; TIEScore "76.37197";
chr4    FANTOM  transcript  56319927    56376198    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000506747.1"; transcript_type "retained_intron"; transcript_name "CLOCK-007"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "62.99050";
chr4    FANTOM  transcript  56329865    56412783    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000435527.2"; transcript_type "protein_coding"; transcript_name "CLOCK-010"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "75.56549";
chr4    FANTOM  transcript  56301189    56413300    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; transcript_type "protein_coding"; transcript_name "CLOCK-001"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "74.60385";

gene と transcript の位置情報が得られましたが、これでは exon の情報がないので、
transcript_id で検索し直します。
CLOCK-001 の transcript_id "ENST00000513440.1" で grep して clock.gtf とでもして保存しましょう。

$ gzcat FANTOM_CAT.lv4_stringent.gtf.gz | grep ENST00000513440.1
chr4    FANTOM  transcript  56301189    56413300    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; transcript_type "protein_coding"; transcript_name "CLOCK-001"; coding_status "coding"; cumulative_support "GENCODE:FANTOM"; TIEScore "74.60385";
chr4    FANTOM  exon    56301189    56301761    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 1;
chr4    FANTOM  exon    56304449    56304704    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 2;
chr4    FANTOM  exon    56308599    56308801    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 3;
chr4    FANTOM  exon    56309854    56310063    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 4;
chr4    FANTOM  exon    56310800    56310952    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 5;
chr4    FANTOM  exon    56314946    56315035    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 6;
chr4    FANTOM  exon    56315563    56315663    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 7;
chr4    FANTOM  exon    56316258    56316399    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 8;
chr4    FANTOM  exon    56319221    56319296    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 9;
chr4    FANTOM  exon    56319844    56319991    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 10;
chr4    FANTOM  exon    56322064    56322170    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 11;
chr4    FANTOM  exon    56322385    56322467    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 12;
chr4    FANTOM  exon    56325059    56325177    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 13;
chr4    FANTOM  exon    56325315    56325428    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 14;
chr4    FANTOM  exon    56329852    56329972    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 15;
chr4    FANTOM  exon    56336884    56336973    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 16;
chr4    FANTOM  exon    56342130    56342221    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 17;
chr4    FANTOM  exon    56344982    56345130    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 18;
chr4    FANTOM  exon    56345807    56345866    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 19;
chr4    FANTOM  exon    56348906    56348995    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 20;
chr4    FANTOM  exon    56355541    56355632    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 21;
chr4    FANTOM  exon    56376079    56376232    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 22;
chr4    FANTOM  exon    56412949    56413300    .   -   .   gene_id "ENSG00000134852.10"; transcript_id "ENST00000513440.1"; exon_number 23;

複数の配列が欲しい時は, テキストファイルに transcript id を書き込んだものを保存しておいて、 grep -f transcript_id.txt annot.gtf
とすれば複数の転写物の位置情報を記述した gtf が作れます。
今回は例として１つの転写物の配列を取得しますが、こうすれば、以下の作業で multi-fasta を取得できます。

gtf から転写物の配列取得

以上で得られた CLOCK ( ENST00000513440.1 )転写物の位置情報だけを切り出した、 "clock.gtf" を用いて、リファレンスから転写物の配列を取得していきます。
以下のようにして、転写物の配列を取得します。

$ gffread -w clock.fa -g genome.fa clock.gtf

-w で出力先のファイル名を指定して、-g でリファレンス配列を指定します。

clock.fa に配列がしっかり保存できているかどうか確認します。

$ head -3  clock.fa
>ENST00000513440.1 gene=ENSG00000134852.10
AGTAACCGGCGCCGTTCCCGGCCGGGGCAGGGAACGGTGCGCCTGCAGAGCCAGATTTCGGCCCAAGGGG
CGCGGGAGTCTCTTCGGGCGTCCGGGATCCCCTGGCGCGGCTCCGTGCTGCCTAACGGGGCAAGTCGCAT

ちゃんと配列が保存されていそうですね。

しかし、切り出されてきた配列が本当に CLOCK 遺伝子の転写物の配列と一致するのかはこれで見ただけでは
不安なので一応 BLAST で確認して見ましょう。