私が普段使っているSTAR, RSEMを使った発現量推定法を紹介します。 あまり最新のアップデート情報などをフォローできていないので、もっと良いやり方をご存知の方はご教示ください。
ここでは例題として、ヒトES細胞とHUVEC細胞をペアエンドで読んだサンプル(各2サンプル)の発現量を比較し、発現変動遺伝子を抽出するというケースを考えます。
リファレンスゲノムは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
長くなりましたので続きは次回。年内更新予定です。