超初心者向け！！RNA-seq解析シリーズ⑤StringTieの使い方

Q: StringTieはbrewとconda、どちらで入れればいい？

どちらでも構いません。Macで環境を分ける必要がなければ brew install stringtie が手軽です。解析環境をプロジェクトごとに分けたい場合や、Linux／Google Colabで作業する場合は conda install -c bioconda stringtie を使います。同じマシンに両方入れると which stringtie でどちらが呼ばれるか分からなくなるので、どちらかに統一するのがおすすめです。

Q: stringtie --merge が動かない／エラーになる

まず --merge のハイフンが2本になっているか確認してください。Webページからコピーすると、ハイフン2つが「—」（ダッシュ記号）に変換されてしまうことがあります。次に mergelist.txt のパスです。作業ディレクトリ以外にgtfファイルがある場合は、mergelist.txt に書いたパスがカレントディレクトリから見て正しいかを確認します。テキストエディットで作った場合は、リッチテキスト（.rtf）になっていないかも確認してください。

Q: -e と -B は何をしているの？

-B はBallgownが読み込む形式のテーブル（.ctab ファイル群）を出力させるオプションです。-e は -G で指定したGTFに含まれる転写産物だけを対象にし、新規の転写産物を推定しないようにします。Ballgownに渡す最終ステップでは、この2つをセットで使うのが定石 です。統合GTFの転写産物セットに固定して各サンプルを測るという意味になります。

Q: samtools sortでメモリ不足になる

-m オプションでスレッドあたりの使用メモリを指定できます（例：-m 1G）。既定は768MBで、-@ で指定したスレッド数の分だけ掛け算で消費される点に注意してください。-@ 4 -m 2G なら合計8GB使います。搭載メモリが少ない場合はスレッド数を減らすか -m を下げます。

Q: StringTieの出力（FPKM/TPM）をそのまま発現量として使っていい？

サンプル内で遺伝子間を比べる分には使えますが、サンプル間の発現差を統計的に検定するには専用のツールが必要 です。このシリーズでは次の記事でBallgownを使います。カウントベースで解析したい場合はDESeq2やedgeRが標準的な選択肢になります。 次の記事ではballgownを用いて発現差解析を行っていきたいと思います。 > [!POINT] > ここまでの samtools → StringTie → merge → Ballgown用ファイル作成を、サンプルごとにコマンドを並べる代わりに一括実行する 方法もあります。nf-core/rnaseq なら定量まで通しで走り、途中で失敗しても -resume で続きから再開できます。 Nextflowとnf-core/rnaseqでbulk RNA-seq解析を自動化する：QC・STAR・Salmonを再現可能に実行

誰でもできるRNA-seq解析シリーズ！
今回はHISAT2によるマッピング結果をアセンブルし、発現差を比較するようにします。

StringTieというソフトを使用します。

sam→bamへの変換

まず、アノテーションを実施する前に、
HISAT2でマッピングしたファイルをソートしてバイナリファイルのbam形式に変換します。

こうすることによってコンピュータが認識しやすくなり、
高速に処理することができるそうです。

ソート&bam変換は**「samtools」**を使います。

samtoolsのインストール

まずは、samtoolsをインストールします。
インストールはbrew installで簡単にできます。

brew install samtools

condaを使っている場合は、bioconda から入れられます。

conda install -c bioconda samtools

samファイルをソートしてbamに変換

samtools sort -@ 4 -O bam -o SRR1571967.sort.bam SRR1571967.sam

samtoolsのコマンド説明

samtools samtoolコマンド実行を指示
sort samtoolでソートする
-@ 4 CPUの使用するスレッド数を指定
-O bam outputのファイル形式をbamに指定
-o SRR1571967.sort.bam outputのファイル名を指定
SRR1571967.sam hisat2で出力したマッピングファイルを指定

すべてのファイルに対して同様に実施しましょう。

samtools sort -@ 4 -O bam -o SRR1571967.sort.bam SRR1571967.sam
samtools sort -@ 4 -O bam -o SRR1571968.sort.bam SRR1571968.sam
samtools sort -@ 4 -O bam -o SRR1571969.sort.bam SRR1571969.sam
samtools sort -@ 4 -O bam -o SRR1571970.sort.bam SRR1571970.sam
samtools sort -@ 4 -O bam -o SRR1571971.sort.bam SRR1571971.sam
samtools sort -@ 4 -O bam -o SRR1571972.sort.bam SRR1571972.sam

これで全てのマッピングファイルをソートしbamに変換できました。

StringTieでassembleする

マッピング→ソート→bam変換が完了しました。
今あるbamデータはRNA-seqの何千万とあるリードが、
ゲノムのどこに相当するものかという情報だけです。
これからこのデータをもとに遺伝子発現量を比較するわけですが、
そのためにはゲノムのどこに何の遺伝子があるのかを知り、
各リードは何の遺伝子の一部であったのか情報を付与する必要があります。

このためのプログラムが**「StringTie」**になります。

StringTieのインストール

以下のコマンドを実行

brew install stringtie

condaを使っている場合は、bioconda から入れられます。

conda install -c bioconda stringtie

アノテーションファイルの取得

次にアノテーションファイルを取得いたします。
リファレンスゲノムと同じデータベースからアノテーションを取ります。
今回はUCSCのmm10をリファレンスゲノムとしましたので、
UCSCからアノテーションファイルを取得します。

UCSCのTable Browserで、次のように設定して get output をクリックしてダウンロードします。

項目	選択する値
clade / genome	Mammal / Mouse
assembly	リファレンスゲノムに対応するもの（今回は mm10）
group	Genes and Gene Predictions
track	NCBI RefSeq
output format	GTF - gene transfer format
output filename	UCSC.mm10.gtf.gz
file type returned	gzip compressed

macOSですと、勝手にファイルが解凍されると思います。
もし解凍されていなかったら、以下のコマンドで解凍しましょう。

#gzipファイルの解凍
gzip -d UCSC.mm10.gtf.gz

StringTieによるアノテーション情報付与

stringtie SRR1571967.sort.bam -o SRR1571967.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571967 -A SRR1571967.gene.abundance.tab

StringTieコマンド説明

stringtie stringtieコマンド実行
SRR1571967.sort.bam samtoolsでsortしたbamファイルを指定
-o SRR1571967.gtf出力ファイル名を指定
-p 4CPUの使用するスレッド数を指定
-G UCSC.mm10.gtf先程UCSCからダウンロードしたアノテーションファイルを指定
-l SRR1571967ラベル名を指定
-A SRR1571967.gene.abundance.tab タブ区切り形式ファイルも出力させ、そのファイル名を指定

StringTieの詳しい説明につきましては、下記の公式マニュアルに詳しい説明があります。

外部リンクStringTie 公式マニュアル（出力形式の説明） ↗https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual#output

これを全てのsort.bamファイルに対して実行します。

stringtie SRR1571967.sort.bam -o SRR1571967.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571967 -A SRR1571967.gene.abundance.tab
stringtie SRR1571968.sort.bam -o SRR1571968.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571968 -A SRR1571968.gene.abundance.tab
stringtie SRR1571969.sort.bam -o SRR1571969.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571969 -A SRR1571969.gene.abundance.tab
stringtie SRR1571970.sort.bam -o SRR1571970.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571970 -A SRR1571970.gene.abundance.tab
stringtie SRR1571971.sort.bam -o SRR1571971.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571971 -A SRR1571971.gene.abundance.tab
stringtie SRR1571972.sort.bam -o SRR1571972.gtf -p 4 -G UCSC.mm10.gtf -l SRR1571972 -A SRR1571972.gene.abundance.tab

実行すると、作業ディレクトリに各サンプルの .gtf と .gene.abundance.tab が作成されます。ls *.gtf で6サンプル分そろっていることを確認してください。

StringTie実行後、各サンプルの.gtfファイルと.gene.abundance.tabファイルが作業ディレクトリに作成された状態

gtfファイルをmergeする

次にStringTieで作成された.gtfファイルを一つにまとめ、比較できるようにします。
各サンプルで同定された転写産物をまとめたgtfファイルが作成されます。

なぜこの工程が必要なのか、先に押さえておきます。StringTieはサンプルごとに独立して新規転写産物を同定するため、サンプルAの MSTRG.1 とサンプルBの MSTRG.1 が別の転写産物を指してしまいます。このままでは発現量を横に並べて比較できません。

mergeは全サンプルのGTFを統合し、共通の転写産物セットを作る工程です。この統合GTFを基準に各サンプルを測り直すことで、サンプル間で同じIDが同じ転写産物を指すようになります。

mergelist.txtファイルの作成

mergeさせるgtfファイル(stringtieでアノテーションをつけた)のpathを全て指定するtxtファイルを作成します。
一番ラクなのは全てのファイルを作業ディレクトリに入れておくことで、その場合はファイル名を並べるだけで済みます。
macであればテキストエディットで作成し、ファイル名をmergelist.txtにして保存しましょう。
作業ディレクトリにファイルがない場合は相対パスを一つづつ書けばOKです。

mergelist.txt の中身は次のようになります。

SRR1571967.gtf
SRR1571968.gtf
SRR1571969.gtf
SRR1571970.gtf
SRR1571971.gtf
SRR1571972.gtf

mergelist.txtに、mergeするgtfファイル名を1行ずつ並べた状態

⚠

注意

テキストエディットで作成する場合は、標準テキスト形式（.txt）で保存してください。リッチテキストのままだと拡張子が .rtf になり、StringTieが読み込めません。

mergeの実行

stringtie --merge -G UCSC.mm10.gtf -o stringtie_merged.gtf mergelist.txt

コマンドの説明

stringtie --mergestringtieでmergeを実行することを指定
-G UCSC.mm10.gtf アノテーションファイル（GTF）を指定
-o stringtie_merged.gtf mergelist.txt 出力されるファイル名を指定、mergeするファイルを指定したtxtファイル名を指定

mergeを実行し、stringtie_merged.gtfファイルが作成されればOKです。

Ballgown用ファイルの作成

作業ディレクトリにballgownというフォルダを作ります。
その後、

stringtie SRR1571967.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571967/SRR1571967.gtf

コマンド説明

stringtiestringtieを実行
SRR1571967.sort.bam sort.bamファイルを指定
-e -Gで指定されたリファレンス転写物にマッチした転写物のみを使用するオプション
-B Ballgown用ファイルを作成
-p 4 使用するCPUのスレッド数を指定
-G stringtie_merged.gtf 先ほどStringtieでmergeしたファイルを指定
-o ballgown/SRR1571967/SRR1571967.gtf 出力先とファイル名を指定

全てのファイルに対して実行

stringtie SRR1571967.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571967/SRR1571967.gtf
stringtie SRR1571968.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571968/SRR1571968.gtf 
stringtie SRR1571969.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571969/SRR1571969.gtf
stringtie SRR1571970.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571970/SRR1571970.gtf
stringtie SRR1571971.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571971/SRR1571971.gtf
stringtie SRR1571972.sort.bam -e -B -p 4 -G stringtie_merged.gtf -o ballgown/SRR1571972/SRR1571972.gtf

よくある質問

Q. StringTieはbrewとconda、どちらで入れればいい？ どちらでも構いません。Macで環境を分ける必要がなければ brew install stringtie が手軽です。解析環境をプロジェクトごとに分けたい場合や、Linux／Google Colabで作業する場合は conda install -c bioconda stringtie を使います。同じマシンに両方入れると which stringtie でどちらが呼ばれるか分からなくなるので、どちらかに統一するのがおすすめです。

Q. なぜgtfファイルをmergeする必要があるの？ StringTieはサンプルごとに独立して転写産物を同定するため、mergeしないとサンプル間で転写産物IDが対応しません。サンプルAの MSTRG.1 とサンプルBの MSTRG.1 が別物になり、発現量を横に並べて比較できなくなります。mergeで共通の転写産物セットを作り、それを基準に測り直すことで初めてサンプル間比較が成立します。

Q. stringtie --merge が動かない／エラーになる まず --merge のハイフンが2本になっているか確認してください。Webページからコピーすると、ハイフン2つが「—」（ダッシュ記号）に変換されてしまうことがあります。次に mergelist.txt のパスです。作業ディレクトリ以外にgtfファイルがある場合は、mergelist.txt に書いたパスがカレントディレクトリから見て正しいかを確認します。テキストエディットで作った場合は、リッチテキスト（.rtf）になっていないかも確認してください。

Q. -e と -B は何をしているの？ -B はBallgownが読み込む形式のテーブル（.ctab ファイル群）を出力させるオプションです。-e は -G で指定したGTFに含まれる転写産物だけを対象にし、新規の転写産物を推定しないようにします。Ballgownに渡す最終ステップでは、この2つをセットで使うのが定石です。統合GTFの転写産物セットに固定して各サンプルを測るという意味になります。

Q. samtools sortでメモリ不足になる -m オプションでスレッドあたりの使用メモリを指定できます（例：-m 1G）。既定は768MBで、-@ で指定したスレッド数の分だけ掛け算で消費される点に注意してください。-@ 4 -m 2G なら合計8GB使います。搭載メモリが少ない場合はスレッド数を減らすか -m を下げます。

Q. StringTieの出力（FPKM/TPM）をそのまま発現量として使っていい？ サンプル内で遺伝子間を比べる分には使えますが、サンプル間の発現差を統計的に検定するには専用のツールが必要です。このシリーズでは次の記事でBallgownを使います。カウントベースで解析したい場合はDESeq2やedgeRが標準的な選択肢になります。

次の記事ではballgownを用いて発現差解析を行っていきたいと思います。

◆

ポイント

ここまでの samtools → StringTie → merge → Ballgown用ファイル作成を、サンプルごとにコマンドを並べる代わりに一括実行する方法もあります。nf-core/rnaseq なら定量まで通しで走り、途中で失敗しても -resume で続きから再開できます。