公開日: 最終更新日:

FeatureCountsでマッピングしたペアエンドデータのカウントを行う

前回の話とのつながりとしては、STARでRNA-seqのリードをリファレンスゲノムにマッピングしたという前提である。

さて、このマッピングしただけのデータ(SAMファイル)では統計的な解析に向かないという問題点がある。もちろん、これをIGVなどのゲノムブラウザに読み込ませてどこにリードが貼り付いているかを”目で見て”確かめることはできる。しかしRNA-seqで網羅的に遺伝子の発現解析を行うような場合はそのようなことをしていたらきりがない。

そこでゲノム上の「遺伝子」にどれだけリードが貼り付いたかどうかをカウント、集計するためのソフトウエアが存在する。それが今回紹介するFeatureCountsである。

インストール

Ubuntuであればaptでインストールすることが可能。

$ sudo apt install featureCounts

使い方

 $ 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. Bioinformatics30(7), 923-930.