13.GEO數(shù)據(jù)集的R語言差異分析和代碼解析

忘仔扣扣糖 2022-04-26

展開全文

一、舉例回顧

本節(jié)下載GSE1009數(shù)據(jù)集，使用limma包進行差異分析舉例。

GSE1009

樣本量：共6個樣本，其中后3個為糖尿病腎病(DN)腎小球樣本，前3個為正常腎小球樣本。

使用芯片：Affymetrix Human Genome U95 Version 2 Array。

平臺：GPL8300。

二、需要準備的文件：

上一節(jié)中保存的GSE1009.Rdata和GSE1009expressionmetrix_GSE.csv（請把這兩個文件放在一個文件夾中。）

三、差異分析代碼解析：

>setwd("D:\\Rfile")

>rm(list = ls())

>options(stringsAsFactors=F)

#老規(guī)矩，先設置工作目錄。

>load(file = 'GSE1009.Rdata')

#加載上一節(jié)保存的R文件，里面包含表達矩陣和分組信息，不記得了可以看看上一節(jié)代碼。

>table(group_list)

#計數(shù)每組樣本的個數(shù)。

>library(limma)

#調用limma包，注意調用包之前，請先安裝。

#limma包是biocondutor中的包，搜索安裝命令的方式見上一節(jié)，這里就不重復了。最新的安裝命令如下：

>logFoldChange=1

>adjustP=0.05

#設置差異基因DEGs的篩選閾值。

#關于篩選閾值：一般有三個篩選條件，log2FC、P、FDR(調整的P),大家可以根據(jù)自己的需要設置。

>rt=read.csv("GSE1009expressionmetrix_GSE.csv")

#讀取整理好的表達矩陣數(shù)據(jù)，讀入R中的矩陣的數(shù)據(jù)結構為列表，將整個列表賦值給變量rt.

>rt=as.matrix(rt)

#將列表rt轉換為矩陣。

#這時表達矩陣如下：

注意行名

>rownames(rt)=rt[,1]

#將第一列數(shù)據(jù)（基因名）作為行名。

再看看行名，已經變了

>exp=rt[,2:ncol(rt)]

#選擇第2列到最后一列，即表達矩陣，并賦值給exp。這樣就得到行名為基因名，列名為樣本名的矩陣。

>dimnames=list(rownames(exp),colnames(exp))

#提取出exp的行名和列名。這里的dimnames就是矩陣創(chuàng)建matrix里面的參數(shù)，在數(shù)據(jù)結構矩陣那一節(jié)有講。

>rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)

#合并整個矩陣，得到一個行名為基因名，列名為樣本名的表達矩陣。

>modType=c(rep("Control",3),rep("DN",3))

#構建分組模型，根據(jù)GSE1009樣本的排列，前3個為正常樣本，后3個為糖腎樣本，構建一個分組的向量模型。

>design <- model.matrix(~0+factor(modType))

>colnames(design) <- c("Control","DN")

#創(chuàng)建線性模型。

#注意這里，在構建模型時，有一些數(shù)據(jù)集的處理組在前，對照組在后，且因為R里面默認的字符存儲順序為字母表順序，如果處理組的首寫字母的順序在對照組的后面，構建模型時前三個又是處理組，那design出來的結果會上下調就完全相反。所以如果是這種情況的話，需要定義因子順序。如：“modType=c(rep("DN",40),rep("Control",40))

design <- model.matrix(~0+factor(modType,levels = c("DN","Control"),ordered = F)).”

>fit <- lmFit(rt,design)

>cont.matrix<-makeContrasts(DN-Control,levels=design)

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

>fit2 <- eBayes(fit2)

>allDiff=topTable(fit2,adjust='fdr',number=200000)

#貝葉斯檢驗。

>diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]

#根據(jù)篩選條件篩選出差異基因。

>diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]

#篩選上調的基因。

>diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]

#篩選下調的基因。

##如果想把這些基因的結果輸出，可以用write.table函數(shù)write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F)，要基因名字就把參數(shù)row.names設置為T.

##繪制熱圖

>group<-modType

#設置熱圖的分組。

>group<-as.data.frame(group)

>rownames(group)<-colnames(rt)

#按照分組，將rt的列名改為DN或control

>heat<-rt[rownames(rt) %in% c(head(rownames(subset(diffSig,diffSig$logFC>0)),20),head(rownames(subset(diffSig,diffSig$logFC<0)),20)),]

>library(pheatmap)

>x <- t(scale(t(heat)))

>pheatmap(x,annotation_col=group)

#繪制前20個基因的熱圖，自己調整基因數(shù)、圖形顏色等等。

這里GSE1009的熱圖就繪制成功，實際處理時大家可根據(jù)自己的需求調整圖形參數(shù)。

##繪制火山圖

>png(file="火山圖.png",width = 600,height = 600)

#設置圖片格式

>yMax=max(-log10(allDiff$adj.P.Val))

#定義y軸最大值

>xMax=max(abs(allDiff$logFC))

#定義x軸最大值

>plot(allDiff$logFC,-log10(allDiff$adj.P.Val), xlab="logFC",ylab="-log10(adj.P.Val)",main="Volcano", xlim=c(-xMax,xMax),ylim=c(0,yMax),yaxs="i",pch=19, cex=1.2)

#繪圖

>diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange)

>points( diffSub$logFC,-log10(diffSub$adj.P.Val), pch=19, col="red",cex=1.2)

#修改上調基因的圖形參數(shù)

>diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange))

>points(diffSub$logFC,-log10(diffSub$adj.P.Val), pch=19, col="green",cex=1.2)

#修改下調基因的圖形參數(shù)

>abline(h=-log10(adjustP),lty=2,lwd=2)

>abline(v=c(-1,1),lty=2,lwd=2)

#拼接圖形

>dev.off()

#返回終端。