STAR-RSEMによる発現量推定その１

私が普段使っているSTAR, RSEMを使った発現量推定法を紹介します。あまり最新のアップデート情報などをフォローできていないので、もっと良いやり方をご存知の方はご教示ください。

ここでは例題として、ヒトES細胞とHUVEC細胞をペアエンドで読んだサンプル（各２サンプル）の発現量を比較し、発現変動遺伝子を抽出するというケースを考えます。リファレンスゲノムはhg38、遺伝子のアノテーションはEnsemblのGRCh38です。
以下のサイトも併せて参考にしてください。

RSEM (A. thaliana, paired-end RNA-Seq) | STAR-RSEM による転写産物の発現量推定

サンプルのダウンロード

RNA-seqのfastqをダウンロードします。

$ fastq-dump --split-files --gzip SRR065504   # Myers_H1-hESC_cell_2x75_200_1
$ fastq-dump --split-files --gzip SRR065495   # Myers_H1-hESC_cell_2x75_200_2
$ fastq-dump --split-files --gzip SRR065509   # Myers_HUVEC_cell_2x75_200_1
$ fastq-dump --split-files --gzip SRR065524   # Myers_HUVEC_cell_2x75_200_2

これにはものすごく時間がかかるので、コマンド末に&をつけて同時に回すのがよいですね。それでも一晩くらいはかかりました。
時間がかかりすぎるので、手元にご自身のサンプルがある場合はそちらを代用してください。

ゲノム配列・遺伝子情報のダウンロード

UCSC hg38のゲノムをダウンロードする方法は以下を参考に。

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する - Palmsonntagmorgen

gtfファイルはEnsemblからGRCh38のデータをダウンロードし、一部修正します。方法は以下を参考にしてください。

Gene annotation データを用意する（gtf形式） - Palmsonntagmorgen

RSEM, STARのインストール

RSEMはGitHubからクローンして、makeでコンパイルしています。

$ git clone https://github.com/deweylab/RSEM.git
$ cd RSEM
$ make

STARはソースコードのアーカイブをダウンロード後解凍、コンパイルします。

$ wget https://github.com/alexdobin/STAR/archive/2.6.1d.tar.gz
$ tar -xzf 2.6.1d.tar.gz
$ cd STAR-2.6.1d/source
$ make STAR              # Linuxの場合
$ make STARforMacStatic  # Macの場合

RSEM, STARの両方にパスを通してください。RSEMはRSEMディレクトリ直下、STARは STAR-2.6.1d/bin/ 以下の Linux_x86_64_static または MacOSX_x86_64 ディレクトリ以下に実行ファイルがあります。
パスの通し方は以下を参考：

環境変数PATHの通し方 - Palmsonntagmorgen

STARのインデックスを作成する

rsem-prepare-reference コマンドを使います。これはRSEMの中でSTARを使うためのインデックス作成コマンドですが、STAR単体で動かす時のインデックスとしても使えます。

$ mkdir rsem-star-index  # index用ディレクトリの作成
$ rsem-prepare-reference --star \    
  -p 12 \      # 使用CPU数
  --gtf Homo_sapiens.GRCh38.94.addchr.gtf \ # 入力gtfファイル
  UCSC-hg38.fa \             # 入力ゲノムファイル
  rsem-star-index/UCSC-hg38  # 出力されるindexファイル名
Parsed 200000 lines
Parsed 400000 lines
Parsed 600000 lines
Parsed 800000 lines
Parsed 1000000 lines
Parsed 1200000 lines
Parsed 1400000 lines
Parsed 1600000 lines
Parsed 1800000 lines
Parsed 2000000 lines
Parsed 2200000 lines
Parsed 2400000 lines
Parsed 2600000 lines
Parsing gtf File is done!
UCSC-hg38.fa is processed!
Warning: Cannot extract transcript ENST00000361681's sequence since the chromosome it locates, chrMT, is absent!
Warning: Cannot extract transcript ENST00000361739's sequence since the chromosome it locates, chrMT, is absent!
Warning: Cannot extract transcript ENST00000361789's sequence since the chromosome it locates, chrMT, is absent!
Warning: Cannot extract transcript ENST00000361453's sequence since the chromosome it locates, chrMT, is absent!
...

おや…何やらWarningが出てきました。"chrMTなんて染色体はないよ"と怒られているようです。 EnsemblではミトコンドリアをMTと表記していますが、UCSCではchrMとしているので、ゲノムと遺伝子ファイルで染色体名が異なっているのですね。

大した数ではないので無視してもよさそうですが、ここではきちんとgtfの染色体名を変換しておきましょう。

$ sed -e 's/chrMT/chrM/g' Homo_sapiens.GRCh38.94.addchr.gtf \
> Homo_sapiens.GRCh38.94.addchr.new.gtf

sedというのはファイル内の文字列を置換してくれるコマンドで、ここではHomo_sapiens.GRCh38.94.addchr.gtfファイル内の"chrMT"を"chrM"に置換してHomo_sapiens.GRCh38.94.addchr.new.gtfに出力しています。 sedコマンドはこんな感じに使える便利コマンドなので、ぜひ覚えましょう。一括置換の際は関係ない部分文字列までヒットして置換されていないか確認が必要ですが、今回は大丈夫でしょう。
(例えば"chr1"で置換したりすると、"chr11,chr12, ..."が全部ヒットしてしまいます)

では新しいgtfファイルを使って再びインデックス作成にチャレンジしましょう。

$ rsem-prepare-reference --star -p 12 \
  --gtf Homo_sapiens.GRCh38.94.addchr.new.gtf \
  UCSC-hg38.fa \
  rsem-star-index/UCSC-hg38
...
Parsing gtf File is done!
UCSC-hg38.fa is processed!
206534 transcripts are extracted.
Extracting sequences is done!
Group File is generated!
Transcript Information File is generated!
Chromosome List File is generated!
Extracted Sequences File is generated!

rsem-preref rsem-star-index/UCSC-hg38.transcripts.fa 1 rsem-star-index/UCSC-hg38
Refs.makeRefs finished!
Refs.saveRefs finished!
rsem-star-index/UCSC-hg38.idx.fa is generated!
rsem-star-index/UCSC-hg38.n2g.idx.fa is generated!

STAR  --runThreadN 12  --runMode genomeGenerate  --genomeDir rsem-star-index  --genomeFastaFiles /home/Database/UCSC/hg38/genome.fa  --sjdbGTFfile Homo_sapiens.GRCh38.94.addchr.new.gtf  --sjdbOverhang 100  --outFileNamePrefix rsem-star-index/UCSC-hg38
...

今度はエラーなくきちんと通っています。一番下にSTARで始まるコマンドが記載されていますが、これがSTARのインデックスを作るコマンドです。つまりRSEMは内部でSTARを起動してSTARのインデックスを作成しているのです。

完了したら、インデックスが正しく生成されているか確認しましょう。

$ ls rsem-star-index/UCSC-hg38/
Genome                UCSC-hg38.ti              exonInfo.tab
SA                    UCSC-hg38.transcripts.fa  geneInfo.tab
SAindex               UCSC-hg38Log.out          genomeParameters.txt
UCSC-hg38.chrlist     chrLength.txt             sjdbInfo.txt
UCSC-hg38.grp         chrName.txt               sjdbList.fromGTF.out.tab
UCSC-hg38.idx.fa      chrNameLength.txt         sjdbList.out.tab
UCSC-hg38.n2g.idx.fa  chrStart.txt              transcriptInfo.tab
UCSC-hg38.seq         exonGeTrInfo.tab

長くなりましたので続きは次回。年内更新予定です。

STAR-RSEMによる発現量推定その2 - Palmsonntagmorgen