【万能代码分享】一键搞定Limma包非配对差异分析+热图+火山图+GO分析+KEGG分析+GSEA分析！

### 公众号：科研部### 公众号：科研部### 公众号：科研部### 该代码分为十部分：### 1.读入探针矩阵数据### 2.数据预处理：NA去除+log2转换+样本芯片矫正### 3.探针ID和基因名symbol转换### 4.基因Symbol表达矩阵获取：探针矩阵转换为基因Symbol矩阵+同名symbol保留表达量的最大值。### 5.使用limma包来做芯片的差异分析### 6.输出差异基因### 7. 作图，差异表达图，热图，火山图。### 8. 差异基因GO分析### 9. 差异基因KEGG分析### 10. 差异基因集的GSEA富集分析。

##################################################################################################### 0. 镜像##################################################################################################options(download.file.method = 'libcurl')options(url.method='libcurl')options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CR

##################################################################################################### 1.读入探针矩阵数据##################################################################################################exprSet <- read.table("1_GSE42872_series_matrix.txt",                      comment.char="!",                      stringsAsFactors=F,                      header=T)### comment.char="!" 意思是！开头的行作为注释
rownames(exprSet) <- exprSet[,1]exprSet <- exprSet[,-1]### 第一列变成行名

##################################################################################################### 2.数据预处理：NA去除+log2转换+样本芯片矫正#################################################################################################### NA的去除，公众号：科研部ex <- exprSetqx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))LogC <- (qx[5] > 100) ||  (qx[6]-qx[1] > 50 && qx[2] > 0) ||  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
## 自动判断是否需要log2转换，需要则进行log2转换if (LogC) {  ex[which(ex <= 0)] <- NaN  ## 取log2  exprSet <- log2(ex)  print("log2 transform finished")}else{  print("log2 transform not needed")}
## 样本芯片矫正library(limma)boxplot(exprSet,outline=FALSE, notch=T, las=2)#矫正之前单个样本值，此时平不平都无所谓。exprSet=normalizeBetweenArrays(exprSet)### 该函数默认使用quntile 矫正样本间差异。（公众号：科研部）boxplot(exprSet,outline=FALSE, notch=T, las=2)#矫正之后单个样本值，此时是平的。exprSet <- as.data.frame(exprSet)

##################################################################################################### 3.探针ID和基因名symbol转换#################################################################################################### 读入platformMap.txt文件，platformMap 中有常见的GPL平台和ID转换需要的R包（一一对应）。platformMap <- data.table::fread("1_platformMap.txt",data.table = F)
## 找到该平台对应的ID转换R包index <- "GPL6244"## 数据储存在bioc_package这一列中paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")
## 安装该R包#if (!require("BiocManager", quietly = TRUE))#  install.packages("BiocManager")#BiocManager::install("hugene10sttranscriptcluster.db")
## 加载该R包，并将探针ID转换成基因symbollibrary(hugene10sttranscriptcluster.db)probe2symbol_df <- toTable(get("hugene10sttranscriptclusterSYMBOL"))

##################################################################################################### 4.基因Symbol表达矩阵获取：探针矩阵转换为基因Symbol矩阵+同名symbol保留表达量的最大值。##################################################################################################library(dplyr)library(tibble)exprSet <- exprSet %>%  rownames_to_column("probe_id") %>%  inner_join(probe2symbol_df,by="probe_id") %>%  select(-probe_id) %>%  select(symbol,everything()) %>%  mutate(rowMean =rowMeans(.[,-1])) %>%  arrange(desc(rowMean)) %>%  distinct(symbol,.keep_all = T) %>%  select(-rowMean) %>%  column_to_rownames("symbol")
save(exprSet,file = "2_exprSet_rmdup.Rdata")

##################################################################################################### 5.使用limma包来做芯片的差异分析#################################################################################################### 1.创建分组group <- c(rep("con",3),rep("treat",3))#意思是在矩阵里面前面3个对照组，后面3个实验组group <- factor(group,levels = c("con","treat"))
## 样本PCA主成分分析：提前预测结果，看两组样本是否分开，分开则说明差异挺大的，可以进行差异分析。res.pca <- prcomp(t(exprSet), scale = TRUE)library(factoextra)fviz_pca_ind(res.pca,col.ind = group)
## 构建比较矩阵design <- model.matrix(~group)## 比较矩阵命名colnames(design) <- levels(group)design
## 2.线性模型拟合fit <- lmFit(exprSet,design)## 3.贝叶斯检验fit2 <- eBayes(fit)## 4.输出差异分析结果,其中coef的数目不能操过design的列数allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)## 这个数据很重要需要保存一下save(allDiff,file = "2_allDiff.Rdata")