典型医学设计实验GEO数据分析(step-by-step)-数据获取到标准化介绍了实验的设计、数据获取、数据标准化和注释（点蓝字，助回忆），下面是如何利用Limma和线性模型鉴定差异基因，并进行GO富集分析。
现在这套代码综合高通量数据中批次效应的鉴定和处理-系列总结和更新，已经迁移到高颜值免费在线绘图工具，新增WGCNA和差异分析。
线性模型为了分析发炎和未发炎组织的差异表达，我们需要构建一个线性模型。线性模式是实验数据分析的常用方法，适用于几乎任何复杂的实验设计。后面我们专门出文介绍，推荐MikeLove和MichaelIrizzary的genomicsclass和可汗学院的线性代数免费公开课。
我们这里主要用limma包构建线性模型进行差异表达分析。这个包可以同时比较很多实验组并且尽量维持其易用性。首先对每个基因的表达拟合一个线性模型，然后用经验贝叶斯(EmpiricalBayes)或其他方法进行残差分析获得合适的t统计量，并针对小样本实验的方差估计进行优化，使得分析结果更加可靠。
在构建线性模型时，采用缩写UC和CD代表疾病类型，non_infl.和infl.代表发炎与否。
#获得所有的个体信息individual-as.character(Biobase::pData(palmieri_final)$Characteristics.individual.)#组织信息的空格替换为下划线tissue-str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.phenotype.,"","_")#简化组织的描述信息tissue-ifelse(tissue=="non-inflamed_colonic_mucosa","nI","I")#疾病信息替换下划线，并简化描述disease-str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.disease.,"","_")disease-ifelse(str_detect(Biobase::pData(palmieri_final)$Factor.Value.disease.,"Crohn"),"CD","UC")
因为要找发炎和未发炎组织的差异表达基因，所以理论上只需要比较这两个变量。但因为每个独立的个体有两套芯片检测结果(发炎和未发炎组织)，这是一个配对样品实验设计。下游分析时需要考虑个体差异的影响。如果一个实验特征对结果可能有系统性影响，需要对引入这个特征作为阻断因子(bolckingfactors)。
为了与文章一致，我们为每个疾病分别构建了一个设计矩阵。(也可以针对完整数据集设计一个联合模型，但两种疾病可能特征差别很大，放在一起可能不太合适，从典型医学设计实验GEO数据分析(step-by-step)-数据获取到标准化文中的PCA结果也可以看出来)
#获得CD疾病的个体i_CD-individual[disease=="CD"]#获得两种组织类型和个体的矩阵#0+表示不设置截距项，所有样品都有自己的回归系数design_palmieri_CD-model.matrix(~0+tissue[disease=="CD"]+i_CD)colnames(design_palmieri_CD)[:]-c("I","nI")rownames(design_palmieri_CD)-i_CDi_UC-individual[disease=="UC"]design_palmieri_UC-model.matrix(~0+tissue[disease=="UC"]+i_UC)colnames(design_palmieri_UC)[:]-c("I","nI")rownames(design_palmieri_UC)-i_UC
检查下设计矩阵:
head(design_palmieri_CD[,:6])##InIi_CD83i_CD4i_CD09i_CD55############head(design_palmieri_UC[,:6])##InIi_UC44i_UCi_UC99i_UC######4##4##0##0
设计矩阵(designmatrix)中行代表每个芯片，列代表囊括入线性模型的变量，包含是否发炎，个体信息。i_UC44是病人44的变量，UC代表病人所患疾病。矩阵中的0和代表所属关系(也称激活状态)。
在线性模型中，第一个个体(设计矩阵第一行)会作为其他个体的基准，不会包含到样品变量中。这里~0是去除截距项，每个样品都计算回归系数。
除了像上面把个体作为blockingfactor的方式，也可以构建混合模型(mixedmodel)，增加randomeffect，以后再详细讲述。
单个基因差异表达分析测试先选择一个基因查看其分布和拟合出的线性模型，这里选择的是PROBEID为，genesymbol为CRAT的基因。
首先看下其表达，整体在未发炎样品中高。而且不同病人间基因的绝对表达丰度差别挺大。如果我们没有在设计矩阵中考虑到这个问题，线性模型就会把这些数据视为一个整体。考虑到每个个体的基准表达水平不同，最终获得的差异倍数会有较高的方差。
tissue_CD-tissue[disease=="CD"]crat_expr-Biobase::exprs(palmieri_final)["",disease=="CD"]crat_data-as.data.frame(crat_expr)colnames(crat_data)[]-"org_value"crat_data-mutate(crat_data,individual=i_CD,tissue_CD)crat_data$tissue_CD-factor(crat_data$tissue_CD,levels=c("nI","I"))ggplot(data=crat_data,aes(x=tissue_CD,y=org_value))+geom_boxplot(aes(fill=tissue_CD))+geom_line(aes(group=individual,color=individual))+ggtitle("ExpressionchangesfortheCRATgene")+theme(legend.position="none")
我们拟合线性模型计算回归系数。
crat_coef-lmFit(palmieri_final[,disease=="CD"],design=design_palmieri_CD)$coefficients["",]crat_coef##InIi_CD83i_CD4i_CD09i_CD55i_CD55i_CD86##6...38-0..-0..-.##i_CDi_CDi_CD3i_CD36i_CDi_CD37i_CDi_CD##-0.-0.6.-0.-0.-0.-0.-0.
设计矩阵(design_palmieri_CD)与回归系数向量(crat_coef)相乘获得拟合后的表达值。
crat_fitted-design_palmieri_CD%*%crat_coefrownames(crat_fitted)-names(crat_expr)colnames(crat_fitted)-"fitted_value"crat_fitted##fitted_value##64_I_.CEL7.9##64_II.CEL6.##83_I.CEL7.36##83_II.CEL6.89##4_I.CEL6.##4_II.CEL6.##09_A.CEL7.##09_B.CEL7.34##55_I.CEL6.##55_II.CEL6.##55_I.CEL7.##55_II.CEL7.##86_I.CEL6.9##86_II.CEL5.##_I.CEL6.##_II.CEL5.##_I.CEL7.##_II.CEL6.##3_I.CEL7.09##3_II.CEL6.##36_I.CEL6.##36_II.CEL6.4##_I.CEL6.##_II.CEL6.##37_I.CEL6.##37_II.CEL6.49##_I.CEL6.##_II.CEL6.##_I.CEL7.##_II.CEL6.
设计矩阵中每一行只有值为的变量用于计算拟合的表达值，crat_fitted的每一项代表每个样品各个变量回归系数的和。
例如对病人样品4_I.CEL6.3:对应的激活变量的回归系数的和“nI”(7.97)and“i_CD4”(-0.5)。
可视化发炎和未发炎组织拟合后的表达值:
crat_data$fitted_value-crat_fittedggplot(data=crat_data,aes(x=tissue_CD,y=fitted_value))+geom_boxplot(aes(fill=tissue_CD))+geom_line(aes(group=individual,color=individual))+ggtitle("ExpressionchangesfortheCRATgene")+theme(legend.position="none")
可以看到，所有病人中发炎组织和未发炎组织的CRAT基因表达差别一致，并且是由变量I(6.)和nI(7.97)的回归系数的差值决定。
这是因为一个病人的两个样本的相同blockingvariable在设计矩阵中都是激活的，只能在病人内比较。这些blockingvariable校正拟合后的组织特异性表达值趋向每个病人的表达值，因此最终的估计是个体内差异的平均值，也就是对应基因的差异倍数(logfoldchange)。(Theseblockingvariablescorrectthefittedtissuespecificexpressionvaluestowardstheexpressionlevelsoftheindividualpatients.Thereforethefinalestimateislikeanaverageofallthewithin-individualdifferences.)
CRAT基因差异表达检测
为了检测基因是否是差异表达，需要执行零假设(nullhypothesis)为基因在发炎和未发炎的样品中表达无差异的t-test。我们的这种设计概念上类似于配对t-test，统计量t的计算如下：
$$t=\frac{\bar{d}}{s/\sqrt{n}}$$
个体间表达值的差异的平均值。配对t-test的方差来源于配对样品。这与标准t-test不同，因此只要配对样品的表达是相关的，配对t检验就有更高的统计检出力(
