????在轉(zhuǎn)錄組數(shù)據(jù)分析過程中，我們最常做的是不同處理方式的樣本之間的比較（Treated vs Control）邪乍，這時(shí)候我們采用“DEG分析+pathway分析”的方式就可基本完成對(duì)數(shù)據(jù)的分析薪韩。但偶爾我們也會(huì)碰到一類特殊的數(shù)據(jù)穷蛹，即同一種處理在不同的時(shí)間點(diǎn)可能呈現(xiàn)不同的生物學(xué)特征偶惠，因此我們需要分析不同基因隨著時(shí)間的動(dòng)態(tài)變化丹莲。這一類分析即所謂的“時(shí)序分析”光坝。
????那么如何對(duì)轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行時(shí)序分析呢？這里小編給大家介紹一款簡(jiǎn)單好用的R包——maSigPro包甥材。
該包的分析流程如下圖所示：
該包最初由Conesa等人于2006年開發(fā)盯另，于2014年經(jīng)Nueda等人進(jìn)一步優(yōu)化。maSigPro采用二步回歸的策略對(duì)數(shù)據(jù)進(jìn)行分析洲赵，并且可以找出不同處理組之間的差異基因鸳惯。并且在第二次回歸模型中可以將具有相似表達(dá)模式的基因進(jìn)行聚類并可視化。
????需要注意的是叠萍，該包不會(huì)對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理芝发，因此需要自己提前進(jìn)行預(yù)處理哦。
其分析代碼如下：
1. 加載maSigPro包及演示數(shù)據(jù)
library(maSigPro)# load maSigPro library
data(data.abiotic)
data(edesign.abiotic)
數(shù)據(jù)格式如下：
> edesign.abiotic
> data.abiotic
2. 定義回歸模型
design <- make.design.matrix(edesign.abiotic, degree = 2)
其中degree參數(shù)與時(shí)間點(diǎn)的數(shù)量有關(guān)苛谷，時(shí)間點(diǎn)較多時(shí)degree數(shù)值應(yīng)相應(yīng)增加（該演示數(shù)據(jù)中為3個(gè)時(shí)間點(diǎn)）辅鲸。design為列表形式，其中g(shù)roup.vector可以看回歸參數(shù)是如何分配的腹殿。
> design$groups.vector
[1] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"?????? "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"?????
[9] "ColdvsControl" "HeatvsControl" "SaltvsControl"
3. 差異基因分析
通過p.vector()函數(shù)對(duì)每一個(gè)基因進(jìn)行回歸分析独悴，同時(shí)給出每一個(gè)基因的p-value及p.adjust（MT.adjust例书、Q等參數(shù)）。
fit <- p.vector(data.abiotic, design, Q = 0.05, MT.adjust = "BH", min.obs = 20)
由此返回以下參數(shù)：
> fit$i # returns the number of significant genes
> fit$alfa # gives p-value at the Q false discovery control level
> fit$SELEC # is a matrix with the significant genes and their expression values
4. 差異分析
通過T.fit()函數(shù)可以進(jìn)一步分析不同組別的差異基因之間的不同特征绵患，該步采用逐步回歸的方式（step.method參數(shù)）雾叭。
tstep <- T.fit(fit, step.method = "backward", alfa = 0.05)
5. 提取差異基因列表
sigs <- get.siggenes(tstep, rsq = 0.6, vars = "groups")
其中rsq為回歸模型R-squared的cutt-off值，vars用于表示分組變量落蝙。
6. Venn圖查看差異基因情況
suma2Venn(sigs$summary[, c(2:4)])
suma2Venn(sigs$summary[, c(1:4)])
7. see.genes()
采用see.genes()查看差異基因結(jié)果
> sigs$sig.genes$SaltvsControl$g
[1] 191
該函數(shù)可以將表達(dá)模式相似的基因進(jìn)行聚類并可視化织狐。其中k參數(shù)為cluster的數(shù)量，cluster.method為聚類的方式（hclust筏勒，kmeans移迫，Mclust）
see.genes(sigs$sig.genes$ColdvsControl, show.fit = T, dis =design$dis, cluster.method="hclust" ,cluster.data = 1, k = 9)
8. PlotGroups()
使用PlotGroups()函數(shù)可以查看某一特定基因的表達(dá)情況。
?STMDE66 <- data.abiotic[rownames(data.abiotic)=="STMDE66", ]
?PlotGroups (STMDE66, edesign = edesign.abiotic)
并且我們可以通過以下方式添加回歸曲線：
PlotGroups (STMDE66, edesign = edesign.abiotic, show.fit = T, dis = design$dis, groups.vector = design$groups.vector)
怎么樣管行，是不是很簡(jiǎn)單呢厨埋？快來(lái)試一試吧。
