rm(list = ls()) ###清除所有的环境变量
options(stringsAsFactors = F)
load(file = '../input.Rdata') 加载已保存的R文件 ../表示上级目录，类似linux的用法
suppressMessages(library(limma))取消加载包时出现的文字信息
design <- model.matrix(~0+factor(group_list))
model.matrix 线性回归时自动生成设计矩阵
image
image
然而实际上，我们在limma学习过程中，并不是真的在进行线性回归，我们在做的工作叫做差异分析！（上述公式并看不懂0-0）
- 差异分析：
分组矩阵（设计矩阵）将数据按照一定的要求进行分组，然后比较不同组之间的差异
比较矩阵：比如说哪一些数据更重要，哪一些不重要的权重划分以及怎么比较等等
因子：limma差异分析需要转化为因子，而不是应用数字。即分组的自变量可以是发病或不发病，自变量转为因子就可以将其作为多维坐标上的一个点来进行画图及后续分析。 -
design <- model.matrix(~0+factor(group_list))
design <- model.matrix(~factor(group_list))
从回归的角度来讲，～group_list是对group_list进行回归，剩下的内容是残差，而～0是对0这个变量进行回归，也就是认为没有变量，剩下的group内容是一个常数项 - limma
三个数据：表达矩阵 分组矩阵（设计矩阵） 差异比较矩阵
三个步骤：lmFit eBayes topTable - 管道函数，类似linux中“ | ”
来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤，可以大量减少内存中的对象，节省内存。符号%>%，这是管道操作，其意思是将%>%左边的对象传递给右边的函数，作为第一个选项的设置（或剩下唯一一个选项的设置），减少中间步骤的赋值
#CLL包中里面有一个sCLLex数据，它是ExpressionSet格式
library(CLL) #or suppressPackageStartupMessages(library(CLL))
sCLLex #这个数据集中有22个样本，12625个基因；使用exprs可以查看这个数据集的表达水平；
data(sCLLex)
class(sCLLex)
exprSet=exprs(sCLLex)#提取出来表达矩阵，赋给data_expression
##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)#查看样本编号
pdata=pData(sCLLex)# 查看分组信息
varMetadata(sCLLex) #查看所有表型变量
data_phenotype=pData(sCLLex)#提取表型信息
featureNames(sCLLex)[1:10] #查看基因芯片编码
featureNames(sCLLex) %>% unique() %>% length() # 查看是否有重复的编码
##上面的%>%这个符号相当于Linux中的管道命令
#直接运行报错，需要安装magrittr or dplyr
#BiocManager::install("dplyr")
#library("dplyr")
data_expression <- exprs(sCLLex) #提取出来表达矩阵，赋给data_expression
View(data_expression)#可以直观看出来这个表达矩阵
group_list=as.character(pdata$Disease) #从数据框中只要表型信息
table(group_list) #统计表型信息
group_list=as.character(pdata[,2])
dim(exprSet)
# [1] 12625 22
exprSet[1:5,1:5]
##### 绘制芯片数据的质量值，如下所示：
library(ggplot2)
library(reshape2)
library(gpairs)
library(corrplot)
y <- melt(as.data.frame(data_expression))
p <- ggplot(data=y,aes(x=variable,y=value))
p <- p + geom_boxplot(aes(fill=variable))
p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
p <- p + xlab("分组信息") + ylab("表达值") + guides(fill=FALSE)
p
本文大范围参考简书天涯清水/ctomclancy 及生信菜鸟团的帖子，感谢
网友评论