30hのデータがおかしいので、それ以外の時間で分析を行う。tophatでのマッピングは其の２で行ったので、cuffdiffの定量だけやり直せばよい。これだけでかなり時間を短縮できる。
cuffdiff -o cuffdiff_except-_for_30h -p 20 Chlamydomonas_reinhardtii.v3.1.34.gtf -L 0h-1,0h-2,2h,4h,8h,12h,18h,24h,45h,60h 0h-1_out/accepted_hits.bam 0h-2_out/accepted_hits.bam 2h_out/accepted_hits.bam 4h_out/accepted_hits.bam 8h_out/accepted_hits.bam 12h_out/accepted_hits.bam 18h_out/accepted_hits.bam 24h_out/accepted_hits.bam 45h_out/accepted_hits.bam 60h_out/accepted_hits.bam
Rを起動
library(cummeRbund)
cuff <- readCufflinks("cuffdiff_out")
genediff <- diffData(genes(cuff))
genediff2 <- genediff[(genediff$significant == 'yes'),]
write.table(genediff2, file="result3.csv", sep=",")
genediffdataID2 <- genediff2[,1]
genediffdataID2 <- t(genediffdataID2)
uesakaGenes2 <-getGenes(cuff, genediffdataID2)
k.means <-csCluster(uesakaGenes2, k=9)
k.means.plot <- csClusterPlot(k.means)
k.means.plot
k=9
k=6
k=4
k=9の 結果をよくよく見ると、30h以外にも24で発現が大きく落ちる遺伝子グループがある（１のウィンドウ）。もしかしたら30hの応答も生物的な応答かもしれない。 とりあえずk=9では十分なグループ分けができていない可能性もある。一度k=12で解析してみた。
k=12 
time 4h、12h、24h、45h、60hでもそれぞれ大きく落ち込む遺伝子グループがある。0hはそれ以前の変動がわからないが0h特異的に発現が少ない遺伝子群かもしれない。
k=20まで上げてみよう。
K=20
まだまだパターンがある気がするし、人間の目で見てかぶったパターンも増えたように思える。nitrate responseのサンプルなので,0hから発現がグンと上がっている遺伝子が気になる。k=20のウィンドウ１２の遺伝子をまとめた。
Nar1;2という硝酸同化系遺伝子が見事に応答している。非常に妥当な結果である。
ここで最初に立ち返ってt=30のデータはどうなのかを考える。例えばRNAの分解が進んだサンプルのrRNA/mRNA比率が代わり、特に低発現の遺伝子の発現パターンにばらつきが出るようになる。まあそれはそうだろう。分解スピードが同じと仮定し細胞内に10コピー存在するmRNAと100コピー存在するmRNAを考えると、先に0になるのは10コピーのmRNAの方だろう。このようなサンプルからRNAをとってシーケンスすると、分解してないサンプルと比べた時の比率が無限大になる。たとえFPKM 0.1と仮定して数値を加えたとしても、分解していないサンプルとの比率を正確に求めることはできない。また、分解しきっていなくても、低発現遺伝子の定量がばらつく要因はいくらでも考えられる。というのも、RNAの分解もシーケンス時の逆転写、PCR、シーケンスも（世の中のこと全て）確率的にしか進行しないため、初発の数値が低いまま進むとブレが大きくなるからだと解釈できる。例えると、サイコロを6回しか振らないと、6の目が３回出て、１は１回も出ない（オレすごい！）ようなことが起こるが、サイコロを6万回ふるとどの目も限りなく1/6に近い値がでる（オレ普通だった....）。それでも1/6に近づかないなら、サイコロが正６面体になっていないか重心がずれている可能性がある。
まずは色々な時間で散布図を書いてみる。
csScatter(genes(cuff),"X2h","X4h",smooth=T)
csScatter(genes(cuff),"X2h","X8h",smooth=T)
csScatter(genes(cuff),"X12h","X18h",smooth=T)
csScatter(genes(cuff),"X18h","X45h",smooth=T)
csScatter(genes(cuff),"X45h","X60h",smooth=T)
csScatter(genes(cuff),"X24h","X45h",smooth=T)
この中では45hが他と一番違うように感じる。さて、では30との比較だとどうなるか？
csScatter(genes(cuff),"X24h","X30h",smooth=T)
csScatter(genes(cuff),"X18h","X30h",smooth=T)
csScatter(genes(cuff),"X2h","X30h",smooth=T)
思ったほどばらつきはでかくなっていない。まあ対数表記だからか。しかし、どの時間と比較しても左下の方に弧を描くように描画されているプロット出る。説明できないけどアーティファクトっぽいなあ。
k=9の結果をヒートマップ化する。
k.means
出てきた数値をMiにコピー。エクセルに読み込んでフィルターをかけ、欲しいウィンドウだけ取り出してMiで保存。Rに再読み込み。
ウィンドウ１、8、3の遺伝子でヒートマップを出してみる。
csv <- read.table("1_3_8.txt" , header=F, sep="\t")
csv2 <- t(csv)
csv2data <-getGenes(cuff, csv2)
csHeatmap(csv2data)
ウィンドウ6、7の遺伝子でヒートマップを出してみる。
csv <- read.table("1_3_8.txt" , header=F, sep="\t")
csv2 <- t(csv)
csv2data <-getGenes(cuff, csv2)
csHeatmap(csv2data)
思ったほどまとまらない。k=20の19と20のウィンドウだとどうだろうか？
数が少ないのもあってさすがにこれは綺麗にまとまっている。
