maSigPro 处理时间序列转录组数据

2020-12-10 本文已影响0人 FengSL

时间序列研究的是基因表达的动态行为，测量的是一系列和时间点之间有强烈相关性的过程。和针对某一时间点的基因表达进行差异分析不同，时间序列更加关注是发现基因表达的趋势，以有助于理解生物学动态变化过程（比如对刺激的反应、发育过程、周期行为等）。也就是说，时间序列关注的是整体变化趋势而不是某特异表达。这样就有几个挑战，一是要分析的数据量会很大，二是实验条件变多，三是要发掘实验动态变化本质，传统的统计学方法比如t-tests就无能为力了，需要运用新的统计学方法，四是样本间的时间间隔并不总是相等。

maSigPro的全称是Microarray Significant Profiles，采用2步回归策略。第一步选择基因，首先，可以使用统计学程序来鉴定显著表达变化的基因；第二步选择变量，把随时间变化发生显著表达变化的基因进行聚类并且可视化。这个可以通过回归来解决，其中，时间被视为数量变量，实验条件视为分类变量，这样就可以分析趋势变化。

下面简要概述maSigPro分析的步骤、参数及原理。

1. 读入表达数据以及实验设计

data("data.abiotic")
head(data.abiotic)
data("edesign.abiotic")
head(data.abiotic)

data.abiotic为一个芯片表达数据集，每一行为一个基因，每一列为一个样品/条件/重复。maSigPro同样可以处理RNA-seq数据，但是maSigPro没有集合数据标准化的算法，数据读入之前需要首先经过标准化的处理，可以将edge R得到的cpm作为maSigPro的读入数据集。
edesign.abiotic为实验设计表，行名为data.abiotic（读入数据）的列名，Time为时间信息，Replicate为生物学重复信息，生物学重复之间使用相同的数字表示。后面的列比较灵活，用于表示实验设计，有3种常见设计：

1. 考察单一时间变量下基因表达变化（无处理）

所在列均为1即可，如

Time = rep(c(1,5,10,24),each=3)
Replicates = rep(c(1:4),each=3)
Group = rep(1,12)
edesign = cbind(Time,Replicates,Group)
rownames(edesign) = paste("Array",c(1:12), sep= "")

2. 不同处理具有相同的时间起点

data("edesignCT")
head(edesignCT)

不同样品在相同时间点均为1即可。

3. 多种处理组合

在用maSigPro进行分析时，一般情况都是两个或两个以上的感兴趣的变量，其中一个典型的就是时间变量，另外一个通常都是分类变量，代表实验组别（比如不同的处理，细胞株，组织等）。

edesign.abiotic

对应处理列为1，非对应为0。

2. 确定回归模型

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

maSigPro采用多项式回归，参数degree设置多项式使用的次数，在本示例中有3个时间点，degree设置为2，多项式次数设置过高会导致过拟合，一般在能够解释自变量和因变量关系的前提下，次数应该越低越好。

3. 寻找显著基因

这时maSigPro两步法的第一步，计算得到每个基因的回归拟合，并挑选显著基因用于后续分析。

fit = p.vector(data.abiotic, d, Q=0.05, MT.adjust = "BH",counts = FALSE)

Q为FDR阈值，MT.adjust为多重检验方法。counts为布尔值，当为False时表示采用正太分布对芯片数据进行回归拟合，当为True时表示采用负二项分布对Counts数据进行回归拟合。
这一步将筛选出至少在一个group下表现出差异的基因。

4. 寻找显著差异

tstep = T.fit(fit,step.method="backward")

T.fit将采用逐步回归法，查找基因在哪些group下具有显著差异。

5. 获得显著基因集

sigs = get.siggenes(tstep, rsq = 0.7, vars = "groups")

rsq是R-squared的阈值，默认是0.7，作者文中认为当有三个生物学重复时，0.5已经可以得到非常收敛的结果。
vars可以接受3个参数:"groups","all"和"each"。
groups：每一个实验设计groups产生一个列表，包括reference groups以及处理和control之间的差异。
all：只有一个列表，包含所有的模型变异。
each：每一种可能的变异为一个列表，包括control和处理之间，以及处理或control和时间尺度之间的差异。
运行masigpro 主要有四步：

5. 结果获取

1. 冷处理条件下，基因表达随时间的变化

data("data.abiotic")
head(data.abiotic)
data.control = data.abiotic[,c(10:18)]  ## 提取冷处理数据
data("edesign.abiotic")
head(edesign.abiotic)
edesign.control = edesign.abiotic[c(10:18),c(1:2,4)]  ##构建冷处理实验设计
edesign.control

d= make.design.matrix(edesign.control,degree = 2)
fit = p.vector(data.control,d,Q=0.05,counts=FALSE)
fit$i  ## 差异基因数
tstep = T.fit(fit)
sigs = get.siggenes(tstep,rsq=0.5,vars="groups")
names(sigs$sig.genes)
pdf("Figure 1.pdf")
cl = see.genes(sigs$sig.genes$Cold,newX11 = F,alfa = 0.05, k = 5)## 对差异结果进行聚类分析
dev.off()
cluster = as.data.frame(cl$cut)  ## 提取聚类信息
head(cluster)

Figure 1

2. 多组合处理

d = make.design.matrix(edesign.abiotic,degree = 2)
fit = p.vector(data.abiotic,d,Q=0.05,counts=FALSE)
fit$i
tstep = T.fit(fit)
sigs = get.siggenes(tstep,rsq=0.5,vars="groups")
names(sigs$sig.genes)
pdf("Figure 2.pdf")
cl = see.genes(sigs$sig.genes$ColdvsControl,newX11 = F,alfa = 0.05,k=5)
dev.off()

Figure 2