Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

記事一覧

解析環境構築

SSH公開鍵の生成・設定の方法

環境変数PATHの通し方

pyenvでPython環境を構築する

バイオインフォマティクスのためのpython環境構築方法を考える

GitHubからプログラムをダウンロード・インストール

データ生成

 

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

2bit genome を作成する

genome tableを作成する

LiftOver: BEDファイルを異なるbuildへ変換

 

fastqデータ取得

SRAからfastqを取得する

ENA,DDBJからfastqを取得する

 

ゲノムマッピング

Readをゲノムにマッピング (その1) (2017/12/19 追記あり)

Readをゲノムにマッピング (その2)

Readをゲノムにマッピング (その3) 圧縮ファイルを入力にする方法

 

マップデータ操作

SAMtoolsとリダイレクト

SAMtoolsワンライナー覚書

マッピング: CRAM形式を試す

 

ChIP-seq解析: DROMPA

DROMPA3: その1 インストール

DROMPA3: その2 parse2wig

DROMPA3: その3 ピーク抽出(peak calling)

DROMPA3: その4 マップリード分布の可視化その1

DROMPA3: その5 シェル変数を使う

DROMPA3: その6 ChIP/Input ratio 及び p値の可視化

DROMPA3: その7 -i オプション詳細

DROMPA3: その8 GVコマンドでのマクロな可視化

DROMPA3: その9 リードプロファイル

DROMPA3: その10 ヒートマップ

ChIP-seq解析: 品質評価

Library complexity (PCR bias)とは何か

S/N比の評価手法 その1

S/N比の評価手法 その2 Cross-correlation profile

S/N比の評価手法 その3 deepTools

S/N比の評価手法 その4 SSP

ChIP-seq解析: ピークを入力とする操作

BEDtoolsワンライナー覚書

2サンプル間ピーク比較

 

RNA-seq解析:

RNA-seqによる発現量解析

HISAT-StringTie-Ballgown を試してみよう

HISAT-StringTie-Ballgown を試してみよう その2

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

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

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

SSH公開鍵の生成・設定の方法

たかが公開鍵、されど公開鍵。「公開鍵 生成」でググるとたくさん出てくるんだけど、やり方が色々ありすぎて新人に適当に検索させるとあれこれ迷ってしまうので、ここにまとめておきます。

参考記事

公開鍵とは

サーバやGitHubなどのリモートサービスにsshでログインする際、通常のパスワード形式を用いると、パスワードが何らかの理由で漏洩した場合に他人にログインされる危険があります。そしてパスワードは常に漏洩の危険があります。
これに対して公開鍵形式の場合は、パスワードに関連付けられた「鍵と錠前」を生成します。パスワードを知っており、かつ鍵を持っている人のみがログイン可能であるため、よりセキュアになります。
この「鍵と錠前」のことを秘密鍵と公開鍵と呼びます。この二つはセットになっており、同時に生成されます。
公開鍵は錠前なので、ログインしたいサーバ上に置きます。公開鍵は他人に見られてもかまいません。
秘密鍵は鍵なので、自分だけが持っておくものです。他人に見られてはいけません。  

公開鍵の作成方法

ここではLinuxでのやり方に統一します。暗号化方式はRSAを用います。Windows10の場合はWSL、Macの場合はMacターミナル上で作業してください。
Windows10でWSLを未インストールの方は以下を参考にインストールしてください。

【WSL入門】第1回 Windows 10標準Linux環境WSLを始めよう:ITの教室 - @IT

既存の公開鍵が存在するかチェック

公開鍵はデフォルトでは~/.sshに生成・保管されます(~/はホームディレクトリ)。 既に作成済みの公開鍵を上書きしないよう、まずは公開鍵を持っているか調べてみましょう。

$ ls ~/.ssh

RSA形式の場合、秘密鍵id_rsa、公開鍵はid_rsa.pubという名前で保存されています。これらのファイルが無いか、または~/.sshディレクトリ自体が存在しない場合は、公開鍵はまだ生成されていません。

公開鍵・秘密鍵の生成

以下のssh-keygenコマンドで公開鍵と秘密鍵が生成されます。 鍵長はデフォルトの2048で良いのですが、ここではよりセキュアな4096を指定することにします。

$ ssh-keygen -t rsa -b 4096
Generating public/private rsa key pair.
Enter file in which to save the key (/home/rnakato/.ssh/id_rsa):
Enter passphrase (empty for no passphrase):
Enter same passphrase again:
Your identification has been saved in /home/rnakato/.ssh/id_rsa.
Your public key has been saved in /home/rnakato/.ssh/id_rsa.pub.
The key fingerprint is:
SHA256:6KfApJsjSs/tj6FPQioHS2MWk4fYlnJXaCp9uowg8Ac rnakato@RyuHomePC
The key's randomart image is:
+---[RSA 4096]----+
|     ..          |
|..o.o.           |
|o*=+.            |
|o+E..  .         |
|oB =. . S        |
|*o*+..           |
|==o+oo. .        |
|=.*o=.oo         |
|o.o=o=o.         |
+----[SHA256]-----+

パスフレーズの入力を求められますが、ここで入力したものがサーバにログインする際のパスワードになります。 サーバ上にあるアカウントのパスワード(sudoなどで利用するもの)とは別であることに注意してください。

パスワードは5文字未満だとエラーになります。しかし空パスワードはOKのようです(謎ですが…)。 当然ですが、空パスワードや容易に推測可能な単純なパスワードは危険ですので、十分複雑なパスワードを指定してください。

終了後再び$ ls ~/.sshを実行し、id_rsaid_rsa.pubが生成されていれば成功です。

公開鍵を使ったログイン

公開鍵の設定自体はサーバサイドで行われるものなので、ユーザの設定は不要です。 サーバ管理者に公開鍵を送付すれば準備OKです。

Linuxシェル上だと、以下のような感じでログインするでしょう。

$ ssh rnakato@hostname

何も指定しなかった場合、秘密鍵~/.ssh/id_rsaが自動的に利用されます。秘密鍵の場所が違ったり、ファイル名が異なる場合は以下のように明示的に指定しましょう。

$ ssh rnakato@hostname -i <秘密鍵>

公開鍵は複数のサーバで共用可能です。その場合すべて同じ秘密鍵とパスワードを用いてログインすることが可能です。 逆に秘密鍵を無くしたりパスワードを忘れて作り直したという場合は、サーバ上の公開鍵もすべて新しいものに更新する必要があります。

秘密鍵パーミッション

秘密鍵は「誰にも見られない」ことが暗黙に想定されているので、秘密鍵パーミッションがオープン過ぎる場合、以下のようなエラーになり、ログインできません。

$ ssh rnakato@hostname
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@         WARNING: UNPROTECTED PRIVATE KEY FILE!          @
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Permissions 0777 for '/home/rnakato/.ssh/id_rsa' are too open.
It is required that your private key files are NOT accessible by others.
This private key will be ignored.
bad permissions: ignore key: /home/rnakato/.ssh/.ssh/id_rsa
Permission denied (publickey).

この場合は以下のようにパーミッションを400に変更すればログインできるようになります。

$ chmod 400 ~/.ssh/id_rsa

マッピング: CRAM形式を試す

マップデータの形式にはSAM, BAMの他にCRAMという形式があります。

https://www.ga4gh.org/news/cram-compression-for-genomics/

CRAMはBAMと比べて更に高圧縮率だそうです。
今まであまり使う機会が無かったのですが、Twitterで Ewan Birneyさんが強く推していたので、今日はCRAMを試してみようと思います。

SAMtoolsのインストール

今はSAMtools でCRAMも扱えるんですね。便利ですね。
SAMtoolsのインストールは↓を参照してください。

SAMtoolsとリダイレクト - Palmsonntagmorgen

手元のSAMtoolsのバージョンは1.9です。

$ samtools --version
samtools 1.9
Using htslib 1.9
Copyright (C) 2018 Genome Research Ltd.

ゲノムデータの準備

CRAMの変換にはゲノム配列データが必要になります。ゲノム配列作成は以下を参照してください。

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

また、ゲノムに indexが付加されていることが推奨なので、faidxコマンドでindexを作成しておきます。

$ samtools faidx genome_hg19.fa

genome_hg19.fa.fai が作成されます。

bamのダウンロード

bamファイルは何でもいいのですが、今回は以下からダウンロードした GSM733643_hg19_wgEncodeBroadHistoneK562Pol2bStdAlnRep2.bam を使います。

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM733643

(リンク先下部の"Supplementary file" からBAMがダウンロードできます)

BAM -> CRAMの変換

SAMtoolsを使ってBAMからCRAMへ変換してみます。上でindexを作成しましたが、-Tで指定するのはゲノム配列であることに注意してください。 ここではtime コマンドを追加して実行時間を計測しています。
bamのファイル名が長いと見づらいので、ここでは仮にPol2.bamとします。

$ time samtools view -C Pol2.bam -T genome_hg19.fa > Pol2.cram

real    1m46.909s
user    1m37.133s
sys 0m5.343s

変換に時間がかかるという噂がありましたが、シングルコアでも1分ほどで終わりました。

$ ls -lh
合計 906M
-rwxrwxrwx 1 rnakato rnakato 618M  410 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 288M  410 17:15 Pol2.cram

確かに、ファイルサイズが半分以下になっています。情報量が同じでサイズが半分以下なのは良いですね。

CRAM -> BAMの変換

ファイルが可逆変換か確認するために、CRAMからBAMに変換してみます。 元のBAMファイルと同一のファイルが生成されることを期待しています。

$ time samtools view -b Pol2.cram > Pol2.bam2
real    1m19.342s
user    1m17.821s
sys 0m1.135s

BAMに戻す時にはゲノムを指定する必要はないようです。

$ ls -lh
合計 1.6G
-rwxrwxrwx 1 rnakato rnakato 618M  410 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 683M  410 17:17 Pol2.bam2
-rwxrwxrwx 1 rnakato rnakato 288M  410 17:15 Pol2.cram

何故だか元bamファイルよりファイルサイズが大きくなっています。

3つともsamファイルに変換して調べたところ、Pol2.bam2はPol2.bamと比べてカラムが2つ増えていました。詳細をきちんと確認していませんが、このあたりは元SAM・BAMファイルをどのツールで生成したかによって結果が変わるかもしれません。 なお、Pol2.bam2とPol2.cramは同一でした。

CRAMのソート

$ samtools sort Pol2.cram -O cram -@ 12 > Pol2.sort.cram

問題なくソートできました。 -O cramで出力をCRAMに指定するのを忘れるとbamになってしまうので注意。

bowtie2からのリダイレクト

最後に、bowtie2でマッピングした結果をリダイレクトして直接ソート済CRAMを生成することにチャレンジ。

まず元fastqのダウンロード(何でもよいです)

$ fastq-dump SRR227359

bowtie2でマッピングしてみます。bowtie2のデフォルト出力はSAM形式です。

# SAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq > SRR227359.sam
# sorted BAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort > SRR227359.sort.bam
# sorted CRAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort -O cram - -T genome_hg19.fa > SRR227359.sort.cram
[bam_sort_core] merging from 5 files and 1 in-memory blocks...
[E::cram_get_ref] Failed to populate reference for id 0
samtools sort: failed to write header to "-"

CRAMはsamtools sortに直接投げるとエラーになりました。 修正して以下のように2段階のリダイレクトにします。

# sorted CRAM(修正)
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq \
  | samtools view -C - -T genome_hg19.fa \
  | samtools sort -O cram \
  >  SRR227359.sort.cram

今度はうまく行きました。

BAMの時と比較しても実行時間はそれほど変わりません。計算時間的な意味でのオーバーヘッドはほとんど気になりません。

考察

SAMtoolsを使うとCRAMはほとんどストレスなく取り扱えます。ゲノム配列を指定するのがちょっと面倒なくらいです。 計算時間はほとんど変わらず、ファイルサイズは小さくなるので、有用性は高いです。
あとは各種解析ツールが入力にCRAM形式をどの程度受け付けるかによりますね。今のところは対応しているツールはそれほど多くないと思うので、「しばらく使う予定は無いが消すことはできない」というような大量のBAMファイルがある場合にCRAM形式にしておくと容量が削減できてうれしい、というところでしょうか。

あとはpaired-endの時にどうなるのかなどチェックした方がいいかもしれません。

【ご報告】PIになりました

この度4月1日付で東京大学定量生命科学研究所の講師となりました。
大規模生命情報解析研究分野という名前で研究室を主宰します。
研究室HP: http://www.iam.u-tokyo.ac.jp/nakatolab/index.html

年度末営業と研究室異動のばたばたでこのブログも全く更新できていませんでしたが、これからも引き続きゆるく続けていきたいと思っています。
なお当研究室では学生を募集していますので、こういったNGS解析に興味がありそうな学生さんがいましたら是非ご連絡くださいませ。

今後とも宜しくお願い致します。

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

前回の記事↓の続きです。

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

前回のコマンドを最後まで完了すると、starディレクトリの中にstar/Myers_HUVEC_cell_2x75_200_1.<genes|isoforms>.results のようなファイルがサンプルごとに生成されているのがわかります。 *.genes.results*.isoforms.resultsはそれぞれ遺伝子単位、transcript単位の発現量(count, TPM, FPKM)が出力されます。 transcript単位の出力は基本的には精度が十分でないことが多く、ノイズがまぎれこみやすいので、特に利用目的がないのであれば遺伝子単位のものを利用しましょう。

RSEMの出力

では、ファイルの中をのぞいてみましょう。

$ head star/Myers_HUVEC_cell_2x75_200_1.genes.results
gene_id transcript_id(s)    length  effective_length    expected_count  TPM FPKM
ENSG00000000003 ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2146.69 1926.19 1086.00 40.20   34.28
ENSG00000000005 ENST00000373031,ENST00000485971 1046.60 826.17  14.00   1.21    1.03
ENSG00000000419 ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1003.23 782.75  355.00  32.34   27.57
ENSG00000000457 ENST00000367770,ENST00000367771,ENST00000367772,ENST00000423670,ENST00000470238 3672.40 3451.91 112.23  2.32    1.98
ENSG00000000460 ENST00000286031,ENST00000359326,ENST00000413811,ENST00000459772,ENST00000466580,ENST00000472795,ENST00000481744,ENST00000496973,ENST00000498289 2167.45 1946.97 501.77  18.38   15.67
ENSG00000000938 ENST00000374003,ENST00000374004,ENST00000374005,ENST00000399173,ENST00000457296,ENST00000468038,ENST00000475472 2021.00 1800.50 9.00    0.36    0.30
ENSG00000000971 ENST00000359637,ENST00000367429,ENST00000466229,ENST00000470918,ENST00000496761,ENST00000630130 2587.83 2367.38 0.00    0.00    0.00
ENSG00000001036 ENST00000002165,ENST00000367585,ENST00000451668 2285.18 2064.69 1079.67 37.29   31.79
ENSG00000001084 ENST00000229416,ENST00000504353,ENST00000504525,ENST00000505197,ENST00000505294,ENST00000509541,ENST00000510837,ENST00000513939,ENST00000514004,ENST00000514373,ENST00000514933,ENST00000515580,ENST00000616923,ENST00000643939,ENST00000650454 2327.40 2106.94 773.00  26.16   22.30

Ensemblからダウンロードしたgtfファイルを用いているので、gene_idとtranscript_idは当然Ensemblのものになっています。 effective length はpoly(A) tailを除いた転写物の長さだそうです。 expected_count は(isoformごとに振り分けられた)マップリード数、FPKMはリード数を全マップリード数と遺伝子長で正規化したスコアで、TPMはFPKMを更にFPKMの総和で正規化した値になります。

発現量の評価にはTPM、edgeRの入力にはcountを使おう

遺伝子発現量を直接見たい場合、例えば box plot や scatter plot を描く場合には TPMを使っておけば基本的に間違いありません。 一方、edgeRやDESeq2のような発現変動比較ツールにはTPMではなくcountデータを利用することが推奨されます(マニュアル参照)。何故なら、edgeRやDESeq2の中で用いている正規分布や負の二項分布はraw countに対してモデル化するものであり、遺伝子長などで正規化してしまうとこの分布を満たさなくなってしまうおそれがあるからです。全ての正規化はedgeRやDESeq2の中で行われますので、前もって正規化しておく必要はありません。更に言えば、サンプル間比較は遺伝子ごとに行われますので、基本的に遺伝子長の正規化は必要ないのです。

なお、DESeq2は整数値しか扱えませんので、DESeq2を利用する時にはcountデータの少数を切り捨てにしておきましょう。

発現量データをマージ (count)

それでは、4つのサンプルの発現量データをマージしてみましょう。RSEMの中にある rsem-generate-data-matrixを使います。 ここでは遺伝子ごとのファイルを用いています。

$ RSEM/rsem-generate-data-matrix   \
    star/Myers_H1-hESC_cell_2x75_200_1.genes.results  \
    star/Myers_H1-hESC_cell_2x75_200_2.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_1.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_2.genes.results  \
   > GeneExpressionMatrix.tsv

得られた結果を見てみましょう。

$ head GeneExpressionMatrix.tsv 
    "star/Myers_H1-hESC_cell_2x75_200_1.genes.results"   "star/Myers_H1-hESC_cell_2x75_200_2.genes.results"   "star/Myers_HUVEC_cell_2x75_200_1.genes.results" "star/Myers_HUVEC_cell_2x75_200_2.genes.results"
"ENSG00000000003"    1086.00 2327.00 1086.00 1891.00
"ENSG00000000005"    14.00   34.00   14.00   0.00
"ENSG00000000419"    355.00  755.00  355.00  875.00
"ENSG00000000457"    112.23  212.28  112.23  145.02
"ENSG00000000460"    501.77  878.72  501.77  201.98
"ENSG00000000938"    9.00    15.00   9.00    1.00
"ENSG00000000971"    0.00    0.00    0.00    7301.00
"ENSG00000001036"    1079.67 2036.77 1079.67 1889.43
"ENSG00000001084"    773.00  1263.00 773.00  608.00

行が遺伝子、列がサンプルの行列データになっていることがわかります。タブ区切りのテキストファイルなので、ExcelやR, Pythonなど好きなもので分析することが可能です。

発現量データをマージ (TPM)

上記で出力されているデータはカウントデータです。rsem-generate-data-matrixはカウントデータにしか対応していません(手抜き?)。

TPM や RPKM のマトリックスを出力したい場合は、僕が作成した rsem-generate-data-matrix-modified を使ってください。 rsem-generate-data-matrix-modifiedは、↓の記事で紹介しているscript_rnakato の中にあります。

GitHubからプログラムをダウンロード・インストール - Palmsonntagmorgen

ではやってみましょう。

$ script_rnakato/rsem-generate-data-matrix-modified \
    TPM \    # count または TPM または FPKM を指定
    star/Myers_H1-hESC_cell_2x75_200_1.genes.results  \
    star/Myers_H1-hESC_cell_2x75_200_2.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_1.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_2.genes.results  \
    > GeneExpressionMatrix.TPM.tsv
# 結果を表示
$ head GeneExpressionMatrix.TPM.tsv 
    Myers_H1-hESC_cell_2x75_200_1.genes.results Myers_H1-hESC_cell_2x75_200_2.genes.results Myers_HUVEC_cell_2x75_200_1.genes.results   Myers_HUVEC_cell_2x75_200_2.genes.results
ENSG00000000003 40.20   41.67   40.20   39.30
ENSG00000000005 1.21    1.91    1.21    0.00
ENSG00000000419 32.34   34.32   32.34   39.18
ENSG00000000457 2.32    2.47    2.32    1.88
ENSG00000000460 18.38   15.50   18.38   4.85
ENSG00000000938 0.36    0.30    0.36    0.02
ENSG00000000971 0.00    0.00    0.00    93.63
ENSG00000001036 37.29   34.21   37.29   36.32
ENSG00000001084 26.16   25.42   26.16   11.66

無事TPMマトリックスを得ることができました。ちなみにgene id の" も必要ないので除去しています。

おわりに

これで無事、fastqファイルから遺伝子発現量マトリックスを生成できるようになりました。 あとはこのtsvファイルをBioconductorなどで解析すればばっちりですね。 細かい点では、gene idを遺伝子シンボルに変更したいなどの要求があると思いますが、機会があれば追記します。

【Ubuntu 18.04】 Rのバージョンを3.5.2にアップグレード

最近知ったツールがRの3.5以上を要求しているのだが、手元のUbuntu 18.04はRが3.4.4 だったので、3.5にアップグレードしました。

以下備忘録。 やってることは、aptで参照されるレポジトリにR3.5を含むレポジトリを追加した上で改めてRをインストールするというものです。シンプル。

既存のRをアンインストール(必要ないかも?)

$ sudo apt remove --purge r-base r-base-core r-recommended
$ sudo apt autoremove

aptのレポジトリを追加

$ sudo add-apt-repository 'deb https://cloud.r-project.org/bin/linux/ubuntu bionic-cran35/'
$ sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E084DAB9
$ sudo apt update

Rをインストール

$ sudo apt install r-base r-base-core r-recommended r-base-dev

Rを起動して

update.packages()  # 既存パッケージのアップデート

余談

最初追加したレポジトリをbionic-cran35でなくxenial-cran35 にしていたせいで libpng12-0がないですとか色々依存関係のエラーになってしまった。Ubuntuのバージョンはちゃんと確認しよう。

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

前回の記事↓の続きです。

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

Stranded/Unstranded RNA-seq

RNA-seqにはunstrandedとstranded の二種類があります。unstrandedの場合はmRNAに対して半々の確率で順鎖と逆鎖が読まれ、stranded の場合は逆鎖のみが読まれます(illumina TruSeq kitの場合)。 最近のRNA-seqはほぼ全てStranded RNA-seqになっていると思いますが、古いデータはunstrandedである場合があります。 それぞれマッピング時のパラメータが異なりますので、事前にチェックしておく必要があります。

どちらかわからんという場合は以下のように簡単にマッピングすることで確かめることができます。 indexは既知のmRNA配列から作成したindexとします。

$ index=NCBI-H_sapiens-RNA
$ bowtie $index <(zcat SRR065495_1.fastq.gz) -p12 | cut -f2 | sort | uniq -c
# reads processed: 50170737
# reads with at least one reported alignment: 34437131 (68.64%)
# reads that failed to align: 15733606 (31.36%)
Reported 34437131 alignments
17168386 +
17268745 -

詳細は割愛しますが、この結果はマップされた34437131リードのうち、17168386リードが順鎖、17268745が逆鎖にマップされたことを表しています。つまりほぼ半々なので、今回のデータはunstranded RNA-seqとわかります。

STARでマッピング

それでは改めてSTARでマッピングします。 RSEMコマンド中でSTARのマッピングもできるのですが、ログの出力が貧弱なため、私は独立でSTARを動かしてマッピングしています。
STARのコマンドは以下のようになります。

STAR --outSAMtype BAM SortedByCoordinate \   # sort済BAM を出力
     --quantMode TranscriptomeSAM \    # RSEMでの発現量推定に必要なオプション
     --runThreadN 12 \            # 使用CPU数
     --outSAMattributes All \     # 出力BAMファイルに含まれる情報を最大に
     --readFilesCommand zcat \    # gzip圧縮のfastqを入力にする場合のみ必要
     --genomeDir rsem-star-index/UCSC-hg38 \     # indexファイルを指定
     --readFilesIn SRR065504_1.fastq.gz SRR065504_2.fastq.gz \  # 入力fastq
     --outSAMstrandField intronMotif \    #   unstranded RNA-seqの場合のみ必要
     --outFileNamePrefix star/Myers_H1-hESC_cell_2x75_200_1.  # 出力ファイル名(最後のピリオドを忘れず)

STARは相当メモリを食いますので注意してください(十数GB程度)。

4サンプルに対して上記を実行しますが、入力ファイル名と出力ファイル名が異なるので、今回はfor文は使いません。 色々やり方はありますが、ここでは以下のようにパラメータ列をまとめることにします。

starparam="--genomeLoad NoSharedMemory \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode TranscriptomeSAM \
     --runThreadN 12 \
     --outSAMattributes All \
     --readFilesCommand zcat \
     --genomeDir rsem-star-index/UCSC-hg38 \
     --outSAMstrandField intronMotif"

STAR $starparam \
     --readFilesIn SRR065504_1.fastq.gz SRR065504_2.fastq.gz \
     --outFileNamePrefix star/Myers_H1-hESC_cell_2x75_200_1.
STAR $starparam \
     --readFilesIn SRR065495_1.fastq.gz SRR065495_2.fastq.gz \
     --outFileNamePrefix star/Myers_H1-hESC_cell_2x75_200_2.
STAR $starparam \
     --readFilesIn SRR065504_1.fastq.gz SRR065504_2.fastq.gz \
     --outFileNamePrefix star/Myers_HUVEC_cell_2x75_200_1.
STAR $starparam \
     --readFilesIn SRR065524_1.fastq.gz SRR065524_2.fastq.gz \
     --outFileNamePrefix star/Myers_HUVEC_cell_2x75_200_2.

starディレクトリの中に以下のようにファイルが生成されていれば成功です。

$ ls star
Myers_H1-hESC_cell_2x75_200_1.Aligned.sortedByCoord.out.bam
Myers_H1-hESC_cell_2x75_200_1.Aligned.toTranscriptome.out.bam
Myers_H1-hESC_cell_2x75_200_1.Log.final.out
Myers_H1-hESC_cell_2x75_200_1.Log.out
Myers_H1-hESC_cell_2x75_200_1.Log.progress.out
Myers_H1-hESC_cell_2x75_200_1.SJ.out.tab
Myers_H1-hESC_cell_2x75_200_2.Aligned.sortedByCoord.out.bam
Myers_H1-hESC_cell_2x75_200_2.Aligned.toTranscriptome.out.bam
Myers_H1-hESC_cell_2x75_200_2.Log.final.out
Myers_H1-hESC_cell_2x75_200_2.Log.out
Myers_H1-hESC_cell_2x75_200_2.Log.progress.out
Myers_H1-hESC_cell_2x75_200_2.SJ.out.tab
Myers_HUVEC_cell_2x75_200_1.Aligned.sortedByCoord.out.bam
Myers_HUVEC_cell_2x75_200_1.Aligned.toTranscriptome.out.bam
Myers_HUVEC_cell_2x75_200_1.Log.final.out
Myers_HUVEC_cell_2x75_200_1.Log.out
Myers_HUVEC_cell_2x75_200_1.Log.progress.out
Myers_HUVEC_cell_2x75_200_1.SJ.out.tab
Myers_HUVEC_cell_2x75_200_2.Aligned.sortedByCoord.out.bam
Myers_HUVEC_cell_2x75_200_2.Aligned.toTranscriptome.out.bam
Myers_HUVEC_cell_2x75_200_2.Log.final.out
Myers_HUVEC_cell_2x75_200_2.Log.out
Myers_HUVEC_cell_2x75_200_2.Log.progress.out
Myers_HUVEC_cell_2x75_200_2.SJ.out.tab

マッピングスタッツを確認してみましょう。*.Log.final.outにスタッツが記載されています。

$ cat star/Myers_HUVEC_cell_2x75_200_2.Log.final.out 
                                 Started job on |   Dec 21 14:48:27
                             Started mapping on |   Dec 21 14:50:50
                                    Finished on |   Dec 21 15:09:14
       Mapping speed, Million of reads per hour |   156.28

                          Number of input reads |   47926048
                      Average input read length | 150
                                    UNIQUE READS:
                   Uniquely mapped reads number |   37351930
                        Uniquely mapped reads % |   77.94%
                          Average mapped length |   147.79
                       Number of splices: Total |   20037714
            Number of splices: Annotated (sjdb) |   19925915
                       Number of splices: GT/AG |   19886116
                       Number of splices: GC/AG |   121436
                       Number of splices: AT/AC |   15391
               Number of splices: Non-canonical |   14771
                      Mismatch rate per base, % |   1.13%
                         Deletion rate per base |   0.01%
                        Deletion average length |   1.44
                        Insertion rate per base |   0.00%
                       Insertion average length |   1.34
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |   3004986
             % of reads mapped to multiple loci |   6.27%
        Number of reads mapped to too many loci |   9073
             % of reads mapped to too many loci |   0.02%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |   0.00%
                 % of reads unmapped: too short |   14.25%
                     % of reads unmapped: other |   1.53%
                                  CHIMERIC READS:
                       Number of chimeric reads |   0
                            % of chimeric reads |   0.00%

マップリード数のほか、スプライスサイトの情報なども載っていて便利です。

RSEMで発現量推定

STARで出力された*.Aligned.toTranscriptome.out.bamを使ってRSEMで発現量推定をします。今度はfor文が使えます。

for prefix in Myers_H1-hESC_cell_2x75_200_1 \
              Myers_H1-hESC_cell_2x75_200_2 \
              Myers_HUVEC_cell_2x75_200_1 \
              Myers_HUVEC_cell_2x75_200_2
do
    rsem-calculate-expression \
               --paired-end --alignments \
               --estimate-rspd \
               --forward-prob 0.5 \  # unstrandedの場合は0.5
               --no-bam-output \
               -p 12 \
               star/$prefix.Aligned.toTranscriptome.out.bam \
               rsem-star-index/UCSC-hg38/UCSC-hg38 \
               star/$prefix
done

indexの指定方法がSTARの時と異なることに注意してください。完了するとstarディレクトリの中に*.genes.results*.isoforms.resultsいうファイルが生成され、それぞれ遺伝子単位、transcript単位の発現量(count, TPM, FPKM)が出力されます。 発現量としてはTPMを使うのが最も望ましいでしょう。
ここではstarディレクトリの中に発現量も出力していますが、場合によっては別にrsemディレクトリを作って発現量をそちらに格納する方がSTARとRSEMの区別がはっきりして良いかもしれません。

また長くなったので、ファイルのマージ方法は次回に回します。 次回は年明けの予定です。