參考學習GEO提供的GEO2R在線差異分析工具使用的代碼流程汉嗽。筆記里的代碼在理解每一步的基礎上谁榜，部分做了修改叉寂。但大致思路不變震糖，即使用limma包對芯片測序數(shù)據(jù)進行差異分析(暫未涉及可視化的代碼)女揭。
示例數(shù)據(jù)
- GSE：GSE43760
- Selected group design：
（1）Disease："GSM1070439","GSM1070441", "GSM1070443", "GSM1070445", "GSM1070447", "GSM1070449"
（2）Normal："GSM1070451", "GSM1070453", "GSM1070455", "GSM1070457", "GSM1070459", "GSM1070461"
分析流程Pipeline
Step 1 ：下載表達矩陣等
library(GEOquery)
GSE=paste0("GSE","43760")
gset=getGEO(GSE, getGPL = F)
class(gset)
#"list"
length(gset)
#1
- 如上有兩個注意點
（1）設置不下載GPL注釋信息蚤假，因為此前以專門整理了該GPL平臺的基因注釋信息，節(jié)省時間吧兔。詳見之前得到筆記http://www.reibang.com/p/d89a25d43549磷仰；
（2）初步下載的GSE對象格式之所以是list，是因為同一個GSE的samples可能分別在不同平臺（gpl）測序結(jié)果境蔼，對應list里的不同元素組成灶平。一般來說都是一個GSE對應一個GPL，但遇到多平臺情況箍土，需要注意自己選擇的樣本測序結(jié)果是什么平臺逢享，再進行選擇。
gset_sub=gset[[1]]
gset_sub
# ExpressionSet (storageMode: lockedEnvironment)
# assayData: 32968 features, 24 samples
# element names: exprs
# protocolData: none
# phenoData
# sampleNames: GSM1070439 GSM1070440 ... GSM1070462 (24 total)
# varLabels: title geo_accession ... treatment:ch1 (41 total)
# varMetadata: labelDescription
# featureData: none
# experimentData: use 'experimentData(object)'
# pubMedIds: 23771909
# Annotation: GPL6244
# Get GPL ID
gpl=unique(gset_sub@phenoData@data[,"platform_id"])
[1] "GPL6244"
#Get expression matrix
expr=gset_sub@assayData$exprs
#[1] 32968 24
#Get sample meta(optional)
meta=gset_sub@phenoData@data[,grep(":ch1",colnames(gset_sub@phenoData@data))]
#[1] "disease state:ch1" "gender:ch1" "individual:ch1"
#[4] "time:ch1" "tissue:ch1" "treatment:ch1"
Step 2：表達矩陣與分組設計
- 首先需要核對下吴藻，表達數(shù)據(jù)是否已經(jīng)經(jīng)過log轉(zhuǎn)換瞒爬；如果沒有則log轉(zhuǎn)換下。GEO2R通過下述代碼實現(xiàn)自動檢測沟堡。
qx <- as.numeric(quantile(expr, 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) { expr[which(expr <= 0)] <- NaN
expr <- log2(expr) }
- 然后根據(jù)選擇的分組取表達矩陣的子集
GSM_case=c("GSM1070439","GSM1070441", "GSM1070443", "GSM1070445", "GSM1070447", "GSM1070449")
GSM_normal=c("GSM1070451", "GSM1070453", "GSM1070455", "GSM1070457", "GSM1070459", "GSM1070461")
index_case=which(colnames(expr) %in% GSM_case)
index_control=which(colnames(expr) %in% GSM_normal)
expr_sle=expr[,c(index_case,index_control)]
- 分組設計
groups=factor(c(rep(1,length(index_case)),rep(0,length(index_control))),
levels=c(0,1),labels = c("normal","disease"))
Step 3：limma包差異分析
library(limma)
#首先設置分組矩陣侧但，指明所有sample的分組情況
design=model.matrix(~0+groups)
colnames(design) <- levels(groups)
# disease normal
# 1 0 1
# 2 0 1
# 3 0 1
# 4 0 1
# 5 0 1
# 6 0 1
# 7 1 0
# 8 1 0
# 9 1 0
# 10 1 0
# 11 1 0
# 12 1 0
#然后根據(jù)分組進行線性擬合
fit <- lmFit(expr_sle, design)
#需要進一步指明差異分析的比較為 disease組 VS normal組
cts <- paste("disease", "normal", sep="-") #注意順序：前者比上后者
cont.matrix <- makeContrasts(contrasts=cts, levels=design) #比較矩陣，指明哪兩組進行比較
fit2 <- contrasts.fit(fit, cont.matrix)
#最后得到差異分析結(jié)果
fit2 <- eBayes(fit2, 0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="p",number=Inf)
head(tT)
# logFC AveExpr t P.Value adj.P.Val B
# 8107769 -0.8016667 9.627500 -6.764851 1.062535e-05 0.1905344 2.556737
# 7892531 -0.6966667 7.201667 -6.207964 2.611828e-05 0.1905344 1.959251
# 7894928 -0.7716667 6.700833 -6.123500 3.004739e-05 0.1905344 1.863787
# 7893992 -0.9983333 7.504167 -6.065703 3.309032e-05 0.1905344 1.797715
# 8012349 -1.3166667 8.798333 -5.915814 4.258809e-05 0.1905344 1.623514
# 8059996 -0.9233333 7.895000 -5.890448 4.445985e-05 0.1905344 1.593624
deg <- data.frame(Probe=rownames(tT),P.Value=tT[,"P.Value"],
logFC=tT[,"logFC"],FDR=tT[,"adj.P.Val"])
Step 4：結(jié)合GPL平臺信息注釋symbo基因名
gpl_anno=read.csv(paste0("/your path/to/GPL_array/",gpl,".csv"))
gpl_anno[,1]=as.character(gpl_anno[,1])
deg=dplyr::left_join(deg,gpl_anno[,-3],by=c("Probe"="ID"))
deg=na.omit(deg)
rownames(deg)=1:nrow(deg)
# Probe logFC P.Value FDR Symbol
# 1 8107769 -0.8016667 -0.8016667 0.1905344 SLC12A2
# 2 8012349 -1.3166667 -1.3166667 0.1905344 PER1
# 3 8059996 -0.9233333 -0.9233333 0.1905344 PER2
# 4 7899018 -0.4683333 -0.4683333 0.1905344 TMEM57
# 5 7937696 0.7750000 0.7750000 0.2157829 KRTAP5-AS1
# 6 7897449 -0.4866667 -0.4866667 0.2164292 SPSB1
- 以上是針對GEO2R里結(jié)合limma對芯片數(shù)據(jù)分析的方法航罗；
- 如果是RNAseq的counts數(shù)據(jù)俊犯，則需要結(jié)合
edgeR包的voom()函數(shù)，總體思路不變（不需要考慮GPL的基因名注釋問題了）伤哺，相關(guān)代碼如下燕侠。 -
示例數(shù)據(jù)：GSE75852
#加載包
library(limma)
library(edgeR)
library(GEOquery)
library(data.table)
#下載數(shù)據(jù)
count_link="https://ftp.ncbi.nlm.nih.gov/geo/series/GSE75nnn/GSE75852/suppl/GSE75852_gene_count_frags.txt.gz"
# https://ftp.ncbi.nlm.nih.gov/geo/series/GSE112nnn/GSE112348/suppl/
counts=as.data.frame(fread(count_link))
colnames(counts)
counts[1:4,1:4]
#group分組
gset=getGEO("GSE75852",getGPL = F)
gsms=gset[[1]]@phenoData@data[,grep(":ch1",colnames(gset[[1]]@phenoData@data))]
index_case=which(gsms[,1]=="induced pluripotent stem cell" &
gsms[,2]=="Ataxia Telangiectasia")
index_control=which(gsms[,1]=="induced pluripotent stem cell" &
gsms[,2]=="Normal")
counts_sub=data.frame(counts[,c(index_case+1,index_control+1)],row.names=counts[,1])
groups=factor(c(rep(1,length(index_case)),rep(0,length(index_control))),
levels=c(0,1),labels = c("normal","disease"))
design=model.matrix(~0+groups)
colnames(design) <- levels(groups)
#limma+voom差異分析
v <- voom(counts_sub,design,normalize="quantile")
fit <- lmFit(v,design)
cts <- paste("disease", "normal", sep="-")
cont.matrix <- makeContrasts(contrasts=cts, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2, 0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="p",number=Inf)
deg <- data.frame(Probe=rownames(tT),P.Value=tT[,"P.Value"],
logFC=tT[,"logFC"],FDR=tT[,"adj.P.Val"],AveExpr=tT[,"AveExpr"])
head(deg)
