爾云間? 一個(gè)專門做科研的團(tuán)隊(duì)

以GSE70493為例進(jìn)行GEO2R分析所使用R代碼的學(xué)習(xí)

代碼如下：

# 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 limma
library(GEOquery)
library(limma)
library(umap)

#此處說明使用limma進(jìn)行差異表達(dá)分析，需要加載的一些R包，包括GEOquery，limma，umap。

# load series and platform data from GEO
gset <- getGEO("GSE70493", GSEMatrix =TRUE, AnnotGPL=FALSE)
if (length(gset) > 1) idx <- grep("GPL17586", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

#此段代碼是下載GEO數(shù)據(jù)集，數(shù)據(jù)集名稱為gset

# make proper column names to match toptable
fvarLabels(gset) <- make.names(fvarLabels(gset))

#此段代碼是生成gset數(shù)據(jù)集的列名，用到了make.names函數(shù)，目的是為了定義有效的名稱

# group membership for all samples
gsms <- "111010111101100000111110011010001000001010000011110111100100101"
sml <- strsplit(gsms, split="")[[1]]

#此段代碼是對數(shù)據(jù)集中所有的樣本進(jìn)行成員關(guān)系分析 ，首先是定義了一個(gè)組，然后使用字符串去進(jìn)行分割并進(jìn)行提取，其中的[[1]]，表示提取某個(gè)列表中第一個(gè)子表中的所有元素，一般我們提取列表其中的元素對象 盡量采用[[]]直接提取，不容易出錯(cuò)。

# log2 transformation
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)
if (LogC) { ex[which(ex <= 0)] <- NaN
? exprs(gset) <- log2(ex) }

#此段代碼是對表達(dá)矩陣ex進(jìn)行l(wèi)og2 對數(shù)值轉(zhuǎn)換，通常使用limma處理時(shí)，需要經(jīng)過log2后的矩陣作為表達(dá)矩陣輸入。

# assign samples to groups and set up design matrix
gs <- factor(sml)
groups <- make.names(c("control","treat"))
levels(gs) <- groups
gset$group <- gs
design <- model.matrix(~group + 0, gset)
colnames(design) <- levels(gs)

fit <- lmFit(gset, design)? # fit linear model

#此段代碼是對樣本進(jìn)行分組，進(jìn)行表達(dá)矩陣設(shè)計(jì)

# set up contrasts of interest and recalculate model coefficients
cts <- paste(groups[1], groups[2], sep="-")
cont.matrix <- makeContrasts(contrasts=cts, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)

#此段代碼是建立比對關(guān)系并重新計(jì)算模型回歸系數(shù) ，具體用到的函數(shù)有paste，makeContrast等

# compute statistics and table of top significant genes
fit2 <- 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")

#此段代碼是統(tǒng)計(jì)差異顯著的基因，這里用到了topTable函數(shù)，通過fdr的方法進(jìn)行矯正排序，篩選了250個(gè)差異基因。