Palmsonntagmorgen

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

記事一覧

解析環境構築

環境変数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ワンライナー覚書

 

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

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の区別がはっきりして良いかもしれません。

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

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

私が普段使っている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

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

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

Gene annotation データを用意する(gtf形式)

RNA-seq解析の記事を書こうとしたらgtfファイルの部分が長くなり過ぎたので単独記事にしました。

Gene annotation

既知遺伝子の情報を記載するファイルにはいくつかの形式があり、gtf形式はそのひとつです。 よく似たgff形式というものもありますが、微妙に異なります。

gtf形式の詳細については以下の記事を参考にしてください(丸投げ)。

GTFファイル | 遺伝子アノテーションファイルの処理

GTFファイルのダウンロード

遺伝子情報はダウンロード元のデータベースによってIDや情報の詳細度がかなり違います。 ここでは私がよく使うUCSC genome browserとEnsemblからダウンロードしてみたいと思います。

UCSC genome browser

Table browserで遺伝子情報を指定し、"get output" ボタンを押します。RefseqのHuman build hg38 であれば、以下のようになります。

Table Browser

Ensembl

Ensemblは各種データを↓のFTP siteからダウンロード可能です。

http://asia.ensembl.org/info/data/ftp/index.html

Human のGRCh38 (hg38に対応)データであれば、以下のコマンドでダウンロードできます。

$ wget ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.gtf.gz

scaffoldなどの遺伝子情報を含んだgtfファイルも別にあるますが、普通は使いません。

ダウンロードした二つのファイルを解凍後、サイズを比較してみましょう。

$ gunzip *.gtf.gz
$ ls -lh *.gtf
-rwxrwxrwx 1 rnakato rnakato 1.1G Dec 13 11:17 Homo_sapiens.GRCh38.94.gtf
-rwxrwxrwx 1 rnakato rnakato  50M Dec 13 11:13 UCSC.hg38.gtf

Ensemblの方がファイルサイズが圧倒的に大きいですね。解凍すると1Gbを超えています。Ensemblの方がprocessed transcriptやpseudogeneなどの情報が多く、アノテーションも詳細に記述されていますので、遺伝子情報をもとに詳細な解析をしたい場合はEnsemblを使う方が良いと思います。

染色体名の違い

UCSCEnsembl はgene IDが違いますが、染色体名も違います。UCSCの方は 1番染色体は"chr1"となっていますが、Ensemblは単に"1"となっています。UCSCのgenomeファイルを用いて解析している場合、この違いを考慮しないツールであれば「染色体名が異なる」とエラーになり、これが地味に面倒です。
ゲノムの染色体名の方をEnsemblに揃えるやり方もありますが、ここではEnsemblの遺伝子名にchrを追加する方法を試してみましょう。

修正前

$ head Homo_sapiens.GRCh38.94.gtf
#!genome-build GRCh38.p12
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.27
#!genebuild-last-updated 2018-07
1       havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
1       havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
1       havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";

冒頭が染色体名になっているので、行頭にchrを付加するだけでよいですね。これは以下のようなワンライナーで実行可能です。

# 冒頭のコメントアウト文を除外し、各行の行頭にchrを追加して出力
$ grep -v \# Homo_sapiens.GRCh38.94.gtf | awk '{print "chr" $0 }' > Homo_sapiens.GRCh38.94.addchr.gtf
# 内容確認
$ head Homo_sapiens.GRCh38.94.addchr.gtf
 head Homo_sapiens.GRCh38.94.addchr.gtf
chr1    havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
chr1    havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
chr1    havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    havana  exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    havana  exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    havana  transcript      12010   13670   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12010   12057   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12179   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12613   12697   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12975   13052   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";

無事追加されました。 冒頭のコメント文を除外したくない場合は後から付け足すか、最初のコマンドでgrep部分を使わずに出力してください(後で手作業で#前のchrを取り除く)。

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

前回の記事の続きです。

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

発現量データ生成

前回の記事で発現量データを生成するのを忘れていたので、ここで生成します。stringtieを使います。

for prefix in ERR188044 ERR188104 ERR188234 ERR188245 ERR188257 ERR188273 ERR188337 ERR188383 ERR188401 ERR188428 ERR188454 ERR204916
do
    stringtie -e -B -p 16 -G stringtie_merged.gtf \
    -o ballgown/${prefix}/${prefix}_chrX.gtf \
    ${prefix}_chrX.bam
done

これでballgown/${prefix} というディレクトリが自動生成され、その中に発現量データが出力されます。
ディレクトリの中身を見てみましょう。

$ ls ballgown/ERR188044/
ERR188044_chrX.gtf  e2t.ctab  e_data.ctab  i2t.ctab  i_data.ctab  t_data.ctab

gtfファイルと、いくつかのテキストファイルができています。 eはexon, iはid, tがtranscriptの情報を表しているようですね。

比較するサンプルの確認

サンプルの情報はchrX_data/geuvadis_phenodata.csvに記述されています。 どんなサンプルなのか確認しましょう。

$ cat chrX_data/geuvadis_phenodata.csv
"ids","sex","population"
"ERR188044","male","YRI"
"ERR188104","male","YRI"
"ERR188234","female","YRI"
"ERR188245","female","GBR"
"ERR188257","male","GBR"
"ERR188273","female","YRI"
"ERR188337","female","GBR"
"ERR188383","male","GBR"
"ERR188401","male","GBR"
"ERR188428","female","GBR"
"ERR188454","male","YRI"
"ERR204916","female","YRI"

サンプルIDごとに、性別とpopulation情報が載っています。 解析内容は、population間のばらつきを除去しながら性別間の遺伝子発現変動を見る、というもののようです。

Ballgownで解析

BallgownはRで動きますので、同じディレクトリでRを起動しましょう。
以降のコマンドはR上で実行するものとします。Nature Protocolのコマンドから少し修正しています。

# ライブラリ読み込み
library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
# サンプル情報の読み込み(通常は自分で用意する)
pheno_data = read.csv("chrX_data/geuvadis_phenodata.csv")
# 先程生成した発現量データの読み込み
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)

# bg_chrXに格納されたものを確認
bg_chrX
ballgown instance with 3520 transcripts and 12 samples

bg_chrX という名前のballgownインスタンスが生成されました。

# 低発現遺伝子の除去
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)

# 発現変動遺伝子の推定(transcript単位)
results_transcripts = stattest(bg_chrX_filt,
                   feature="transcript",
                   covariate="sex",
                   adjustvars = c("population"),
                   getFC=TRUE,
                   meas="FPKM")

# 得られた発現変動遺伝子の確認
head(results_transcripts)
     feature id        fc      pval      qval
1 transcript  1 0.9347423 0.7043147 0.9373013
2 transcript  2 1.2006990 0.8716451 0.9723579
3 transcript  3 1.0173174 0.9896415 0.9990549
4 transcript  4 0.3814504 0.5200179 0.9192672
5 transcript  5 0.6014996 0.3216161 0.9192672
6 transcript  6 0.6540123 0.3264784 0.9192672

同様に、遺伝子単位の発現変動情報を抽出します。

results_genes = stattest(bg_chrX_filt,
             feature="gene",
             covariate="sex",
             adjustvars = c("population"),
             getFC=TRUE,
             meas="FPKM")

head(results_genes)
  feature         id        fc      pval      qval
1    gene    MSTRG.1 1.0824443 0.6847873 0.9196130
2    gene   MSTRG.10 0.9335310 0.8132130 0.9599962
3    gene  MSTRG.100 0.6126762 0.5687947 0.8823732
4    gene MSTRG.1000 1.2044223 0.4500794 0.8360675
5    gene MSTRG.1001 1.3799856 0.6140598 0.8936018
6    gene MSTRG.1003 0.9087684 0.5909100 0.8867180

以下ではresults_transcripts に遺伝子名の情報を付加します。

results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
                 geneIDs=ballgown::geneIDs(bg_chrX_filt),
                 results_transcripts)

head(results_transcripts)
  geneNames geneIDs geneNames.1 geneIDs.1    feature   id          fc
1         . MSTRG.2           . MSTRG.528 transcript 1668 0.030602170
2    PLCXD1 MSTRG.2        XIST MSTRG.528 transcript 1667 0.002990002
3         . MSTRG.2           . MSTRG.528 transcript 1666 0.015949119
4         . MSTRG.2           . MSTRG.528 transcript 1669 0.028359333
5         . MSTRG.3        TSIX MSTRG.527 transcript 1665 0.077438632
6    PLCXD1 MSTRG.2           . MSTRG.612 transcript 1862 7.335431029
          pval         qval
1 1.474474e-10 2.824911e-07
2 2.497711e-10 2.824911e-07
3 6.564345e-10 4.949516e-07
4 6.871521e-08 3.885845e-05
5 1.911435e-06 8.647331e-04
6 1.331176e-05 5.018534e-03

遺伝子名はid順にソートされていますが、発現変動遺伝子が見たいので、p値の小さい順にソートしなおします。

results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)

head(results_transcripts)
  geneNames   geneIDs    feature   id          fc         pval         qval
1         . MSTRG.528 transcript 1668 0.030602170 1.474474e-10 2.824911e-07
2      XIST MSTRG.528 transcript 1667 0.002990002 2.497711e-10 2.824911e-07
3         . MSTRG.528 transcript 1666 0.015949119 6.564345e-10 4.949516e-07
4         . MSTRG.528 transcript 1669 0.028359333 6.871521e-08 3.885845e-05
5      TSIX MSTRG.527 transcript 1665 0.077438632 1.911435e-06 8.647331e-04
6         . MSTRG.612 transcript 1862 7.335431029 1.331176e-05 5.018534e-03

head(results_genes)
  feature        id          fc         pval         qval
1    gene MSTRG.528 0.002409209 2.041411e-11 2.086323e-08
2    gene MSTRG.527 0.082689539 6.110906e-06 2.409310e-03
3    gene MSTRG.143 3.158100078 7.719104e-06 2.409310e-03
4    gene MSTRG.612 7.327271932 1.073484e-05 2.409310e-03
5    gene MSTRG.763 0.278924648 1.178723e-05 2.409310e-03
6    gene MSTRG.618 9.109365087 4.555636e-05 7.759766e-03

以下のコマンドを実行すると、得られたresults_transcripts, results_genesをテキストファイルに出力します。(ファイル名は任意)

write.csv(results_transcripts, "chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv", row.names=FALSE)

可視化

以下は可視化のコマンドです。得られた発現変動遺伝子についてさらに調査する時に用います。

# 色の設定
tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)

# 箱ひげ図
fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

ラベルが見切れていますが、サンプルごとの遺伝子発現量の分布が性別ごとに色分けされて箱ひげ図で表示されました。 縦軸のFPKMというのは遺伝子ごとにマップされたリード数を遺伝子長と全リード数で正規化したものです。

次に、特定の遺伝子について可視化してみましょう。 ここでは例として12番めの遺伝子について可視化します。

ballgown::transcriptNames(bg_chrX)[12]
         12 "NR_027231" 
ballgown::geneNames(bg_chrX)[12]
         12 "LINC00685" 

おや、、、何故だかNature Protocolの例と表示される遺伝子名が違います。順番が変わったのでしょうか。
ともかく表示してみましょう。

plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
     main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
        ballgown::transcriptNames(bg_chrX)[12]),
     pch=19,
     xlab="Sex",
     ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
       col=as.numeric(pheno_data$sex))

サンプルごとの発現量のばらつきを性別で分けてプロットしているのですが、全てのサンプルで発現量が0のため、差が見えません。
この遺伝子ではない遺伝子が例になっていたはずなのですが。僕のスクリプトにどこか間違っているところを見つけたら教えてください。

plotTranscripts(ballgown::geneIDs(bg_chrX)[1729],
        bg_chrX,
        main=c('Gene XIST in sample ERR188234'),
        sample=c('ERR188234'))

遺伝子の簡単なisoform情報も表示することができます。

plotMeans('MSTRG.56', bg_chrX_filt, groupvar="sex", legend=FALSE)

ここではおそらく性別間で差がある遺伝子構造を可視化しているのですが、差がない遺伝子を指定していまっているのでよくわかりません。

おわりに

最後ミスったためしまりがないことになってしまいましたが、ともかく一通りの作業ができるようになりました。
新規転写物を探索したい時には有効なワークフローなので、機会があれば利用されると良いと思います。