image.png
转录组数据差异基因的分析是转录组的核心，针对不同的数据类型有不同的分析方法，之前我们也提到过。这次要做的是芯片数据的分析，采用limma包，方法是通用的（适用于mRNA以及小RNA的芯片数据）。
一般数据挖掘分析都是下载表达矩阵进行分析，当然也可以从原始数据开始，感兴趣的可以查看limma包的详细介绍。
图片
读入表达矩阵，对其log转化（至于为和要log化，可以百度下，如果作者上传的数据是已经log过的，则无需处理），简单QC一下，用boxplot看看数据样本之间的一致性。
setwd("E:/生物信息学")
A <- read.table("GSE125424_series_matrix.txt",header = T,
row.names = 1,sep = "\t",comment.char = "!")
A <- log2(A)
boxplot(data.frame(A))
image.gif
这个数据还可以，一致性挺好。但是通常时候，样本之间的差异还是挺大的，就需要用如下的代码进行标准化再进行差异分析。使用wateRmelon包进行标准化的处理。
#library("wateRmelon")
#A<- betaqn(A)
#boxplot(A)
#write.table(A,file="norm.txt",sep="\t",quote=F)
limma包分析差异基因还是比较简单的，我们首先需要构建分组，并构建比较矩阵，这里只有两组，虽然默认就是两组比较，但是还是需要定义，明确是哪个比哪个。
library(limma)
library(dplyr)
f <- c(rep("control", 3), rep("Treat", 3)) %>% as.factor()
desigN <- model.matrix(~ 0 + f)
colnames(desigN) <- levels(f)
fit = lmFit(A, desigN)
contrast.matrix <- makeContrasts(Treat - control, levels = colnames(coef(fit)))
contrast.matrix
# Contrasts
#Levels Treat - control
# control -1
# Treat 1
差异分析：
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tempOutput <- topTable(fit2,coef=1,n=Inf,adjust="BH")
head(tempOutput)
# logFC AveExpr t P.Value adj.P.Val B
#1446537_at 3.848140 2.164047 22.71094 2.092108e-06 0.09435616 1.339931
#1440714_at -3.426720 1.887167 -17.83508 7.254608e-06 0.09706875 1.172877
#1429297_at -3.198857 1.599429 -17.80127 7.325561e-06 0.09706875 1.171274
#1443604_at -3.069908 1.616981 -16.22950 1.176281e-05 0.09706875 1.086811
#1421100_a_at -2.234104 1.150563 -16.19503 1.189145e-05 0.09706875 1.084715
#1438524_x_at 2.943123 1.471561 15.71083 1.388776e-05 0.09706875 1.054005
最后将差异分析文件导出：
write.csv(tempOutput, file="DEGs.csv")
导出的差异分析结果就可以进行下一步的筛选了，设定不同的阈值筛选差异基因，以及基因功能富集等下游分析了。
网友评论