转录组数据分析之时序分析(maSigPro包)

    在转录组数据分析过程中,我们最常做的是不同处理方式的样本之间的比较(Treated vs Control),这时候我们采用“DEG分析+pathway分析”的方式就可基本完成对数据的分析。但偶尔我们也会碰到一类特殊的数据,即同一种处理在不同的时间点可能呈现不同的生物学特征,因此我们需要分析不同基因随着时间的动态变化。这一类分析即所谓的“时序分析”。

    那么如何对转录组数据进行时序分析呢?这里小编给大家介绍一款简单好用的R包——maSigPro包。

该包的分析流程如下图所示:

该包最初由Conesa等人于2006年开发,于2014年经Nueda等人进一步优化。maSigPro采用二步回归的策略对数据进行分析,并且可以找出不同处理组之间的差异基因。并且在第二次回归模型中可以将具有相似表达模式的基因进行聚类并可视化。

    需要注意的是,该包不会对数据进行标准化处理,因此需要自己提前进行预处理哦。

其分析代码如下:

1. 加载maSigPro包及演示数据

library(maSigPro)# load maSigPro library

data(data.abiotic)

data(edesign.abiotic)

数据格式如下:

> edesign.abiotic

> data.abiotic

2. 定义回归模型

design <- make.design.matrix(edesign.abiotic, degree = 2)

其中degree参数与时间点的数量有关,时间点较多时degree数值应相应增加(该演示数据中为3个时间点)。design为列表形式,其中group.vector可以看回归参数是如何分配的。

> design$groups.vector

[1] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"       "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control"     

[9] "ColdvsControl" "HeatvsControl" "SaltvsControl"

3. 差异基因分析

通过p.vector()函数对每一个基因进行回归分析,同时给出每一个基因的p-value及p.adjust(MT.adjust、Q等参数)。

fit <- p.vector(data.abiotic, design, Q = 0.05, MT.adjust = "BH", min.obs = 20)

由此返回以下参数:

> 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()函数可以进一步分析不同组别的差异基因之间的不同特征,该步采用逐步回归的方式(step.method参数)。

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()查看差异基因结果

> sigs$sig.genes$SaltvsControl$g

[1] 191

该函数可以将表达模式相似的基因进行聚类并可视化。其中k参数为cluster的数量,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()函数可以查看某一特定基因的表达情况。

 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)

怎么样,是不是很简单呢?快来试一试吧。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,383评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,522评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,852评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,621评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,741评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,929评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,076评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,803评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,265评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,582评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,716评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,395评论 4 333
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,039评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,027评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,488评论 2 361
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,612评论 2 350

推荐阅读更多精彩内容