macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

RNA seq 非モデル生物の解析

ゲノム情報がない場合、まず最初にアセンブリから始めることになる。RNAアセンブリはTrinityが有名なのでtrinityを使ってみる。練習データとして、SRAに登録されているレタスのRNA seqデータを使う。

 

http://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA004542

36サンプルある。スクリプト組むほどではないので、手動でfastqを順番にダウンロードする。ペアでなくフラグメントデータぽい。

 

1、アセンブル

trinityでアセンブルを試みる。 

フラグメントデータなら

Trinity --seqType fq --single 1.fq --max_memory 10G --CPU 18

 --single;ペアードエンドシーケンスでなくフラグメントシーケンスデータの時

 --seqType;fastqならfq

 --max_memory ピークメモリー使用量の上限を指定

 --CPU : 使用するCPU数

 

 

ペアードエンドデータなら

Trinity --seqType fq --left R1.fq --right R2.fq --max_memory 10G --CPU 18

という感じに行う。R1.fqとR2.fqがペアードエンドデータになる。

 

複数のペアードエンドデータなら

(例えば下の3データ)

1_R1.fq x 1_R2.fq

2_R1.fq x 2_R2.fq

3_R1.fq x 3_R2.fq

fastqを, (カンマ) でつないで以下のようにランする。

Trinity --seqType fq --left 1_R1.fq,2_R1.fq,3_R1.fq --right 1_R2.fq,2_R2.fq,3_R2.fq --max_memory 35G --CPU 18

 

ただし、今回のSRAデータを対象に実行したところ、ファイル名が認識できずスタートできなかった。原因を探ったところ以下のような記事が出てきた。

ヘッダーの仕様が古いせいっぽい。試しに割と新しいMiSeqのシーケンスデータを使ったところ、ちゃんとアセンブルできた。まあ下の記事のようにシェルかperlなどを使ってヘッダーを切り替えることのはルールさえわかれば簡単だろうけど...

How to change fastq reads header for running Trinity on them?

例えば

ヘッダーをオリジナルの

@DRR057317.1 MG00HS09:691:C8JNWACXX:7:1101:1476:1994 length=51

とかから

@M02077:53:000000000-ANJF9:1:1101:18340:000$a 1:N:0:15";

に変え($aカウンタで数値を1リードずつ増やすすスクリプトで処理)、さらに3行目は新しいfastqの定義に従い+に変えて実行した。すると正常にランできるようになった。

 

 ただし今回はde novo RNA seqのワークフローを確認するために進めているだけなので、Trinityを使わずアセンブラのrnaspadesを使うことにする。まずはアセンブル前に解凍する。

bunzip2 -c DRR057309.fastq.bz2 > DRR057309.fastq

全サンプル順次解凍する。

 

rnaspadesでアセンブルする。以下のコマンドを打った。

rnaspades.py -t 30 -k auto -o all_data -s DRR057059.fastq -s \
DRR057060.fastq -s DRR057061.fastq -s DRR057062.fastq -s \
DRR057063.fastq -s DRR057064.fastq -s DRR057065.fastq -s \
DRR057066.fastq -s DRR057067.fastq -s DRR057068.fastq -s \
DRR057069.fastq -s DRR057070.fastq -s DRR057071.fastq -s \
DRR057072.fastq -s DRR057073.fastq -s DRR057074.fastq -s \
DRR057075.fastq -s DRR057076.fastq -s DRR057077.fastq -s \
DRR057078.fastq -s DRR057079.fastq -s DRR057080.fastq -s \
DRR057081.fastq -s DRR057082.fastq -s DRR057083.fastq -s \
DRR057084.fastq -s DRR057085.fastq -s DRR057086.fastq -s \
DRR057087.fastq -s DRR057088.fastq -s DRR057089.fastq -s \
DRR057090.fastq -s DRR057091.fastq -s DRR057092.fastq -s \
DRR057093.fastq -s DRR057094.fastq

 

 

全部で36サンプル。n=1、サイズは全て1GBくらい。メモリは200GBほどしかない。うまくできるか?と思いながら待っていたがダメだった。k-merを自動から指定値に変更する。リードは50bpの長さがある。rnaspadesは複数のk-merを指定できないので、奇数のk=21にして実行した。

rnaspades.py -t 30 -k 21 -o all_data_k21 -s DRR057059.fastq -s \
DRR057060.fastq -s DRR057061.fastq -s DRR057062.fastq -s \
DRR057063.fastq -s DRR057064.fastq -s DRR057065.fastq -s \
DRR057066.fastq -s DRR057067.fastq -s DRR057068.fastq -s \
DRR057069.fastq -s DRR057070.fastq -s DRR057071.fastq -s \
DRR057072.fastq -s DRR057073.fastq -s DRR057074.fastq -s \
DRR057075.fastq -s DRR057076.fastq -s DRR057077.fastq -s \
DRR057078.fastq -s DRR057079.fastq -s DRR057080.fastq -s \
DRR057081.fastq -s DRR057082.fastq -s DRR057083.fastq -s \
DRR057084.fastq -s DRR057085.fastq -s DRR057086.fastq -s \
DRR057087.fastq -s DRR057088.fastq -s DRR057089.fastq -s \
DRR057090.fastq -s DRR057091.fastq -s DRR057092.fastq -s \
DRR057093.fastq -s DRR057094.fastq

 

今度は無事解析が終わった。解析時間は6時間くらいだった。できたtranscripts.fastaの中を見てみる。

tail transcripts.fasta 
>NODE_134906_length_32_cov_12.4528_g133128_i0 CCCCTCGCCCATCCTCCTTTCTGTCGAATGTC >NODE_134907_length_31_cov_11.2115_g133129_i0 TTCCAAAGAATCCGAATCTTCGAATTCAAGT >NODE_134908_length_31_cov_4.61538_g133130_i0 CGACACAAAGACAGTAACCAATCAATGGACT >NODE_134909_length_31_cov_3.28846_g133131_i0 TACTGGTCTCATCTCTGTATGAAGGTGGGTG >NODE_134910_length_31_cov_2.63462_g133132_i0 TGATCGAAAAAATGGTGTTAGTCTCATCTGT

 

 

13万4910分子種ある。解析が重くなるので短かすぎるリードを除く。grep -n "length_99" transcripts.fasta

で99bpの行番号を特定し、

head -n transcripts.fasta > transcripts100over.fasta

として100bp以上を保存。nは99bpが始まる1つ前の行番号。spadesの出力はサイズ順なのでこれて過不足ない。

100bp以上に限定すると分子種は134910から86450に減った。流れを確認するトレーニングなので、とりあえずこれでマッピング解析を進める。

 

2、マッピング

#index作成

bowtie2-build -f transcripts_more_than_100bp.fa transcripts_100bp

#mapping

tophat -p 8 -o tophat_DRR057059 --no-novel-juncs transcripts_100bp DRR057059.fastq &
.
.
.
tophat -p 8 -o tophat_DRR057094 --no-novel-juncs transcripts_100bp DRR057094.fastq &

 

省略するが、上のように合計36サンプルを順次解析した。

コア数12、メモリ64GBのmac proでは4ジョブ投げるとCPUが100%に張り付いた。やはり重い。4つずつ実行したらそれだけで半日かかった。

 

 

3、cufflinksによるFPKM計算

cufflinksを使う。cufflinksは各遺伝子のリード数をカウントしてFPKMを算出したり、データを統合して比較することができるプログラム。

まずtophatのマッピングデータを使い発現量をカウントする。

cufflinks -o  cufflinks_output_0h tophat_DRR057059/accepted_hits.bam & cufflinks -o   cufflinks_output_2h tophat_DRR057060/accepted_hits.bam & cufflinks -o   cufflinks_output_4h tophat_DRR057061/accepted_hits.bam & cufflinks -o   cufflinks_output_6h tophat_DRR057062/accepted_hits.bam & cufflinks -o   cufflinks_output_8h tophat_DRR057063/accepted_hits.bam & cufflinks -o  cufflinks_output_10h tophat_DRR057064/accepted_hits.bam & cufflinks -o  cufflinks_output_12h tophat_DRR057065/accepted_hits.bam & cufflinks -o  cufflinks_output_14h tophat_DRR057066/accepted_hits.bam & cufflinks -o  cufflinks_output_16h tophat_DRR057067/accepted_hits.bam & cufflinks -o  cufflinks_output_18h tophat_DRR057068/accepted_hits.bam & cufflinks -o  cufflinks_output_20h tophat_DRR057069/accepted_hits.bam & cufflinks -o  cufflinks_output_22h tophat_DRR057070/accepted_hits.bam & cufflinks -o  cufflinks_output_24h tophat_DRR057071/accepted_hits.bam & cufflinks -o  cufflinks_output_26h tophat_DRR057072/accepted_hits.bam & cufflinks -o  cufflinks_output_28h tophat_DRR057073/accepted_hits.bam & cufflinks -o  cufflinks_output_30h tophat_DRR057074/accepted_hits.bam & cufflinks -o  cufflinks_output_32h tophat_DRR057075/accepted_hits.bam & cufflinks -o  cufflinks_output_34h tophat_DRR057076/accepted_hits.bam & cufflinks -o  cufflinks_output_36h tophat_DRR057077/accepted_hits.bam & cufflinks -o  cufflinks_output_38h tophat_DRR057078/accepted_hits.bam & cufflinks -o  cufflinks_output_40h tophat_DRR057079/accepted_hits.bam & cufflinks -o  cufflinks_output_42h tophat_DRR057080/accepted_hits.bam & cufflinks -o  cufflinks_output_44h tophat_DRR057081/accepted_hits.bam & cufflinks -o  cufflinks_output_46h tophat_DRR057082/accepted_hits.bam & cufflinks -o  cufflinks_output_48h tophat_DRR057083/accepted_hits.bam & cufflinks -o  cufflinks_output_50h tophat_DRR057084/accepted_hits.bam & cufflinks -o  cufflinks_output_52h tophat_DRR057085/accepted_hits.bam & cufflinks -o  cufflinks_output_54h tophat_DRR057086/accepted_hits.bam & cufflinks -o  cufflinks_output_56h tophat_DRR057087/accepted_hits.bam & cufflinks -o  cufflinks_output_58h tophat_DRR057088/accepted_hits.bam & cufflinks -o  cufflinks_output_60h tophat_DRR057089/accepted_hits.bam & cufflinks -o  cufflinks_output_62h tophat_DRR057090/accepted_hits.bam & cufflinks -o  cufflinks_output_64h tophat_DRR057091/accepted_hits.bam & cufflinks -o  cufflinks_output_66h tophat_DRR057092/accepted_hits.bam & cufflinks -o  cufflinks_output_68h tophat_DRR057093/accepted_hits.bam & cufflinks -o  cufflinks_output_70h tophat_DRR057094/accepted_hits.bam &

 

 

計算結果をcuffmergeでmergeしてgtfファイルを作る。そのために、上で行なった36個の解析結果のパスを書いたテキストファイルを準備し、cuffmergeにそれを読み込ませる。ファイル名はmerge.txtとし、以下のように記載した。

 

>cat merge.txt #中身をチェック
./cufflinks_output_0h/transcripts.gtf ./cufflinks_output_2h/transcripts.gtf ./cufflinks_output_4h/transcripts.gtf ./cufflinks_output_6h/transcripts.gtf ./cufflinks_output_8h/transcripts.gtf ./cufflinks_output_10h/transcripts.gtf ./cufflinks_output_12h/transcripts.gtf ./cufflinks_output_14h/transcripts.gtf ./cufflinks_output_16h/transcripts.gtf ./cufflinks_output_18h/transcripts.gtf ./cufflinks_output_20h/transcripts.gtf ./cufflinks_output_22h/transcripts.gtf ./cufflinks_output_24h/transcripts.gtf ./cufflinks_output_26h/transcripts.gtf ./cufflinks_output_28h/transcripts.gtf ./cufflinks_output_30h/transcripts.gtf ./cufflinks_output_32h/transcripts.gtf ./cufflinks_output_34h/transcripts.gtf ./cufflinks_output_36h/transcripts.gtf ./cufflinks_output_38h/transcripts.gtf ./cufflinks_output_40h/transcripts.gtf ./cufflinks_output_42h/transcripts.gtf ./cufflinks_output_44h/transcripts.gtf ./cufflinks_output_46h/transcripts.gtf ./cufflinks_output_48h/transcripts.gtf ./cufflinks_output_50h/transcripts.gtf ./cufflinks_output_52h/transcripts.gtf ./cufflinks_output_54h/transcripts.gtf ./cufflinks_output_56h/transcripts.gtf ./cufflinks_output_58h/transcripts.gtf ./cufflinks_output_60h/transcripts.gtf ./cufflinks_output_62h/transcripts.gtf ./cufflinks_output_64h/transcripts.gtf ./cufflinks_output_66h/transcripts.gtf ./cufflinks_output_68h/transcripts.gtf ./cufflinks_output_70h/transcripts.gtf  

mergeコマンドの実行。

cuffmerge -o all_merged merge.txt

-o 出力ディレクトリ名

gtfファイルが出力される。

 

 

4、発現変動遺伝子の分析

cufflinksのコマンドcuffdiffを使い発現変動を調べる。

 

cuffmergeの出力結果をinputとし、以下のようにコマンドを実行。

> cuffdiff -o cuffdiff_out -p 20 all/merged.gtf -L 0h,2h,4h,6h,8h,10h,12h,14h,16h,18h,20h,22h,24h,26h,28h,30h,32h,34h,36h,38h,40h,42h,44h,46h,48h,50h,52h,54h,56h,58h,60h,62h,64h,66h,68h,70h \
tophat_DRR057059/accepted_hits.bam tophat_DRR057060/accepted_hits.bam \
tophat_DRR057061/accepted_hits.bam tophat_DRR057062/accepted_hits.bam \
tophat_DRR057063/accepted_hits.bam tophat_DRR057064/accepted_hits.bam \
tophat_DRR057065/accepted_hits.bam tophat_DRR057066/accepted_hits.bam \
tophat_DRR057067/accepted_hits.bam tophat_DRR057068/accepted_hits.bam \
tophat_DRR057069/accepted_hits.bam tophat_DRR057070/accepted_hits.bam \
tophat_DRR057071/accepted_hits.bam tophat_DRR057072/accepted_hits.bam \
tophat_DRR057073/accepted_hits.bam tophat_DRR057074/accepted_hits.bam \
tophat_DRR057075/accepted_hits.bam tophat_DRR057076/accepted_hits.bam \
tophat_DRR057077/accepted_hits.bam tophat_DRR057078/accepted_hits.bam \
tophat_DRR057079/accepted_hits.bam tophat_DRR057080/accepted_hits.bam \
tophat_DRR057081/accepted_hits.bam tophat_DRR057082/accepted_hits.bam \
tophat_DRR057083/accepted_hits.bam tophat_DRR057084/accepted_hits.bam \
tophat_DRR057085/accepted_hits.bam tophat_DRR057086/accepted_hits.bam \
tophat_DRR057087/accepted_hits.bam tophat_DRR057088/accepted_hits.bam \
tophat_DRR057089/accepted_hits.bam tophat_DRR057090/accepted_hits.bam \
tophat_DRR057091/accepted_hits.bam tophat_DRR057092/accepted_hits.bam \
tophat_DRR057093/accepted_hits.bam tophat_DRR057094/accepted_hits.bam  

 

-o write all output files to this directory 

-L  comma-separated list of condition labels

-p number of threads used during quantification

gtfファイルも相対パスで指定している。これが一番重くて4日かかった。クラミドモナスのtime courseサンプルのデータでは1時間くらいで終わったが、サンプル数が増えると処理時間が指数関数的に増えるのかもしれない。仕方がないので別のことをしながら待った。

 

 

f:id:kazumaxneo:20170406155216j:plain

解析が終わると指定したcuffdiff_out/に上のようなファイルができる。データベースのdbファイルは25GBもある。36サンプルもあると流石にでかい。

 

Rを起動して前回と同じくcummeRbundで解析してみる。

Rをターミナルで起動

library(cummeRbund) #ライブラリ読み込み
cuff <- readCufflinks("cuffdiff_out") #cufflinks解析結果のディレクトリ全体を読み込み
cuff #cuffコマンドを打って正常に読み込まれたか確認
 CuffSet instance with:  
 36 samples  
44858 genes  
44858 isoforms  
44858 TSS  
0 CDS  
28260540 promoters  
28260540 splicing  
0 relCDS 

 

 

 

 

csBoxplot(genes(cuff)) #genesの分布をbox plotで表示

f:id:kazumaxneo:20170406154047j:plain

 

genediffオブジェクトに差のあるgenesを読み込む。

genediff <- diffData(genes(cuff)) 

統計的に優位なものだけ取ってくる。(significant列がyesのgene)

genediff2 <- genediff[(genediff$significant == 'yes'),]

続けてk.means法で分類。

genediffdataID2 <- genediff2[,1]
genediffdataID2 <- t(genediffdataID2)
uesakaGenes2 <-getGenes(cuff, genediffdataID2) #ここが重たい
k.means <-csCluster(uesakaGenes2, k=4) #k=300-400を越えるとRがクラッシュした
csClusterPlot(k.means)  #グラフを描画

 

f:id:kazumaxneo:20170406190327j:plain

  

genediff3 <- genediff[((genediff$status == 'OK')& (genediff$significant == 'yes')),]
genediffdataID3 <- genediff3[,1]
genediffdataID3 <- t(genediffdataID3)
uesakaGenes3 <-getGenes(cuff, genediffdataID3) 
k.means <-csCluster(uesakaGenes3, k=20) 
csClusterPlot(k.means)  

 

f:id:kazumaxneo:20170406185428j:plain

 

大量のデータがあるのでk.means法では適切な分類は難しいように思える。

やり方を変え、特定の遺伝子と同じ発現パターンをもつ遺伝子を探すことを考える。

 

そのため、まず全遺伝子のfpkmを出力し、変動パターンを確認して自分が興味ある変動パターンを示している遺伝子を見つけ、それと似た発現パターンの遺伝子をサルベージすることを考える。

 

全遺伝子のデータだけを取り出す

genes <- genes(cuff) 
is(genes)
fpkm <- fpkmMatrix(genes)  
write.table(fpkm, file="fpkm.csv", sep=",") #csv形式で保存

excelcsvを開き、グラフを見て興味のある変動パターンを示す遺伝子を目視で調べていく。

f:id:kazumaxneo:20170406203848j:plain

例えばグラフの青の遺伝子は振動が徐々になくなっていっている。これは明暗条件からLL条件に移したサンプルと論文に書いてあるので、明暗の条件がなくなったことで周期を失う遺伝子と推測される。これと似た発現パターンをもつ遺伝子をcummeRbundで提供されているコマンドを使って絞り込んでいく。

 

 

my.similar <- findSimilar(cuff, "XLOC_000003", n = 20)

上のコマンドの意味だが、マニュアルにはFor example, if you wanted to find the 20 genes most similar to "PINK1", you could do the following: と書いてある。XLOC_000003と似た発現パターンを示すtop20遺伝子をfindSimilarコマンドで探してmy.similarに代入している。

 

散布図を描く

my.similar.expression <- expressionPlot(my.similar, logMode = T, showErrorbars = F) 
my.similar.expression

 

f:id:kazumaxneo:20170406204333j:plain

振動が徐々に消えていくようなパターンの発現パターンの遺伝子top20を絞り込めた(あくまでXLOC_000003に似た発現パターンの遺伝子を探しただけなことに注意)。

 

 

 

非モデル生物のRNA seqについては、様々なアプローチが考えられると思います。下のリンクでもう少し考えています。


結論としては、近縁種のゲノム情報が利用できてもde novo assemblyの方がorf検出率ははるかに上です。de novo assemblyしてcDNA配列を得て進めて行くことになるのではないでしょうか?

私が試みた限り、同じ属のゲノム情報を使いcufflinksでゲノムガイドアセンブルしても、ハウスキーピング遺伝子が10ほどしか検出できませんでした。