专栏首页生物信息云GEO数据库表达谱差异基因分析

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

关于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)。


本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-06-02

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 基于肿瘤微环境基因构建I-III期结肠癌免疫预后模型(IF=6.68)

    今天解读一篇发表在EBioMedicine上的文章,影响因子6.68。文章通过建立了一个肿瘤微环境风险评分panel,并证明其可作为I-III期结肠癌患者生存预...

    DoubleHelix
  • 肿瘤免疫细胞浸润与临床相关性分析

    Profiles of immune infiltration in colorectal cancer and theirclinical significa...

    DoubleHelix
  • R语言数据分析与挖掘(第八章):判别分析(3)——费歇尔(Fisher)判别分析

    我们之前介绍了判别分析中,因为判别准则的不同,可分为多种判别分析法。常用的有费歇尔(Fisher)判别分析、贝叶斯(Bayes)判别分析和距离判别分析。在上2篇...

    DoubleHelix
  • Linux 命令(89)—— less 命令

    版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。 ...

    Dabelv
  • 11 Confluent_Kafka权威指南 第十一章:流计算

    kafka 传统上被视为一个强大的消息总线,能够处理事件流,但是不具备对数据的处理和转换能力。kafka可靠的流处理能力,使其成为流处理系统的完美数据源,Apa...

    冬天里的懒猫
  • CAS+SSO配置单点登录完整案例

    CAS+SSO配置单点登录完整案例

    Java架构师必看
  • Android实现语音播放与录音功能

    本文实例为大家分享了Android实现语音播放与录音的具体代码,供大家参考,具体内容如下

    砸漏
  • 中诚信征信闫文涛:个人征信和企业征信未来将走向融合

    图丨中诚信征信总经理 闫文涛 “金融科技的价值在于利用大数据、人工智能等技术把风险识别出来,利用区块链技术对风险进行披露,再利用传统评级、传统金融模型把风险缓释...

    数据猿
  • Nature子刊:卒中的可塑性调控:一种新的神经功能恢复模型

    非侵入性脑刺激(NIBS)技术可以用来监测和调节皮层神经回路的兴奋性。长期的皮层刺激可对大脑功能产生持久影响,这为NIBS在慢性神经疾病中的治疗奠定了...

    用户1279583
  • Dubbo里面线程池的拒绝策略

    luozhiyun

扫码关注云+社区

领取腾讯云代金券