GEO数据库表达谱差异基因分析

DoubleHelix

发布于 2019-08-07 11:05:31

8.7K6

发布于 2019-08-07 11:05:31

关于GEO数据库表达谱差异基因分析，网上有很多教程，但很多都不系统，几乎千篇一律，而且都是直接使用整理好的矩阵文件来操作的。大家都知道，GEO数据库只负责用户上传数据，而不负责对数据质量的控制，因此，有小伙伴也会发现，自己下载好的矩阵文件里面基因表达量数值特别大而且数据不集中，究其原因就是GEO数据库的数据参差不齐，不能确定上传者是否对整理好的数据进行了标准化处理。我们之前也讲过芯片数据的处理和分析流程，不了解的小伙伴们先读一下之前的文章：基因芯片数据挖掘分析表达差异基因。今天公众号：BioInfoCloud将从GEO芯片的原始数据进行分析，为大家详细的讲解。

我们选择了宫颈癌的表达芯片“GSE89657”来分析。

点击芯片的标题，就能看到芯片的全部信息了！包括实验设计，类型样本信息，作者，文献以及数据等等！直接看都能看懂！

将页面下拉至底部，第1个是矩阵文件(GEO分析最常用的)，第2个是原始文件(数据最精确的)。虽然说矩阵文件分析最简单，但是因为GEO不对芯片数据做质量控制，因此矩阵文件在某些时候并不是十分准确的。

下面开始下载数据了，首先我们需要下载原始文件，也就是格式为TAR(OF CEL)的文件，点击http下载原始文件是一个压缩包，同时也要下载平台文件。

下载平台文件是点击GPL6244进入的页面，底部有下载按钮。

下载后是下面这2个文件。

解压后的文件如下图，18个压缩包，对应18个样本的数据（第一个R文件是本例的代码，忽略）

再将所有的压缩包文件解压，讲原来的压缩包删除，或者移动到其他文件夹。这些cel文件就是原始数据了。

我们看GEO详情页里面的18个样本信息，有3个正常组织，其余都是肿瘤。

我们需要将文件进行分类，在工作目录建立一个cancer文件夹和一个normal文件夹，将相应的cel文件复制到相应文件夹中。注意，是复制，我们还要在当前文件夹里用所有的数据演示查看数据质量等操作。

接下来我们就开始分析演示

利用affy包处理原始数据

安装和加载affy包，如果已经安装，就直接加载！

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("affy")
library(affy)

其实，随着R版本的不同，加载该包时也需要很多基础包，需要先加载，而且每个人已经安装的包也不同和R版本的不同，这一过程可能会出错，反正在加载时出错，一般都是缺包或者需要加载一下包，缺什么补什么就行了！

library(Biobase)
library(affy)

读取数据，然后对数据集进行回归计算。

ReadAffy是一个用于读取数据的函数。允许用户读取MIAME信息和CEL文件的affybatch。如果在没有参数ReadAffy()的情况下调用该函数，那么将读取工作目录中的所有CEL文件并将其放入AffyBatch中。这些参数为用户提供了很大的灵活性。

setwd("F:\\BioInfoLab\\GEO_DATA")
cervical_carcinomas.Data<-ReadAffy()  
Pset <-fitPLM(cervical_carcinomas.Data)

绘制灰度图，查看数据质量，灰度图中颜色明显偏白的数据代表质量不好的数据。

image(cervical_carcinomas.Data[,1])

用权重图查看数据在整体中的重要程度，可以看出，本芯片数据较好，芯片质量较高。

image(Pset,type = "weights",which=1,main="weights")

残差图看数据的分布情况，主要是在回归分析中。

image(Pset,type = "resids",which=1,main="Residuals")

还有一种残差图，符号残差图，和残差图差不多的意义，只是图片在色彩上看着比较绚丽。

image(Pset,type = "sign.resids",which=1,main="Residuals.sign")

质量控制:相对对数表达（RLE），指一个探针组在某个样品的表达值除以该探针组在所有样品中表达之的中位数后取对数。反映平行实验的一致性。

colors <- brewer.pal(12,"Set3")
Mbox(Pset,col=colors,main = "RLE",las = 3)

质量控制：相对标准差（NUSE),指一个探针组在某个样品的PM值的标准差除以该探针组在各样品中的PM值标准差的中位数后取对数。反映平行实验的一致性。比RLE更为敏感。

boxplot(Pset,col=colors,main="NUSE",LAS=3)

质量控制：RNA降解图，它的原理是RNA降解从5’端开始，因为芯片结果5端荧光强度要远低于3’端。

data.deg <- AffyRNAdeg(cervical_carcinomas.Data)
plotAffyRNAdeg(data.deg,col = colors)
legend("topleft",sampleNames(cervical_carcinomas.Data),col = colors,lwd=1,inset = 0.05,cex=0.2)

分割线

接下来我们就开始处理数据和分析数据了！

对正常组数据进行预处理

#设置工作目录为正常组cel文件路径为工作目录
setwd("F:\\BioInfoLab\\GEO_DATA\\normal")
Data<-ReadAffy()  
sampleNames(Data) 
N <-length(Data)
eset.rma <- rma(Data)
normal_expre <- exprs(eset.rma)
probeid <- rownames(normal_expre)
normal_expre <- cbind(probeid,normal_expre)
write.table(normal_expre,file = "normal.expres.txt",sep = "\t",quote = F,row.names = F)

代码运行结束后会在工作目录产生一个normal.expres.txt文件。

一样的方法处理cancer组

######对cancer这样进行处理
setwd("F:\\BioInfoLab\\GEO_DATA\\cancer")
cancer_Data<-ReadAffy()  
sampleNames(cancer_Data) 
cancer_N <-length(cancer_Data)
cancer_eset.rma <- rma(cancer_Data)
cancer_expre <- exprs(cancer_eset.rma)
cancer_probeid <- rownames(cancer_expre )
cancer_expre <- cbind(cancer_probeid,cancer_expre)
write.table(cancer_expre,file = "cancer.expres.txt",sep = "\t",quote = F,row.names = F)

再新建一个文件夹命名为cel,将上述用RMA法处理的得到的两个txt文件放在cel文件下面。

然后将两组文件合并，得到cancer.probeid.expres.txt的文件。

setwd("F:\\BioInfoLab\\GEO_DATA\\cel")
normal_expres <- read.table("normal.expres.txt",header = T,sep = "\t")
cancer_expres <- read.table("normal.expres.txt",header = T,sep = "\t")
probe_expres <- merge(normal_expres,cancer_expres,by = "probeid")
write.table(probe_expres,file = "cancer.probeid.expres.txt",sep = "\t",quote=F,row.names = F)

将平台文件GPL6244-17930也放入cel文件夹里面。对平台文件与刚刚得到的标准化文件进行整合。

probe_exp <- read.table("cancer.probeid.expres.txt",header =T,sep="\t",row.names = 1)
# 读取探针文件
probeid_geneid <- read.table("GPL6244-17930.txt",header = T,sep="\t",blank.lines.skip = F)
write.csv(probeid_geneid,file = "probeid_geneid.csv",quote = T)
probe_name <- row.names(probe_exp)
# probe进行匹配
loc <- match(probeid_geneid[,1],probe_name)
# 确定能匹配上的probe表达值
probe_exp<- probe_exp[loc,]
#每个probeid应对的geneid
write.csv(probe_exp,file = "probe_exp.csv",quote = T)

raw_geneid<-as.numeric(as.matrix(probeid_geneid[,3]))
# 找出有geneid的probeid,并建立索引
index <- which(!is.na(raw_geneid))
# 提取geneid的probe
geneid <- raw_geneid[index]
# 找到每个geneid的表达值
exp_matrix <- probe_exp[index,]
geneidfactor <- factor(geneid)
#多个探针对应一个基因的情况，取均值
gene_exp_matrix <- apply(exp_matrix,2,function(x)tapply(x,geneidfactor,mean))
# geneid作为行名
rownames(gene_exp_matrix)<- levels(geneidfactor)
geneid <- rownames(gene_exp_matrix)
gene_exp_matrix2 <- cbind(geneid,gene_exp_matrix)
write.table(gene_exp_matrix2,file="cervical.cancer.geneid.exprs.txt",sep ="\t",quote = F,row.names = F)
loc <- match(rownames(gene_exp_matrix),probeid_geneid[,3])
rownames(gene_exp_matrix)=probeid_geneid[loc,2]
genesymbol <- rownames(gene_exp_matrix)
gene_exp_matrix3<-cbind(genesymbol,gene_exp_matrix)
write.table(gene_exp_matrix3,file="cervical.cancer.genesyb.exprs.txt",sep = "\t",quote = F,row.names = F)

对genesyb这个文件我们需要补充缺失值，这里采用KNN法，依照表达谱相似性加权来填充缺失值。

gene_exp_matrix <- read.table("cervical.cancer.genesyb.exprs.txt",header = T,sep="\t",row,names=1)
gene_exp_matrix<- as.matrix(gene_exp_matrix)
imputed_gene_exp <- impute.knn(gene_exp_matrix,k=10,rowmax = 0.5,colmax = 0.8,maxp=3000,rng.seed= 36242069)
GeneExp <- imputed_gene_exp$data
write.table(GeneExp,file = "cervical.cancer.exprs.txt",sep = "\t",quote = row.names= F)

通过以上方法，就可以整理出一个真正属于我们自己的矩阵文件，最后，对自己的矩阵文件求差异基因——使用R语言“limma”包。

rt <- read.table("cervical.cancer.exprs.txt",header = T,sep = "\t",row.names = "genesymbol")
class <- c(rep("normal",14),rep("cancer",4))
design<- model.matrix(~factor(class))
colnames(design)<- c("normal","cancer")
fit <- lmFit(rt.design)
fit2<- eBayes(fit)
allDiff<- topTable(fit2,adjust = "fdr",coef = 2,number = 200000)
write.table(allDiff,file = "limmaTab.xls",sep = "\t",quote = F)

可以看到，差异基因已经输出在cel文件下面了（limmaTab.xls）。