一文拿捏，GEO数据下载处理和差异分析

BiocManager::install("GEOquery")install.packages("tidyverse")install.packages("ggplot2")BiocManager::install("sva")BiocManager::install("limma")BiocManager::install("preprocessCore")

library(GEOquery)library(tidyverse)library(ggplot2)library(sva)library(preprocessCore)library(limma)

#数据下载，只下载表达矩阵参数：GSEMatrix=TGSEdata <- getGEO("GSE1919", GSEMatrix=T, AnnotGPL=FALSE)eset <- GSEdata[[1]]#提取表达矩阵,通过exprs函数提取，行名为探针ID，列名为样本名exprSet <- as.data.frame(exprs(eset))# 提取GSM ID作为样品名，pData函数来获取样本信息GSMID <- pData(eset)[,c(1,2,8,10,12)]#利用fData函数提取基因注释信息annotation<-fData(eset)[,c(1,11,12)]#有些探针组对应多个基因，使用“///”分割，保留第一个基因annotation[,2]<-data.frame(sapply(annotation[,2],function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]annotation[,3]<-data.frame(sapply(annotation[,3],function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]# 把整理好的探针组ID、Gene symbol、ENTREZ ID对应关系保存到文件里colnames(annotation) <- c("probe_id", "gsym", "entrez_id")write.csv(annotation, "gplTOgene.csv", row.names = F)#把表达矩阵的探针ID换成Gene Symbolannotation <- read.csv("gplTOgene.csv")#将探针信息和表达矩阵合并exprSet$probe_id <- rownames(exprSet)express <- merge(x = exprSet, y = annotation[,c(1:2)], by = "probe_id", all.x = T)#去除probe_id这列信息express$probe_id <- NULL# 有些gene对应多个探针组，此处保留其中表达量最大值rowMeans <- apply(express, 1, function(x) mean(as.numeric(x), na.rm = T))express <- express[order(rowMeans, decreasing =T),]express <- express[!duplicated(express[, ncol(express)]),]  #express的最后一列为gene nameexpress <- na.omit(express)#将最后一列作为行名rownames(express) <- express[,ncol(express)]#删除最后一列express <- express[,-(ncol(express))]#筛选正常样本和RA样本表达矩阵，去除了OA样本。GSE1919<-select(express,-c("GSM34393","GSM34394","GSM34395","GSM34396","GSM34397"))

#查看表达量数量级,通常log2的值应该在0-18之间GSE1919<-log2(GSE1919+1)#均一化处理GSE1919.norm <- normalize.quantiles(as.matrix(GSE1919))colnames(GSE1919.norm) <- colnames(GSE1919)rownames(GSE1919.norm) <- rownames(GSE1919)##limma进行差异分析#设置分组group <- rep(c("normal","RA"),           times=c(5,5))#生成模型设计矩阵design <- model.matrix(~0+factor(group))#设置行名colnames(design) <- levels(factor(group))#设置列名rownames(design) <- colnames(GSE1919.norm)#构建比较矩阵contrast.matrix <- makeContrasts(paste0(unique(group),collapse = "-"),                              levels = design)#线性模型拟合fit <- lmFit(GSE1919.norm,design)#计算estimated coefficients和标准误fit2 <- contrasts.fit(fit,contrast.matrix)#进行差异分析fit2 <- eBayes(fit2)#提取所有差异基因，通过pvalue值来排序temp <- topTable(fit2,coef = 1,n=Inf,sort.by="P")#去除空值DEG <- na.omit(temp)##将行名转化为列名DEG$symbol <- rownames(DEG)#保存差异分析结果文件write.table(DEG, "normal_RA_DEG.txt",sep = '\t', quote = FALSE, col.names = T, row.names = FALSE)