表达谱数据的整理与差异分析

点击蓝字 关注我们

最近，在《GEO数据库介绍与数据检索》和《如何从GEO下载数据辅助文章发表？》两篇教程中已经介绍过GEO数据库中表达谱数据的查找和下载，而在《要是早知道有这样的分析工具就好了！》一文中介绍了如何用GEO2R对表达谱数据进行差异分析，并提到了GEO2R所使用的R脚本（见R script选项窗口）。

那么，具体如何使用R脚本进行分析呢？下面就将这些脚本“掰碎”了，仍以《要是早知道有这样的分析工具就好了！》一文中的数据集（GSE111044）为例，详细为大家演示下如何使用R script 窗口的脚本进行分析吧。

一、数据下载与处理

#下载安装相关的R包；
#library(BiocManager)
#install("GEOquery")
#加载所需R包；
library(GEOquery)
library(limma)
library(umap)

#使用getGEO函数下载GEO表达量数据和平台注释数据；
gset <- getGEO("GSE111044", GSEMatrix =TRUE, AnnotGPL=TRUE)
#查看基因集对象中的数据集个数；
num <- length(gset)
num
#[1] 1

#获取gset列表中与平台数据匹配的数据集；
if (num > 1) {
idx <- grep("GPL570", attr(gset, "names"))
}else {
idx <- 1
}

#从list中提取匹配的数据集；
gset <- gset[[idx]]

#保存数据；
save(gset,file = "GSE111044_data.Rdata")
#载入数据；
load("GSE111044_data.Rdata")

#查看数据标签；
fvarLabels(gset)

#提取并调整变量标签便于用作矩阵列名，比如将"："、空格等转换成"."；
fvarLabels(gset) <- make.names(fvarLabels(gset))
fvarLabels(gset)

#使用featureData函数获取基因注释信息；
#主要包含varMetadata 和 data 两个文件，
fD <- featureData(gset)
#使用fData函数获取基因注释，即fd2中的data文件；
fd <- fData(gset)
#查看部分基因注释信息；
fd[1:8,1:3]

#使用pData函数获取样本信息；
pd <- pData(gset)
#预览样本信息；
pd[1:2]

#使用exprs函数提取基因表达量数据；
ex <- exprs(gset)
#预览数据；
head(ex)

#计算表达量数据的分位数；
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
qx
#[1] 2.083543 3.927326 5.468796 7.344907 11.905448 14.718111

#确定数据是否需要对数转换；
LogC <- (qx[5] > 100) || (qx[6]-qx[1] > 50 && qx[2] > 0)
LogC
#[1] FALSE

#对表达量数据进行log2转换（log2 transformation）；
#根据上面的结果，这里并没有进行log2转换；
if (LogC) { ex[which(ex <= 0)] <- NaN
exprs(gset) <- log2(ex) }

二、差异基因分析（基于limma）

#根据样本信息指定分组组合（1为Tumor组，0为Normal组）；
gsms <- "101010"
#分割成单个字符；
sml <- strsplit(gsms, split="")[[1]]
#[1] "1" "0" "1" "0" "1" "0"
#指定样本分组；
gs <- factor(sml)
gs
#[1] 1 0 1 0 1 0
#Levels: 0 1

#生成分组标签；
groups <- make.names(c("Normal","Tumor"))
groups
#[1] "Normal" "Tumor"
#设置分组水平；
levels(gs) <- groups
gs
#[1] Tumor Normal Tumor Normal Tumor Normal
#Levels: Normal Tumor
#在基因集对象中添加分组信息；
gset$group <- gs
#[1] Tumor Normal Tumor Normal Tumor Normal
#Levels: Normal Tumor

#生成模型设计矩阵；
design <- model.matrix(~group + 0, gset)
#替换列名；
colnames(design) <- levels(gs)
#线性模型拟合 （fit linear model）；
fit <- lmFit(gset, design)
#查看设计矩阵；
design

#生成比较标签；
cts <- c(paste(groups[1],"-",groups[2],sep=""))
#cts
#[1] "Normal-Tumor"

#构建比较矩阵；
cont.matrix <- makeContrasts(contrasts=cts, levels=design)
cont.matrix

#计算estimated coefficients和标准误;
fit2 <- contrasts.fit(fit, cont.matrix)

#差异分析(Empirical Bayes Statistics,经验贝叶斯)；
fit2 <- eBayes(fit2, 0.01)

#提取TOP200基因表格;
tT <- topTable(fit2, adjust="fdr", sort.by="p", number=200)
#选择相关的列；
tT <- subset(tT, select=c("ID","adj.P.Val","P.Value","logFC","Gene.symbol"))
head(tT)
#查看表格行列数；
dim(tT)

#提取所有基因表格；
tT2 <- topTable(fit2, adjust="fdr", sort.by="p", number=Inf)
#保存表格；
write.csv(tT2, file=paste0(cts,"_diff.csv"), row.names=F)

三、数据的可视化

#绘制直方图；
hist(tT2$adj.P.Val, col = "gold", border = "white", xlab = "P-adj",
ylab = "Number of genes", main = "P-adj value distribution")
#可见大多基因是不显著的；

#鉴定差异表达基因；
dT <- decideTests(fit2, adjust.method="fdr", p.value=0.05)
#1对应"up", -1对应"down" ，0对应 "not expressed"；
head(dT)

#参看比较组名称；
colnames(fit2)
#选择对照组；
ct <- 1
#绘制火山图；
volcanoplot(fit2, coef=ct, main=colnames(fit2)[ct], pch=20,
highlight=length(which(dT[,ct]!=0)), names=rep('+', nrow(fit2)))

#Mean-difference plot (MD-plot);
plotMD(fit2, column=ct, status=dT[,ct], legend=F, pch=20, cex=1)
abline(h=0,col="blue")

综上，GEO2R主要用GEOquery包下载数据，用limma R包进行差异分析的。好啦，本次的分享就到这里啦！

# SCIPainter

基迪奥旗下绘图公众号

分享科研绘图技能与工具

欢迎关注与转发~

你的好友拍了拍你

并请你帮她点一下“分享”~