表紙(FrontPage) | 編集(管理者用) | 差分 | 新規作成 | 一覧 | RSS | 検索 | 更新履歴

RNA-seqデータの分析について勉強する -

目次

RNA-seq は当たり前に行われるようになった。様々なデータ解析用のソフトウェア、わかりやすい日本語の解説も多数公開されている。

そこで、公開されているデータを、自分の研究に役立てることを目標とする。 https://www.ncbi.nlm.nih.gov/geo/  Gene Expression Omnibus から多くのデータが得られる。  https://www.omicsdi.org/ は、RNA-seq などのデータを検索しやすい形にまとめてあり、興味深いデータを見つけやすい。 DDBJ (DNA Data Bank of Japan) もデータを収集し検索できるようにしている。  https://www.ddbj.nig.ac.jp/dra/index-e.html  European Nucleotide Archive (ENA)    https://www.ebi.ac.uk/ena/browser/home   も使いやすい。最近は ENA からデータを入手することが多い。 

やってみると様々な問題が出てくる。それらをどう解決したかを記載してみる。 また、ソフトウェアの世界はつねにものすごい速度で進歩してアップデートされている。ここに書いたこともあっという間に古くなっている。何かをアップデートしたら記録してみる。

RNAseq データを処理するソフトウェアは何種類も開発されている。一種類だけ使っていると、そのソフトウェアの性質、癖によって知らないうちに結果の解釈に偏りが出たりする可能性もある。またソフトウェアに与えるパラメーターが自分の用途には不適切な値なことに気がつかずに、ずっと使ってしまうこともあり得る。複数のソフトウェアを使えるようにすればそういう問題が起きにくくなるだろう。だから同じ働きをするソフトウェアでも何種類も使えるようにすることも目標にする。


https://bioinformatics.uconn.edu/rnaseq-arabidopsis   ではモデル植物であるシロイヌナズナの RNA-seq データ分析を親切に解説している。 これをそのまま真似てみる。

見本として取り上げられているデータは https://www.ncbi.nlm.nih.gov/sra?term=SRX1756762 に概要が書かれている。 Illumina HiSeq 2500 で分析している。


GEO database を見てみると、イルミナ社のシーケンサーを使用している例が多い。この場合、サンプル調製には「スタンダードmRNA キット」「Total RNA キット」「Small RNA キット」の3種類がある。たいていのものはスタンダードキットを用いていて、核ゲノムから転写されポリAが付いている標準的な mRNA のみが分析される。 その場合、オルガネラゲノムから転写される mRNA などは分析しにくい。しかしオルガネラゲノムから転写される mRNA の分析をしている人もいる。   http://bfg.oxfordjournals.org/content/12/5/454   RNA-Seq data: a goldmine for organelle research    David Roy Smith   Briefings in Functional Genomics Volume 12, Issue 5Pp. 454-456   オルガネラゲノムからの転写産物は AT-rich なことが多いので、ポリA 選択で混入することがあるらしい。 またオリゴ dT をプライマーにする場合、ポリ A がない RNA でも A リッチな部分にプライマーがついて DNA 合成が起きる可能性がある。 そういう現象を利用するので、オルガネラゲノムからの転写産物にたいしては定量的なことはわからないだろう。定性的に、こういう転写産物があるということを知ることしかできない。それでも、オルガネラゲノムの配列が知られていない生物種の場合はとても役に立つだろう。またモデル生物やヒトでもオルガネラゲノムに生じている変異を検出するには十分使えるのではないか。

上に書いた Smith DR 博士らによって、オルガネラ転写産物を分析するためのソフトウェア Chloroseq が開発された。   http://github.com/BenoitCastandet/chloroseq   論文としても発表されている。   https://www.ncbi.nlm.nih.gov/pubmed?linkname=pubmed_pubmed&from_uid=27402360    ChloroSeq, an Optimized Chloroplast RNA-Seq Bioinformatic Pipeline, Reveals Remodeling of the Organellar Transcriptome Under Heat Stress.   Castandet B, Hotto AM, Strickler SR, Stern DB.   G3 (Bethesda). 2016 Sep 8;6(9):2817-27. doi: 10.1534/g3.116.030783.   PMID: 27402360

A Guide to the Chloroplast Transcriptome Analysis Using RNA-Seq.   Michel EJS, Hotto AM, Strickler SR, Stern DB, Castandet B.   Methods Mol Biol. 2018;1829:295-313. doi: 10.1007/978-1-4939-8654-5_20.   PMID: 29987730

現在では様々な non-coding RNA の機能が解明されたことによりポリA 鎖がついていない RNA の解析の重要性が増している。そのためポリA 鎖をもつ RNA 以外も標的にしたデータも増えている。

High-throughput m6A-seq reveals RNA m6A methylation patterns in the chloroplast and mitochondria transcriptomes of Arabidopsis thaliana.   Wang Z, Tang K, Zhang D, Wan Y, Wen Y, Lu Q, Wang L.   PLoS One. 2017 Nov 13;12(11):e0185612. doi: 10.1371/journal.pone.0185612. eCollection 2017.   PMID: 29131848 という論文では、GSE72706 という公開データを分析している。

ArrayExpress では、「Type」の部分をクリックするとタイプで並び替えられて、「RNA-seq of non coding RNA」のタイプのデータを集めることができる。miRNA を分析しているデータが多い。


https://bioinformatics.uconn.edu/rnaseq-arabidopsis   ではモデル植物であるシロイヌナズナの RNA-seq データ分析を親切に解説している。これを真似てみる。データを自分のマシンにダウンロードするには、SRA Toolkit という専用のソフトウェアを使う。

SRA とはどんなものか。http://www.ncbi.nlm.nih.gov/books/NBK47540/   が説明である。Sequence Read Archive を略して SRA という。 http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software   SRA Toolkit などのソフトウェアが公開されている。これらのものを開発して頂いているおかげで、データを誰でも操作できるようになっている。

前処理、分析のためのソフトウェア

SRA Toolkit (WSL2 + Ubuntu)

Ubuntu 20.04 を入れてあるマシンに SRA Toolkit をインストールしてみる。

WSL2 という便利な仕組みが使えるようになり、WINDOWS10 のマシンに WSL2 を介して Ubuntu をインストールして使用できるようになった。 WINDOWS と組み合わせて使うことができて使いやすい。

BIOCONDA https://bioconda.github.io/ という便利なソフトウェアセンターが公開されている。 ソフトウェアを一つずつばらばらに入れていると管理しにくい。スマートフォンでアプリを入れるのと同じように、どんなソフトウェアもソフトウェアセンターから入れる方が管理しやすい。 BIOCONDA は、Miniconda, Anaconda という管理システムを使って、生物学関連のソフトウェアを集めて簡単にインストールできるようになっている。 ソフトウェアをアップデートするときは、「conda update sra-tools」のようにコマンドを打って行う。


しかし、RNA-seq データ分析ソフトウェアの内には conda でインストールすると他のソフトウェアに影響を与えてうまく動かなくなるものがある。Python という言語で書かれている古いソフトウェアはバージョン 2.7 を使うものがある。バージョン 3 の Python と conflict を起こすことがある。このことについては、この文章よりもっと詳細でわかりやすい解説が公開されている。   http://imamachi-n.hatenablog.com/entry/2017/01/14/212719 「biocondaを利用してNGS関連のソフトウェアを一括でインストールする」 Imamachi-n 氏による解説にあるように、Python のバージョンが異なる仮想環境を作ることができる。私は解説を見て以下のようにした。

 $ conda create --name py27 python=2.7
 py27に入る     (base) $ conda activate py27
 py27から抜ける (py27) $ conda deactivate

Python2.7 を用いる古いソフトウェアは [py27] がプロンプトに出た状態で「conda install .....」とする。py27 を activate した状態では Python2.7 を使う環境で動くようになり、他のソフトウェアに影響を与えなくなる。

ソフトウェアは新しいものが次々と開発され、バージョンアップもされるので今では Python 2.7 の必要は小さくなった。しかし「この古いソフトウェアは開発者が忙しくなって Python3 に対応される見込みがないが代わりがないので、どうしても使いたい」ということもありうるので、そういう事態に対処できるようにしておく方がよい。例: The Molecular Modeling Toolkit   http://dirac.cnrs-orleans.fr/MMTK.html


sickle-trim という RNA-seq データの品質をチェックするソフトウェアがある。これと全然別の機能を持つ sickle というソフトウェアも bioconda に入っていて間違えやすい。こういうこともあるので注意しないといけない。bioconda のページを見るとそれぞれのソフトウェアについて解説されているので、それらをよく見ないといけない。


SRA Toolkit は、コマンドラインから使うプログラムの詰め合わせである。これも BIOCONDA に入っていた。

http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=toolkit_doc&f=std#s-3   SRA Toolkit Installation and Configuration guide にあるように、テストしてみる。

 fastq-dump --split-files -X 5 -Z SRR5330620

と打つと、データの5行分が表示された。このように一部を見るだけなら待ち時間は少ない。この方法で出てくるデータは、fastq というフォーマットで記述されている。細かい型式はシーケンサーが異なる会社の機種だと変わる。

DDBJ Sequence Read Archive (DRA)、European Nucleotide Archive (ENA)

fastq-dump は NCBI (SRA) からデータを fetch するが速度が出ないことがある。DDBJ (DNA Data Bank of Japan) というものもある。   この場合は  https://www.ddbj.nig.ac.jp/dra/index-e.html   をブラウザで開き、Search -> Accession number で検索する。 Accession number は NCBI GEO database で調べられる。検索するとファイルの情報は出てくるが、SRR で始まる Accession number の fastq ファイルは存在しなかった。DRRで始まるデータは入っている。

European Nucleotide Archive   https://www.ebi.ac.uk/ena/browser/home  からもデータを取得することができる。この場合は Edge などの web browser を使ってデータを選びダウンロードでき、使いやすい。

fastq について

一つの read は、4行で表現される。一行目は、@ から始まり、説明文になっている。次の行は塩基配列である。3行目は + で始まり、1行目と同じことが書いてあることもあれば、+ しか書かれていないこともある。 4行目は、データ(塩基)の信頼性を表す情報が入っている。

fastq-dump のオプションについて

fastq-dump --help で見ると、たくさんのオプションがある。 例:

 fastq-dump --split-files --gzip SRR5330620

最後にアクセッション番号を書いて、取り出したいデータを指定する。

データの種類が Layout: PAIRED の場合、オプションに --split-files を付けると AAAAA_1, AAAAA_2 のようにペアとしてデータを取得できる。

データが SINGLE の場合でも、このオプションをつけたままで問題ないので、いつもつけておくのがよい。

いろいろなファイル

(multi-) fasta 形式ファイルは、大量の塩基配列とそれらの簡単な注釈(一行ずつ)を含むデータで、様々な分析に供せられる。

fastq 形式ファイルは、大量の塩基配列とそれらの信頼性とそれらの簡単な注釈(一行ずつ)を含むデータで、様々な分析に供せられる。

SRAファイルのアクセッション番号が分かれば、SRA Toolkit の fastq-dump を用いて、直接 fastq ファイルを取り込める。

fai ファイルは、fasta ファイルから生成したインデックスを含む。染色体の名前、染色体の長さ、そのデータのスタート位置(ファイルの先頭から数える)、一行の文字数、一行のバイト数(文字+改行コードで数える) が、データごとに書かれている。

SAM 型式ファイルは、HISAT2 のようなマッピングソフトウェアの出力である。それをさらに BAM 形式に変換して、カウント情報の取得などを行うことがある。SAMtools http://samtools.sourceforge.net/ というソフトウェアで変換できる。 featureCounts のように SAM 形式のままで処理してくれるソフトウェアもある。

SAM ファイルはテキストデータになっていて、人間が一応読むことができる。BAM ファイルは SAM ファイルと同じ情報を持っているが、ソフトウェアで処理しやすいバイナリファイルになっていて、人間が見てもわからない。SAMtools は、BAM ファイルが持っている情報の必要な部分だけを出力するためにも用いられる。その出力をエクセルや他のプログラムで処理することもできる。

BED 形式ファイルは、ゲノム上の特定の領域を示す情報が書かれているファイルで、転写因子の結合部位などを示すのに使われる。BAM ファイルと対になっていることが多い。ChIP 実験の結果のファイルには、BAM ファイルと BED ファイルが同梱されていることがある

GSM861508_PM1_m1_btb_chrom.bed というファイルでは、8601636 行のデータが納められていた。目で眺めても全く意味をつかむことはできない。BED 形式ファイルは、他のプログラムの入力ファイルとして使われる。 他のプログラムによって、人間が見てわかるような形式で表示する。

GFF/GTF形式ファイル は、遺伝子アノテーションの情報を含む。http://ccb.jhu.edu/software/tophat/index.shtml から、「Index and annotation downloads」へ進んで、それぞれの生物用のファイルを入手することができる。

GFF/GTF形式ファイルには、GTF2 フォーマットと GFF3 フォーマットが存在する。GTF2 でないとうまく動かないソフトウェアが存在する。GFF3 を GTF2 へ gffread というソフトウェアでコンバートできた。参考にしたページ: http://ccb.jhu.edu/software/stringtie/gff.shtml    gffread は Bioconda からインストールできる。 > conda install gffread

refFlat形式ファイルは、遺伝子アノテーションの情報を含む。

データの品質チェック

https://bioinformatics.uconn.edu/rnaseq-arabidopsis   で解説されているように、まずデータの品質チェックを行う。

データを fastqc に与えて品質を可視化する。結果が書かれた html ファイルが作られるので、それをブラウザで開いて見ることができる。

 fastqc /ファイルの場所/trimmed_SRR3498212.fastq   

SRR3498212 をチェックすると、× マークがついている評価点が三つもあった。Per base sequence content, Sequence duplication levels, Adapter content に、× が付いている。それぞれの配列の 30bp よりも後は、ほとんどがアダプターの配列になってしまっている質の悪いデータだった。実際にこのデータは hisat2 で処理してもうまくマッピングされる配列がほとんどなかった。データの質をチェックすることが大切だと言うことがよくわかった。trimmomatic というソフトウェアを使うとアダプター配列を除去できる。これも Bioconda からインストールできる。  http://www.usadellab.org/cms/?page=trimmomatic がオリジナルの場所なので、使い方などはここを見ないといけない。

SRR3229130 を sickle で処理したもので行うと、× マークは一つもなかった。アダプターの配列は全く含まれていない。実際にこのデータは hisat2 で 99.47 % が align される質のよいデータだった。

マッピング

HISAT2 などのソフトウェアで、RNAseq のデータをゲノム配列にマッピングできる。 マッピングプログラムには、二種類の入力ファイルを与える。一つは fastq データで、もう一つはマッピングされる側のゲノム配列である。

マッピングされる側のデータは、どういうものか。ゲノム全体配列が入っているデータを元にして、それをインデックスファイルにしたものを使う。 インデックスファイルは、マッピングプログラムごとに別のものを用意する。

HISAT2

hisat2-build というコマンドを使ってインデックスファイルを用意する。

 hisat2-build -f /ファイルの場所/ファイル名.fa /出力するディレクトリー/インデックス名

このとき一つ問題があった。私が入手したゲノム全体の配列ファイル Arabidopsis.thaliana.TAIR10.dna.chromosome.1.fa というファイルでは、染色体の名前が 1, 2, 3, 4, 5, Mt, Pt と記入されていた。しかし後のステップで使用する別の種類のファイルである Athaliana_167_TAIR10.gene.gff3 や TAIR10_GFF3_genes.gff は、それぞれの染色体の名前が Chr1, Chr2, Chr3, Chr4, Chr5, ChrM, ChrC と記入されていた。このように染色体の名前が異なっているとソフトウェアがうまく働かなくなる。 この問題に対処するためにゲノム配列のファイルをエディタに読み込んで、染色体の名前を変更(1 -> Chr1, 2 -> Chr2, ... 行のはじめから >1 や >2 のように書かれていたので、>Chr1 のように直した)して他のファイルに合わせたものをセーブして hisat2-build によるインデックスの作成に使用した。

インデックスができたらマッピングをする。Manual に、次のように書かれている。

 hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]

オプションはもっと多数存在する。

同じ Illumina のマシンによるデータでも、リードの長さが違ったり、ペアでなかったりするのでデータに関する説明をよく読んで適切なオプションをつける必要がある。

fastQC で良い品質であることが確かめられた SRR3229130 を、以下のように処理してみた。

 $ hisat2 -p 4 -x /場所/index/TAIR10 --dta -q -U /場所/trimmed_SRR3229130.fastq -S SRR3229130.sam

すぐに結果が出る。

 SRR3229130.sam
 8808444 reads; of these:
   8808444 (100.00%) were unpaired; of these:
     46532 (0.53%) aligned 0 times
     8564820 (97.23%) aligned exactly 1 time
     197092 (2.24%) aligned >1 times
 99.47% overall alignment rate

質のよいデータだと、ほとんどが align された。

ペアエンドのデータの場合の例:

 $ hisat2 -p 4 --dta -x /TAIR10/index/TAIR10 -q -1 SRR5905046_1.fastq  -2 SRR5905046_2.fastq -S SRR5905046.sam

 22636976 reads; of these:
   22636976 (100.00%) were paired; of these:
     1925817 (8.51%) aligned concordantly 0 times
     20410609 (90.16%) aligned concordantly exactly 1 time
     300550 (1.33%) aligned concordantly >1 times
     ----
     1925817 pairs aligned concordantly 0 times; of these:
       267792 (13.91%) aligned discordantly 1 time
     ----
     1658025 pairs aligned 0 times concordantly or discordantly; of these:
       3316050 mates make up the pairs; of these:
         2359882 (71.17%) aligned 0 times
         938919 (28.31%) aligned exactly 1 time
         17249 (0.52%) aligned >1 times
 94.79% overall alignment rate

samtools の使い方

マッピングソフトウェアの出力である sam ファイルを bam ファイルに変換する、ファイルから人間が読める形で情報を取り出すには、samtools を使う。

HISAT2 で出力された SRR3229130.sam を、「ソートされた sorted BAM files」に変換する。このステップは、Stringtie という発現量カウントのためのソフトウェアで処理するために必要になる。bam ファイルは圧縮されてサイズが小さくなり、他のソフトウェアで処理しやすい。

 samtools sort -@ 4 -o SRR3229130.bam SRR3229130.sam

Stringtie で Transcript assembly and quantification

読み出された配列を遺伝子ごとに集計し、カウントすることでそれぞれの遺伝子の発現量を測定することができる。そのために、遺伝子領域の情報を含むファイルである gff3 または gtf ファイルが必要になる。

Athaliana_167_TAIR10.gene.gff3 というファイルを検索して、https://github.com/k821209/BAMVIS-GENE というサイトから download した。 サイズは 38.4 MB (38412591 バイト)、このような内容から始まっていた。

 ##gff-version 3
 ##annot-version TAIR10
 Chr1	phytozomev10	gene	3631	5899	.	+	.	 ID=AT1G01010.TAIR10;Name=AT1G01010
 Chr1	phytozomev10	mRNA	3631	5899	.	+	.	ID=AT1G01010.1.TAIR10;Name=AT1G01010.1;pacid=19656964;longest=1;  Parent=AT1G01010.TAIR10
 Chr1	phytozomev10	five_prime_UTR	3631	3759	.	+	.	ID=AT1G01010.1.TAIR10.five_prime_UTR.1;Parent=AT1G01010.1.TAIR10;pacid=19656964
 Chr1	phytozomev10	CDS	3760	3913	.	+	0	ID=AT1G01010.1.TAIR10.CDS.1;Parent=AT1G01010.1.TAIR10;pacid=19656964

https://www.arabidopsis.org/download/index-auto.jsp?dir=%2Fdownload_files%2FGenes%2FTAIR10_genome_release%2FTAIR10_gff3 から、TAIR10_GFF3_genes.gff というファイルも入手した。 Athaliana_167_TAIR10.gene.gff3, TAIR10_GFF3_genes.gff のどちらでも結果は同じだった。


https://www.arabidopsis.org/download/index-auto.jsp?dir=%2Fdownload_files%2FGenes%2FAraport11_genome_release には、 Araport11_GFF3_genes_transposons.201606.gff.gz 17,839 KB 2019-07-11  が公開されている。


次に stringtie を実行した。https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual にマニュアルがある。よく読まないといけない。

 $ stringtie SRR3229130.bam -p 4 -G /場所/Athaliana_167_TAIR10.gene.gff3 -o SRR3229130.gtf 

入手した gff, gff3 ファイルは、それぞれの染色体の名前が Chr1, Chr2, Chr3, Chr4, Chr5, ChrM, ChrC と記入されていた。しかし、私が入手したゲノム全体の配列ファイル Arabidopsis.thaliana.TAIR10.dna.chromosome.1.fa というファイルでは、それぞれが 1, 2, 3, 4, 5, Mt, Pt と記入されていた。このように染色体の名前が異なっているとStringtie は Gene ID をうまく取り出せなくなる。「Please make sure the -G annotation file uses the same naming convention for the genome sequences. 」という warning message が出る。

この問題に対処するためにゲノム配列のファイルをエディタに読み込んで、染色体の名前を変更(1 -> Chr1, 2 -> Chr2, ...)して他のファイルに合わせたものをセーブして hisat2-build によるインデックスの作成に使用した。 そうすると warning は出なくなった。この問題の解決に参考になったページ: 

https://wiki.cyverse.org/wiki/display/DEapps/Evolinc+in+the+Discovery+Environment

https://github.com/griffithlab/rnaseq_tutorial/wiki/Annotation#important-notes

https://github.com/igvteam/igv.js/issues/507

オプションに -e をつけるとマイクロアレイのデータ(遺伝子ごとに発現量が計算される)と似たデータを出力できる。

RNA-seq データはほとんどの場合複数のデータがセットになっている。それぞれについて上のように処理して gtf ファイルができる。次のステップの処理のためにそれらの gtf ファイルを結合 merge する必要がある。

まず、mergelist.txt というファイルを作る。 これは 「ls *.gtf > mergelist.txt」のようにすると作成できる。次に stringtie を 今度は --merge オプションをつけて実行する。すぐに終了する。

 $ stringtie --merge -p 4 -G /場所/Athaliana_167_TAIR10.gene.gff3 -o merged_data.gtf mergelist.txt

次に、ballgown というカウントソフトウェアの入力ファイルを、マージした gtf ファイルから生成するために、stringtie をオプションを変えて (-B) 実行する。

 
まず、ballgown のためのデータをまとめておくためのディレクトリ ballgown を作っておく。

 $ stringtie サンプル1.bam -e -B -p 4 -G merged_data.gtf -o /場所/ballgown/サンプル1/ball_サンプル1.gtf

ballgown に読ませるための gtf ファイルが生成する。同時にいくつかの ctab ファイルが生成する。これらはセットになっているのでサンプルごとに同じディレクトリにまとまっていないといけない。 bam ファイルの名前、出力ファイルの名前を変えて、サンプルごとに実行する。

R で ballgown を実行

ここからは R 言語のパッケージである ballgown というソフトウェアを使う。 私は Rstudio をインストールして、そこから使っている。


ソフトウェアはつねにアップデートされている。2020年5月に ballgown を使おうとしたら、「biocManager という package のインストールが必要」というメッセージが出た。そこで Rstudio から biocManager パッケージをインストールした。そして

 biocManager::install("ballgown")

で ballgown をインストールし直すことができた。その際には必要になる様々なソフトウェアが自動でインストールされる。

以下に書いてあることも、つねに少しずつ変化している


https://bioinformatics.uconn.edu/rnaseq-arabidopsis   に書かれているように必要なパッケージをインストールする。また、http://rnakato.hatenablog.jp/entry/2018/11/26/145847   Ryuichiro Nakato 氏による解説もとてもわかりやすく参考になった。以下のようにスクリプトを作って実行する。

 install.packages("devtools")
 devtools::install_github("vqv/ggbiplot")
 source("http://www.bioconductor.org/biocLite.R")
 biocLite(c("alyssafrazee/RSkittleBrewer","ballgown","genefilter","dplyr","devtools"))

このときうまくインストールできないことがあった。メッセージをよく見ると「libcurl4-openssl-dev」というソフトウェアをインストールすることが必要だった。これは R とは関係なく、

 sudo apt install libcurl4-openssl-dev 

とすることでインストールできた。そうした後ではうまくいった。

https://bioinformatics.uconn.edu/rnaseq-arabidopsis   に書かれているように必要なライブラリをロードする。以下のようにスクリプトを作って実行する。

 library(ballgown)
 library(RSkittleBrewer)
 library(genefilter)
 library(dplyr)
 library(ggplot2)
 library(gplots)
 library(devtools)
 dir <- "ballgown ディレクトリがある場所" # ここの内部に ballgown ディレクトリがあり、その内部にサンプルごとのディレクトリがある
 setwd(dir)

まず、ballgown のために用意したサンプルの名前などを設定したテキストファイルを用意する。phenodata.csv という名前のファイルにする。上に書いた dir で設定したディレクトリに置く。http://rnakato.hatenablog.jp/entry/2018/11/26/145847   Ryuichiro Nakato 氏による解説に説明されている。例:

 "ids", "part"
 "SRR2932182", "Col_root"
 "SRR2932183", "phb3_root"

ids の部分はサンプルごとのディレクトリにつけた名前で、注釈を "part" などの名前でつけられる。注釈はもっと多くできる。"part" の部分が同じデータは、繰り返し実験した一群のデータと解釈され、平均や分散が計算される。

Rstudio のコンソールから、以下のように入力して設定ファイルを読み込む。

 pheno_data <- read.csv("phenodata.csv")

コンソールで > pheno_data と打って、内容を確認する。

次に

 bg <- ballgown(dataDir = "ballgown", pData=pheno_data, samplePattern = "SRR")

と打った。データは ballgown というディレクトリに置いた。設定ファイルは pheno_data に読み込んだ。私は ballgown ディレクトリの内側にデータとして SRR2932182, SRR2932183 という二つのディレクトリを置いたので、それらの頭に共通する SRR をパターンとした。

 Reading linking tables
 Reading intron data files
 Merging intron data
 Reading exon data files
 Merging exon data
 Reading transcript data files
 Merging transcript data
 Wrapping up the results

とメッセージが出て終了する。その結果が、ballgown データを保持するオブジェクト(この場合、オブジェクトは結果となるデータを収納するコンテナのようなものとして具体化する)である bg に収納される。データには「整数型」「文字列型」など様々な種類があるが、この場合オブジェクト bg は ballgown 型のデータを保持することになる。このことを「オブジェクト bg は、ballgown 型のインスタンス(実体)である」と表現する。インスタンスという用語は、必ずそのデータ型と組み合わさって使われる。

オブジェクト bg に保持された、ballgown の結果に補正を加え発現変化などを検出する。

texpr(bg) と打つと、結果を保持するオブジェクト bg の FPKM 値が出力される。

       FPKM.SRR2932182 FPKM.SRR2932183
 1           15.012596       13.118468
 2            0.000000        0.000000
 3            0.791032        0.404203
 4            9.948609        6.854258
 5            0.535599        0.748147
 6            3.567534        7.780320
 7            0.000000        0.000000
 8           11.617834        0.230532

texpr(bg, 'all') と打つと、結果を保持するオブジェクト bg のすべての内容が出力される。これらの結果に ID などを付加してデータフレームとして出力することもできる。

 results_transcripts <- texpr(bg, 'all')         または results_transcripts = data.frame(geneNames = ballgown::geneNames(bg), geneIDs = ballgown::geneIDs((bg)), texpr(bg))

データフレームをファイルとして保存する。

 write.table(results_transcripts, "result_GSE96041.txt", quote=F, sep="\t")

これによって生成したファイルをエクセルなどで読み込んで、アノテーションをつけるなどの処理をすることができる。

結果を保持するオブジェクト bg をさらに処理することもできる。

出現回数が小さい、ほとんど発現していない遺伝子のデータは排除する:

 bg_filt = subset(bg, "rowVars(texpr(bg))>1" ,genomesubset=TRUE)

 > texpr(bg_filt)
       FPKM.SRR2932182 FPKM.SRR2932183
 1           15.012596       13.118468
 4            9.948609        6.854258
 6            3.567534        7.780320
 8           11.617834        0.230532
 9           67.906143       85.427742
 11           1.377149        5.713247

stattest で、統計的に発現変化が検出される遺伝子を抽出する:これを実行するには、phenodata.csv の "part" で指定されるデータ群のメンバーが2個以上ないといけない(そうでないとデータ群に対する分散を計算できない)

 results_transcripts = stattest(bg_filt, feature="transcript" , covariate = "part" , getFC = TRUE, meas = "FPKM")

無理に R で処理しなくても、エクスポートした結果をエクセルや他の言語で書いたソフトウェアに読み込んで処理したほうがやりやすいかもしれない。

Ballgown was not really designed for *gene*-level differential expression analysis − it was written specifically to do *isoform*-level DE.

私の場合、RNAseq のデータをマイクロアレイのデータと同じように扱いたい。しかし Ballgown はそのようには設計されていないと書かれていた。  https://support.bioconductor.org/p/107011/#110717  「DESeq2 vs Ballgown results」 バイオインフォマティクスのソフトウェアは、特定の目的を念頭に置いて作られていることが多く、しかもそれがはっきり書かれていなかったりするらしい。マニュアルをよく読み、そのソフトウェアに関する論文があれば目を通して、気をつけないといけない。

「 Using DESeq2 with FeatureCounts is a much better-supported operation if your main interests are in gene-level DE.」 と書かれている。

featureCounts (subread) で sam または bam ファイルの結果を遺伝子ごとのカウント数にする

RNAseq データを処理するソフトウェアは何種類も開発されている。一種類だけ使っていると、そのソフトウェアの性質、癖によって知らないうちに結果の解釈に偏りが出たりする可能性もある。 東大の門田先生が公開されている資料に、ソフトウェアのマニュアルを丁寧に読み解くことで、そのソフトウェアが必要とするデータの形式や、特別な条件が生じた際の挙動を正しく理解する具体例が紹介されていた。   農学生命情報科学特論I 第1回 2018年の資料   

例として「HTSeqは、複数個所にマップされるリードは無視するポリシーなのですね」と書かれていた。そういう特別な条件においてどう判断するかはソフトウェアによって異なり、生物学者にとって都合の悪い判断をすることもありうる。例えば植物のミトコンドリアゲノムから転写された mRNA は

Complete Sequence of a 641-kb Insertion of Mitochondrial DNA in the Arabidopsis thaliana Nuclear Genome Genome Biol Evol. 2022 May 3;14(5):evac059. doi: 10.1093/gbe/evac059. Peter D Fields など PMID: 35446419 PMCID: PMC9071559

というように、特別な条件を備えているものが多い。ソフトウェアを作成した研究者はそんなことは考慮していないだろうから、どんな結果が出るのかは調べないとわからない。

またソフトウェアに与えるパラメーターが自分の用途には不適切な値なことに気がつかずに、ずっと使ってしまうこともあり得る。複数のソフトウェアを使えるようにすれば、そういう問題がおきにくくなるだろう。

そこで stringtie と同様な機能を持つ、subread というソフトウェアパッケージを試してみる。

 $ conda install subread

でインストールした。これは本体がインストールされるだけなので、ライブラリーが古いものに入れ替わって困るようなことはない。

丁寧に説明をしていただいている解説記事がネットから見つかるので、それを真似る。 またマニュアルを読む。 https://bi.biopapyrus.jp/rnaseq/analysis/expression/featurecounts.html   上坂先生による解説 http://kazumaxneo.hatenablog.com/entry/2017/07/11/114046 

subread パッケージの一員として featureCounts というコマンドがある。 これに対して、bam ファイルと gtf ファイルを与えて処理させる。

このとき使用する gtf ファイルは、GTF2 フォーマットでなければならない。Stringtie などで使用した、TAIR から入手したファイルは GFF3 フォーマットだったので、そのままでは使えなかった。 これは gffread というソフトウェアでコンバートできた。参考にしたページ: http://ccb.jhu.edu/software/stringtie/gff.shtml

gffread は Bioconda からインストールできる。 > conda install gffread

 > gffread -E /場所/TAIR10_GFF3_genes.gtf -T -o- > TAIR10_GTF2_genes.gtf

解説の通りに行ってみる。以前、上に書いた方法でつくっておいた bam ファイルをそのまま使ってみる。 featureCounts は sam ファイルを入力にしてもよいので bam に変換しなくてもよい。

 > featureCounts -p -T4 -t exon -g gene_id -a TAIR10_GTF2_genes.gtf -o 〜〜.txt  SRR〜〜〜〜〜.bam

画面には「ペアエンドデータ」「87.4 %が assign された」などの情報が出る。このソフトウェアは大変優秀で、処理はすぐに終了する。 結果は、ただ一つのテキストファイルにまとめられて出力される。ファイルが多数生成すると何かと面倒なので、まとめられている方が都合がよい。

出力されたテキストファイルを見ると、

Gene ID (AGI コード)、染色体の名前、遺伝子の開始位置、終了位置、塩基鎖の+−、遺伝子領域の長さ、 と並んでいて、最後にその領域にマップされたリードのカウント数が記入されている。 この結果は生のカウント数なので、全リード数の違い、高発現している遺伝子のカウント数が変化することによる偏りなどの影響を補正しないといけない。それらは、featureCounts の出力を DEseq2 などのソフトウェアで処理することで行う。

オプション -M, -O による違い

featureCounts には「複数回アライメントされたpaired-endを全部カウントするオプション(-M)」というオプションがある。 これをつけるとどうなるかを調べてみる。

これをつけると、カウントされる率が 87.4 % から 92.4 % に高くなった。 シロイヌナズナではミトコンドリアゲノムにコードされた遺伝子 (ATMGxxxxx) の配列と、とてもよく似た遺伝子が核ゲノムにも存在する。オプションをつけないと、ATMG にマッピングされるカウント数はとても少ない。-M をつけると数が増えた。

「-O 1 つの領域が複数の feature id で定義されている場合がある。featureCounts は、このような領域にマッピングされるリードを集計していない。-O を付けることで、このような複数の feature で定義されている領域にマッピングされるリードも集計する。」これをつけるとどうなるかを調べてみる。

これをつけると、カウントされる率が 87.4 % から 89.3 % に高くなった。アンチセンス RNA ができるような遺伝子や、二つの遺伝子の末端がオーバーラップしている場合だと、一つの領域が二つの遺伝子に対応する可能性がある。

「-M -O」ではどうか。これをつけると、カウントされる率が 95.4 % に高くなった。 ミトコンドリアゲノムにコードされた遺伝子 (ATMGxxxxx) にマッピングされるカウント数の数が、-M だけより少し増えた。

このようにソフトウェアのマニュアルをよく読んで適切なオプションを選ぶことが結果に影響することがわかった。

featureCounts の結果を DESeq2 で後処理する

DESeq2 は R言語で書かれたソフトウェアである。R 言語で書かれたスクリプトは Rstudio などから実行できるが、

 > Rscript ファイル名.R

で実行することもできる。Rscript を使うと自動化しやすくなる。

インストールは Rstudio から行う。2020/01 に、R version 3.6.3 から、BiocManager::install("DESeq2") を実行した。Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3 (2020-02-29) とコンソール画面に出力された。 一回目はエラーが出た。もう一回実行すると、一回目にうまく追加されなかった必要なライブラリが追加され、エラーは出なくなった。

DESeq2 についてはたくさんの優れた解説が見つかるので、それらを真似てみる。

入力ファイルには、上で作った featureCounts の結果を使う。結果はテキストファイルになっているので、まず読み込む。

このときに問題があった。featureCounts の結果は、行の最初から文字列である gene_id などの情報が書かれていて、同じ行の最後にカウント数が書かれている。これを R 言語で読み込むとカウント数も文字列になる。そのままでは処理しにくいのでカウント数だけを取り出して数値に変換して、行列にしないといけない。

R 言語では、あるデータがどんな種類かを mode(〜) で判別できる。データフレーム、リスト、行列、文字列、数値などたくさんあり、変換が必要になることがある。

例: 

 mat_data <- matrix( unlist(リスト), ncol = 12 )
 mat_data <- matrix( as.numeric(mat_data), ncol = 12 )

データが行列になっていれば、test <- test[ c(-2, -3, -4, -5) ] のように必要ない列を取り除いたり、必要な列だけを取り出すことが簡単にできる。

データを行列の数値にできれば、後は多数の優れた解説が公開されているのでそれらに従ってスクリプトを作る。 DESeq2 で FPKM 値を簡単に計算するためには featureCounts の出力に含まれている遺伝子領域の長さを使うとよい。その値をオブジェクト dds にカラム basepairs として追加する。

 # 遺伝子領域の長さだけを取り出した行列を作る(この場合 2 番目のカラムに入っていた)
 genelength <- matrix( unlist(sampledata [2]), ncol = 1) 
 # Data を set する
 dds <- DESeqDataSetFromMatrix(countData = mat_data, colData = group, design =~ con )
 mcols(dds)$basepairs <- as.numeric(genelength)

オブジェクト dds は S4 クラスというデータ型の実体(インスタンス)になっている。S4 型のオブジェクトにカラムを追加するには、mcols() という関数を使う必要がある。カラム basepairs を設定すると fpkm(dds) で FPKM 値を計算できる。

 group <- data.frame(con = factor(c("WT", "EtOH", "DCMC", "DBmib", "WT", "EtOH", "DBMib", "DCMC" ) ) ) ## GSE96041 ラベルが間違っている
 condition1 <- "DBmib"        # 処理区
 condition2 <- "EtOH"        # 対照区
 num_of_sample <- nrow(group)

 # 二グループの Wald 検定
 dds <- DESeq(dds)
 res <- results(dds, contrast = c("con", condition1, condition2 ) )
 # geneid, genelength, count, FPKM を追加
 res$geneid <- geneid
 res$genelength <- genelength
 res$count <- counts(dds)
 res$FPKM <- fpkm(dds)
 #
 write.table(res, file="DESeq2_result.txt", sep = "\t", row.names = F, col.names = T, quote = F )

データにつけられているラベルに間違いがあることがある

GEO データベースなどで公開されているデータは、作成した研究者が間違ったラベルをつけたとしてもエラーが出るわけではない。たまに、「SRR5330630 と SRR5330631 はどう見てもラベルが入れ違っている」というようなことがある。それに気がつくには FPKM 値、または生のカウント数を出力に含めて、自分の目で見て確認することが有効である。

Stringtie と featureCounts の比較:featureCounts のほうが良い

私が手元のデータで試した結果では、どちらも同じような結果になる。featureCounts はファイルやディレクトリを設定する手間が少なくてすみ、高速に処理してくれるので featureCounts の方が良いような気がする。

Stringtie の結果を DESeq2 で後処理することもできることが紹介されている。   https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual   

この場合 Stringtie のオプションを raw reads 生のカウント数を出力するようにしないといけないらしい。

RNAseq のデータをマイクロアレイのデータと混ぜて区別せずに扱いたい

シロイヌナズナでは Affymetrix などのマイクロアレイのデータが膨大に蓄積され、きわめて有用なリソースになっている。それらは RMA 法などで処理され、遺伝子ごとに発現量が対数で表されている。また対照と処理区の間の発現変化の比率を対数で表した値を使うことが多い。

RNAseq で得られる値には FPKM など何種類かある。それぞれどんな特徴があるか。

(つづく)

The Integrative Genomics Viewer (IGV)

https://software.broadinstitute.org/software/igv/

https://igv.org/app を開くと、オンライン版の IGV を使用できる。

左上の「Genome」で、分析の対象になるゲノムを選ぶ。私は「Arabidopsis TAIR10」にした。

その右に「Tracks」がある。ここで BAM ファイルを読み込む。単なる BAM ファイルではダメで、「位置でソートされ、インデックスがついている」ファイルにしないといけない。それには samtools を使うことができる。SAM ファイルは非常にサイズが大きいが、変換すると 1/4 くらいに小さくできる。

 samtools view -bS 〜〜〜.sam  >  〜〜〜.bam

 samtools sort 〜〜〜.bam  -o  〜〜〜.sorted.bam

 samtools index 〜〜〜.sorted.bam

これで 〜〜〜.sorted.bam ファイルと 〜〜〜.sorted.bam.bai ファイルが生成する。

これらの二つのファイルを同時に「Tracks」で読み込む。読み込まれてから処理が進むのに少し時間がかかる。左上で染色体を選べる。左右に伸びるバーがあり、染色体上の位置を選択できる。右上にズームバーがある。ある程度ズームすると、遺伝子の配置と、それらにマップされたリードが図示される。

RNAseq 受託分析を利用してみた

サンプル: シロイヌナズナ芽生え シャーレ寒天培地で育成 処理区と対照で、3条件、2回ずつ  DEseq でデータを処理する際には、一つの条件で最低2回測定値が必要になる。

RNA の抽出: セパゾール RNA Super G (ナカライテスク)を用いて抽出して、イソプロパノール沈殿を RNeasy Plant mini (キアゲン)で精製して調製した。ナノドロップ(と同じような機器)で吸光度のスペクトルを確認した。

DNase I 処理: 各サンプル RNA を約 15 マイクログラム、チューブに取り DNase I 処理を行った。DNase I はインビトロジェンの、バッファーと EDTA が付属されているものを用いた。反応は全量が 40 マイクロリットルで、DNase を 2 マイクロリットル加えた。反応は 25 ℃で 10 分間にした。植物では RNA サンプルへの DNA の混入は少ないので、これで問題なかった。

EDTA を加えて反応を止めて、純水を加えて 100 マイクロリットルにした。すぐに RNeasy カラムで精製して DNase を除去して RNA を回収した。この方法で RNA は 10 マイクログラムくらい回収できた。

DNase I 処理に使う RNA 量が少ないと回収率が悪くなってしまい、よくない。

チェック: DNase 処理した RNA をナノドロップ(と同じような機器)で定量、スペクトルをチェックした。一部を使って RT-qPCR を行い、PCR がうまく増え DNase 処理する前と同じ結果が得られることを確認した。

受託解析: 生物技研という会社のサービスを利用した。RNAseq では「RNA のクオリティチェックが大切です」と必ず書かれているがバイオアナライザなどが必要になる。この会社では RNA を送ればクオリティチェックもしてくれるので都合がよかった。優れたサービスを行う企業は他にもたくさんある。見積、支払いやサンプルに関する連絡などはメールで行えた。送付などに関する注意事項は会社の Web ページに丁寧に解説されていたのでその通り行った。ドライアイスは 3 kg を発泡スチロールのボックスに詰めた。

結果: 結果は USB メモリに書き込まれて送られてきた。私は上に書いたように fastq データを処理できるようにしたので fastq のデータを送ってもらった。受託解析を行う各社で処理してもらうこともできる。さまざまなデータ分析サービスも提供されているので今後予算があれば利用してみたい。fastq データを fastqc, HISAT2, featureCounts, DEseq で処理したところ、質のよいデータが得られている(約 98 % がマッピングされていた、RT-qPCR で得ていた結果とよく一致した)ことが確認できた。 シーケンシングに使われている機器はイルミナ社のものではないが、fastq ファイルのフォーマットはイルミナ社のものと同じなので同様に処理できる。

早朝の会津若松駅ライブカメラから 

vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12: