2021/05/25 追記: Ubuntuでのインストール方法を更新
前回の話とのつながりとしては、STARでRNA-seqのリードをリファレンスゲノムにマッピングしたという前提である。
さて、このマッピングしただけのデータ(SAMファイル)では統計的な解析に向かないという問題点がある。もちろん、これをIGVなどのゲノムブラウザに読み込ませてどこにリードが貼り付いているかを”目で見て”確かめることはできる。しかしRNA-seqで網羅的に遺伝子の発現解析を行うような場合はそのようなことをしていたらきりがない。
そこでゲノム上の「遺伝子」にどれだけリードが貼り付いたかどうかをカウント、集計するためのソフトウエアが存在する。それが今回紹介するFeatureCountsである。
インストール
Ubuntuであればaptでインストールすることが可能。以前まではfeatureCountsというパッケージ名だったが、2021年5月25日現在、subreadに変更されている。
$ sudo apt install subread
使い方
$ featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2]
基本的には上記の使い方になる。今回はペアエンドの解析に話を絞る。今現在RNA-seq自体の価格が安くなり、4万円を切るところが出てきていることからも今後シングルエンドで解析をする機会はまずないだろう。そのため今回のマッピングされたデータもペアエンドのものである前提である。
ペアエンドの解析ではゲノム上に貼り付いたリードをどのように処理するかのオプションを指定する必要がある。例えば私は以下のコードによってカウントしている。
$ featureCounts -p -B -T 24 -t exon -g gene_id -a RhiBivv1.1.gtf -o FeatureCounts.txt Aligned.sort.bam
-pオプションでは貼り付いたリードをカウントするのではなく、フラグメント(読み取った元のRNA断片)をカウントする。
また、もともとのSAM(BAM)ファイルには片方のリードの相方となるリードの情報(ペアの情報)が含まれている。更にその両方のペアリードがきちんとマッピングされているかの情報も含まれている(SAMファイルの2列目)。そのためその情報を利用して更に厳しく、「両方のリードがマッピングしている場合のみカウントする」というオプションもあり、それが -B オプションである。
ペアエンドのカウントに関連のあるオプションは以下の通り。
-p If specified, fragments (or templates) will be counted instead of reads. This
option is only applicable for paired-end reads.
-B Only count read pairs that have both ends aligned.
-P Check validity of paired-end distance when counting read pairs. Use -d and -D to
set thresholds.
-d <int>
Minimum fragment/template length, 50 by default.
-D <int>
Maximum fragment/template length, 600 by default.
-C Do not count read pairs that have their two ends mapping to different chromosomes
or mapping to same chromosome but on different strands.
--donotsort
Do not sort reads in BAM/SAM input. Note that reads from the same pair are required
to be located next to each other in the input.
その他のオプションは、例えば -T であれば使用するスレッド数。
-t はGFF/GTFファイルから使用する領域(デフォルトでexon)。
-g はGFF/GTFファイルの最後のカラム(列)のどの情報を利用するか、今回はgene_id で遺伝子のidを取得している。これは各GFF/GTFファイルに合わせてやる必要があり、形式が異なるとエラーを吐くことが多い。
-o はoutputファイル名。
最後はインプットファイルだが、復数のデータがある場合スペース区切りで入力してあげる。
例としては以下の通り。
$ featureCounts -p -B -T 24 -t exon -g gene_id -a RhiBivv1.1.gtf -o FeatureCounts.txt Aligned.sort1.bam Aligned.sort2.bam Aligned.sort3.bam
これでアウトプットには復数のカウントデータがまとめて出力されるため、その先のDEG解析などには有効に使うことができる。
引用文献
Liao, Y., Smyth, G. K., & Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7), 923-930.
