StringTie

選択的スプライシングにより 1 つの遺伝子から複数種類の転写産物が生成される場合がある。そのため、マッピング結果から遺伝子発現量を見積もるときに、実験目的に応じて転写産物の isoform を無視して遺伝子ごとに計上する場合と、isoform を考慮して isoform 毎に計上する場合がある。前者の場合は、featureCounts が一般的に使われている。これに対して、後者の場合は StringTie や RSEM を利用するのが一般的である。

StringTie を使用した発現量定量の流れとして、(1) BAM ファイルから新規 isoform と思われるもののアノテーションを作成し、(2) 既存のアノテーションと新しいアノテーションをマージしてから発現量の定量を行う、の 2 ステップからなる。ここで計算される遺伝子発現量は TPM および FPKM として出力される。この発現量に対して引き続き発現変動遺伝子を検出したい場合は Ballgown を使用する。

StringTie

StringTie を利用することで、HISAT2 などのマッピング結果結果を利用して、トランスクリプトのアセンブリを行い、転写産物 isoform 毎の発現量を推定できるようになる。このとき、既知のアノテーションをアセンブリのガイドとして与えることができる。アセンブリは、StringTie に BAM ファイルを与えることで、そのアセンブリの結果が GTF ファイルで出力される。この処理を個々のサンプルに対して行う必要がある。例えば、SRR1792924, SRR1792925, SRR1792926, SRR1792927, SRR1792928, SRR1792929 の 6 サンプルであれば、次のように 6 回 StringTie を実行する。

for libid in 24 25 26 27 28 29
do
  stringtie -p 4 -G Arabidopsis_thaliana.TAIR10.40.gff3 \
          -o SRR17929${libid}_tari10.gtf \
          -l SRR17929${libid} \
          SRR17929${libid}.bam
done

StringTie により、サンプル毎のアノテーションが作られる。次に、StringTie の --merge モードを使用して、これらのアノテーションをマージしていく。まず、上で作った 6 つのアノテーションファイル(.gtf)へのパスをテキストファイルに記載し、次にこのテキストファイルを StringTie に与えてマージを行う。


echo -e "" > merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
echo "SRR1792924_tari10.gtf" >> merge_list.txt
​
stringtie --merge -p 16 -G Arabidopsis_thaliana.TAIR10.40.gff3 \
                  -o tair10_stringtie_merged.gtf \
                  merge_list.txt

プログラムが正しく実行されると、tair10_stringtie_merged.gtf ファイルが作られる。このファイルの中に、今回のデータから観測できた新しい isoform と思われるアノテーションも含まれている。実際に、アセンブリ前とアセンブリ後のアノテーション数の変化を調べるには、gffcompare を使用する。

gffcompare -r Arabidopsis_thaliana.TAIR10.40.gff3 -G \
           -o stringtie_merge tair10_stringtie_merged.gtf
# # gffcompare v0.11.2 | Command line was:
#
# = Summary for dataset: tair10_stringtie_merged.gtf
#     Query mRNAs :   58352 in   32246 loci  (46241 multi-exon transcripts)
#            (11591 multi-transcript loci, ~1.8 transcripts per locus)
#  Reference mRNAs :   54235 in   32619 loci  (42815 multi-exon)
#  Super-loci w/ reference transcripts:    32070
#  -----------------| Sensitivity | Precision  |
#         Base level:   100.0     |    98.8    |
#         Exon level:    97.8     |    96.5    |
#       Intron level:   100.0     |    97.8    |
# Intron chain level:   100.0     |    92.6    |
#   Transcript level:    99.5     |    92.5    |
#        Locus level:    99.7     |    99.3    |
#      Matching intron chains:   42815
#        Matching transcripts:   53983
#               Matching loci:   32536
#           Missed exons:       0/193121  (  0.0%)
#            Novel exons:     900/192441  (  0.5%)
#         Missed introns:       1/132525  (  0.0%)
#        Novel introns:     931/135466  (  0.7%)
#           Missed loci:       0/32619  (  0.0%)
#            Novel loci:     176/32246  (  0.5%)
#
# Total union super-loci across all input datasets: 32246
# 58352 out of 58352 consensus transcripts written in stringtie_merge.annotated.gtf (0 discarded as redundant)

アノテーションのマージが完了してから、次にこの新しいアノテーションを利用して、BAM ファイルを処理し、発現量の推定を行う。発現量の推定は StringTie に -e オプションを付けて実行する。

for libid in 24 25 26 27 28 29
do
  stringtie -e -B -p 4 -G tair10_stringtie_merged.gtf \
            -o ballgown/SRR17929${libid}/SRR17929${libid}_stringtie.gtf \
            SRR17929${libid}.bam
done

プログラムが正しく実行されると、各サンプルに対して 1 つのフォルダが作られる。それぞれのフォルダの中に 7 個のファイルが含まれている。これらのファイルは、isoform id から transcript id に変換するテーブル、exon id から transcript id に変換するテーブルなどが含まれている。また、この中にある gtf ファイルには、発現量として FPKM と TPM が記載されている。

cd  ballgown/SRR1792924
# total 84M
# -rw-r--r-- 1 biopapyrus biopapyrus 4.0M Aug  4 07:24 e2t.ctab
# -rw-r--r-- 1 biopapyrus biopapyrus  14M Aug  4 07:24 e_data.ctab
# -rw-r--r-- 1 biopapyrus biopapyrus 3.3M Aug  4 07:24 i2t.ctab
# -rw-r--r-- 1 biopapyrus biopapyrus 5.1M Aug  4 07:24 i_data.ctab
# -rw-r--r-- 1 biopapyrus biopapyrus  54M Aug  4 07:24 SRR1792924_stringtie.gtf
# -rw-r--r-- 1 biopapyrus biopapyrus 5.0M Aug  4 07:24 t_data.ctab
# drwx------ 2 biopapyrus biopapyrus 4.0K Aug  4 07:21 tmp.XXehOX9J
​
head SRR1792924_stringtie.gtf
# 1 StringTie transcript  3631  5899  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; cov "19.784952"; FPKM "4.125814"; TPM "8.077307";
# 1 StringTie exon  3631  3913  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "1"; cov "19.091873";
# 1 StringTie exon  3996  4276  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "2"; cov "24.274023";
# 1 StringTie exon  4486  4605  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "3"; cov "30.108334";
# 1 StringTie exon  4706  5095  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "4"; cov "20.966667";
# 1 StringTie exon  5174  5326  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "5"; cov "20.751635";
# 1 StringTie exon  5439  5899  1000  + . gene_id "MSTRG.1"; transcript_id "transcript:AT1G01010.1"; exon_number "6"; cov "13.466377";
# 1 StringTie transcript  6788  9130  1000  - . gene_id "MSTRG.2"; transcript_id "transcript:AT1G01020.1"; cov "13.410775"; FPKM "2.796588"; TPM "5.475017";

発現量の推定がおわると、必要に応じて発現量を取り出してモデリングを行なったり、あるいは発現変動遺伝子を検出したりすることができる。StringTie から算出された FPKM を利用した発現変動遺伝子の検出は、Ballgown の使用が推奨されている。