GEO2R分析R代码学习之差异分析
收录于话题
GEO2R
以GSE70493为例进行GEO2R分析所使用R代码的学习：
# Version info: R 3.2.3, Biobase 2.30.0, GEOquery 2.40.0, limma 3.26.8#介绍使用的软件的版本，例如R是版本是3.2.3，Biobase的版本是2.30.0, GEOquery的版本是2.40.0, limma 的版本是3.26.8。
# Differential expression analysis with limmalibrary(GEOquery)library(limma)library(umap)
#此处说明使用limma进行差异表达分析，需要加载的一些R包，包括GEOquery，limma，umap。
# load series and platform data from GEOgset <- getGEO("GSE70493", GSEMatrix =TRUE, AnnotGPL=FALSE)if (length(gset) > 1) idx <- grep("GPL17586", attr(gset, "names")) else idx <- 1gset <- gset[[idx]]
#此段代码是下载GEO数据集，数据集名称为gset
# make proper column names to match toptablefvarLabels(gset) <- make.names(fvarLabels(gset))
#此段代码是生成gset数据集的列名，用到了make.names函数，目的是为了定义有效的名称
# group membership for all samplesgsms <- "111010111101100000111110011010001000001010000011110111100100101"sml <- strsplit(gsms, split="")[[1]]
#此段代码是对数据集中所有的样本进行成员关系分析 ，首先是定义了一个组，然后使用字符串去进行分割并进行提取，其中的[[1]]，表示提取某个列表中第一个子表中的所有元素，一般我们提取列表其中的元素对象 尽量采用[[]]直接提取，不容易出错。
# log2 transformationex <- exprs(gset)qx <- 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)if (LogC) { ex[which(ex <= 0)] <- NaNexprs(gset) <- log2(ex) }
#此段代码是对表达矩阵ex进行log2 对数值转换，通常使用limma处理时，需要经过log2后的矩阵作为表达矩阵输入。
# assign samples to groups and set up design matrixgs <- factor(sml)groups <- make.names(c("control","treat"))levels(gs) <- groupsgset$group <- gsdesign <- model.matrix(~group + 0, gset)colnames(design) <- levels(gs)fit <- lmFit(gset, design) # fit linear model
#此段代码是对样本进行分组，进行表达矩阵设计
# set up contrasts of interest and recalculate model coefficientscts <- paste(groups[1], groups[2], sep="-")cont.matrix <- makeContrasts(contrasts=cts, levels=design)fit2 <- contrasts.fit(fit, cont.matrix)
#此段代码是建立比对关系并重新计算模型回归系数 ，具体用到的函数有paste，makeContrast等
# compute statistics and table of top significant genesfit2 <- eBayes(fit2, 0.01)tT <- topTable(fit2, adjust="fdr", sort.by="B", number=250)tT <- subset(tT, select=c("ID","adj.P.Val","P.Value","t","B","logFC","SPOT_ID"))write.table(tT, file=stdout(), row.names=F, sep="t")
#此段代码是统计差异显著的基因，这里用到了topTable函数，通过fdr的方法进行矫正排序，筛选了250个差异基因。
扫码关注我们
shengxinguoer
生信果
生信硬核知识解答
和小果一起学生信