[step by step]小白做网络分析（3）

本帖的目的是指出GEO-2R的R scripts中哪些代码需要做改动。见标红部分。

##   Differential expression analysis with limma

setwd("F:\\4.本科毕业课题\\2.DEGs") ##设置路径

library(Biobase)

library(GEOquery)

library(limma)

## load series and platform data from GEO ##不需改动

gset <- getGEO("GSE66660", GSEMatrix =TRUE, AnnotGPL=TRUE)

if (length(gset) > 1) idx <- grep("GPL11532", attr(gset, "names")) else idx <- 1

gset <- gset[[idx]]

## make proper column names to match toptable   ##不需改动

fvarLabels(gset) <- make.names(fvarLabels(gset))

## group names for all samples  ##不需改动

gsms <- "11001100"

sml <- c()

for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }

## log2 transform（RMA方法） ##采用RMA方法则不需改动，可改为mas5等方法

ex <- 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) ||

          (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

if (LogC) { ex[which(ex <= 0)] <- NaN

  exprs(gset) <- log2(ex) }

## set up the data and proceed with analysis

sml <- paste("G", sml, sep="")    ## set group names

fl <- as.factor(sml)

gset$description <- fl

design <- model.matrix(~ description + 0, gset)

colnames(design) <- levels(fl)

fit <- lmFit(gset, design)

cont.matrix <- makeContrasts(G1-G0, levels=design)

fit2 <- contrasts.fit(fit, cont.matrix)

fit2 <- eBayes(fit2)  ##原文为eBayes（fit2,0.01） 

tT <- topTable(fit2, adjust="BH", sort.by="B", number=dim(gset)[1])  ##adjust.method是P.value成为adj.P.value所采用的方法，number是tT中包含的基因数目（可自行设置，我这个就是显示所有）

tT <- subset(tT, select=c("ID","adj.P.Val","P.Value","t","B","logFC","Gene.symbol","Gene.title")) ##不要改动顺序

write.table(tT, paste("all genes(T).xls"),quote=FALSE, sep="\t") ##原文为write.table(tT, file=stdout(), row.names=F, sep="\t")，两者均可，个人喜欢改动后的，简单易懂

#######一定要注意：上面保存的是全部基因！！下面我进行差异基因筛选：p<0.05,|logFC|>2。是原文里没有的！

degup<-tT[tT[,"adj.P.Val"]<0.05,] 

degup<-degup[degup[,"logFC"]>2,]   ##上调基因

write.table(degup,paste("GSE66600DEGs(upregulate).xls"),quote=FALSE,sep="\t")

degdown<-tT[tT[,"adj.P.Val"]<0.05,] 

degdown<-degdown[degdown[,"logFC"]<（-2）,] ##下调基因

write.table(degdown,paste("GSE66600DEGs(downregulate).xls"),quote=FALSE,sep="\t")