kuroの覚え書き

96の個人的覚え書き

RNAseqデータの解析 RでcummeRbundその２

Science R

サンプル間で有意な遺伝子発現の差（p<0.05）を持つ遺伝子数をグラフィカルにあらわす

> mySigMatrix <- sigMatrix(cuff, level = 'genes', alpha = 0.05)
> mySigMatrix


> library("psych")
> setwd("/Users/hogehoge/expression/cuffdiff_results")
> cuff_data <- read.table('/Users/hogehoge/expression/cuffdiff_results/genes.fpkm_tracking', header = T)
> cuff_data2 <- cuff_data[c(10,14,18,22,26,30)]
> names(cuff_data2) <- c("P0","P2","P5","P7","P10","Adult")
> cuff_data2$P0 <-cuff_data2$P0 + 1
> cuff_data2$P2 <-cuff_data2$P2 + 1
> cuff_data2$P5 <-cuff_data2$P5 + 1
> cuff_data2$P7 <-cuff_data2$P7 + 1
> cuff_data2$P10 <-cuff_data2$P10 + 1
> cuff_data2$Adult <-cuff_data2$Adult + 1
> log10_data <- cuff_data2
> log10_data$P0 <- log10(cuff_data2$P0)
> log10_data$P2 <- log10(cuff_data2$P2)
> log10_data$P5 <- log10(cuff_data2$P5)
> log10_data$P7 <- log10(cuff_data2$P7)
> log10_data$P10 <- log10(cuff_data2$P10)
> log10_data$Adult <- log10(cuff_data2$Adult)
> library("psych")
> pairs.panels(log10_data[1:6], ellipses = FALSE, scale = FALSE, density = FALSE, smooth = FALSE, pch = ".")
> pairs.panels(log10_data[1:6], ellipses = TRUE, scale = TRUE, density = TRUE, smooth = TRUE, pch = ".")


上のマトリクスで表した遺伝子発現に差のある遺伝子を集めてgene setを作る

> mySigGeneIds<-getSig(cuff,alpha=0.05,level='genes')
> head(mySigGeneIds)
[1] "XLOC_000297" "XLOC_001707" "XLOC_010198" "XLOC_014130" "XLOC_014726" "XLOC_014731"
#できたgene setを覗いてみると・・・
> length(mySigGeneIds)
[1] 454
#できたgene setの中に含まれる遺伝子の数
> mySigGenes <- getGenes(cuff, mySigGeneIds)
Getting gene information:
	FPKM
	Differential Expression Data
	Annotation Data
	Replicate FPKMs
	Counts
Getting isoforms information:
	FPKM
	Differential Expression Data
	Annotation Data
	Replicate FPKMs
	Counts
Getting CDS information:
	FPKM
	Differential Expression Data
	Annotation Data
	Replicate FPKMs
	Counts
Getting TSS information:
	FPKM
	Differential Expression Data
	Annotation Data
	Replicate FPKMs
	Counts
Getting promoter information:
	distData
Getting splicing information:
	distData
Getting relCDS information:
	distData
 警告メッセージ: 
1: Closing open result set, pending rows 
2: Closing open result set, pending rows 

> mySigGenes
CuffGeneSet instance for  454  genes
 
Slots:
	 annotation
	 fpkm
	 repFpkm
	 diff
	 count
	 isoforms	 CuffFeatureSet instance of size 1066 
	 TSS		 CuffFeatureSet instance of size 693 
	 CDS		 CuffFeatureSet instance of size 386 
	 promoters		 CuffFeatureSet instance of size 454 
	 splicing		 CuffFeatureSet instance of size 693 
	 relCDS		 CuffFeatureSet instance of size 454

このgene setを使ってヒートマップを作ると

> h <- csHeatmap(mySigGenes, cluster = "both")
Using tracking_id, sample_name as id variables
No id variables; using all as measure variables
> print(h)


左の遺伝子名のフォントサイズとか変更する方法は要調査

サンプル間の相似性をグラフィカルに表す

> myDistHeat<-csDistHeat(genes(cuff))
> myDistHeat


Adultが一つだけ他と違うということは明らかだが、どういう発現パターンがあるのかをk-meanクラスタリングで解析する

> ic<-csCluster(mySigGenes,k=6)
Using tracking_id, sample_name as id variables
> icp <- csClusterPlot(ic)
> icp

試しに6つのクラスターに分類してみる

  1. P0で発現していて時間経過とともに低下し、Adultではほぼ発現していない
  2. P0~P10ではほとんど発現しておらず、Adultのみで発現している
  3. P2で一過的に発現し、すぐに低下する
  4. 徐々に発現が上昇しAdultでは発現している
  5. 徐々に発現が上昇するがAdultではほとんど発現していない
  6. 初期のみ発現していて低下していく（発現のピークはP2)

k=の値を変えると違ったパターンが見えてくるかもしれない。