選択的スプライシングにより 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
doneStringTie により、サンプル毎のアノテーションが作られる。次に、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 の使用が推奨されている。