2024-06-17

数据探索之RMA预处理、limma差异分析，揭秘基因表达奥秘！热图与火山图带你直观感受差异表达魅力！

小果生信果 2024-04-02 19:00:09

大家好！我是小果！今天小果将为大家带来一场数据处理的盛宴！从数据的预处理开始，历经差异表达分析，最终得到结果可视化的精彩呈现。在这个过程中，我们探寻实验数据内在规律的同时，也能够挖掘出隐藏在基因差异表达背后的故事。跟随小果的脚步，让我们一同学习对实验数据进行差异表达分析及其可视化的方法吧！十年的生信之路，小果已经练就了一身扎实的本领，现在准备用这身本事为小伙伴们服务啦！如果你在生信分析上遇到了难题，那就来找小果吧！小果会用自己的专业知识和技能，为你解决困扰，助你一臂之力。期待你的联系哦~

公众号后台回复“111”

领取本篇代码、基因集或示例数据等文件

文件编号：240402

需要租赁服务器的小伙伴可以扫码添加小果，此外小果还提供生信分析，思路设计，文献复现等，有需要的小伙伴欢迎来撩~

(一)基本概念介绍

小果先和小伙伴们总述一下基本思路，首先利用RMA法预处理法对两组样本数据进行处理，然后对两组数据进行合并，利用limma R包对其进行差异分析，最后再对分析的结果进行可视化！

1.RMA算法

RMA算法基于20组探针的信号分布采用随机模型来评估表达值，对于低噪号的实验适用性更好一点。RMA算法不直接利用MM（Mismatch，不匹配）的信号去除背景噪音，而是通过特定的计算方式来评估和处理数据。

2.Limma R包

limma包是R语言中的一个广泛使用于差异表达分析的包，特别是处理来自微阵列和高通量测序的数据。limma包提供了丰富的统计方法，以识别在两组或多组之间表达水平显著不同的基因或转录本。本次的R包操作占用内存比较大，建议使用服务器，欢迎联系小果租赁性价比居高的服务器哦！

3.差异分析结果可视化

1)热图

主要通过颜色的深浅来反映基因在不同样本或条件下的表达水平。每个基因或样本在热图中占据一个单元格，颜色深浅代表表达量的大小。通过热图，我们可以直观地观察基因在不同样本间的表达差异，以及样本间的相似性。

2)火山图

侧重于展示基因表达差异的显著性。在火山图中，每个点代表一个基因，横坐标表示基因在不同样本或条件下的表达差异倍数（logFC），纵坐标表示差异表达的显著性水平（-log10 p值）。通过火山图，可以快速地筛选出表达差异显著的基因，并观察它们在整体基因表达谱中的分布情况。

(二)代码实操

废话不多说，小果直接带小伙伴们实操！

#RMA法预处理normal样本setwd("C:\Users\sfksfj\Desktop\GSE123924_RAW\N")  library(affyPLM)library(affy)Data<-ReadAffy()#**************sampleNames(Data)#***************N=length(Data)#用RMA预处理数据eset.rma<-rma(Data)#获取表达数据并输出到表格normal_exprs<-exprs(eset.rma)    probeid<-rownames(normal_exprs)normal_exprs<-cbind(probeid,normal_exprs)       write.table(normal_exprs,file="normal.expres.txt",sep='t',quote=F,row.names=F)

#RMA法预处理tumor样本setwd("C:\Users\sfksfj\Desktop\GSE123924_RAW\T")  #自己设定library(affyPLM)library(affy)Data<-ReadAffy()#**************sampleNames(Data)#**************N=length(Data)#用RMA预处理数据eset.rma<-rma(Data)#获取表达数据并输出到表格normal_exprs<-exprs(eset.rma)probeid<-rownames(normal_exprs)normal_exprs<-cbind(probeid,normal_exprs)        write.table(normal_exprs,file="tumor.expres.txt",sep='t',quote=F,row.names=F)

#合并N和T的数据setwd("C:\Users\sfksfj\Desktop\差异分析源文件\GSE123924_RAW\C")normal_exprs<-read.table("normal.expres.txt",header=T,sep="t")tumor_exprs<-read.table("tumor.expres.txt",header=T,sep="t")#T和N合并probe_exprs<-merge(normal_exprs,tumor_exprs,by="probeid")write.table(probe_exprs,file="cancer.probeid.exprs.txt",sep='t',quote=F,row.names=F)           #Probe ID转换为Gene symbol（先对平台文件进行整理）*************************setwd("C:\Users\sfksfj\Desktop\差异分析源文件\GSE123924_RAW\C")#读取基因表达文件probe_exp<-read.table("cancer.probeid.exprs.txt",header=T,sep="t",row.names=1)#读取平台探针文件probeid_geneid<-read.table("platform.txt",header=T,sep="t")probe_name<-rownames(probe_exp)#probe进行匹配loc<-match(probeid_geneid[,1],probe_name)#确定能匹配上的probe表达值probe_exp<-probe_exp[loc,]#每个probeid应对的geneidraw_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)#多个探针对应1个基因的情况，取平均值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="Gastric.cancer.geneid.exprs.txt",sep='t',quote=F,row.names=F)#将gene id 转换为gene symbolloc<-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="Gastric.cancer.genesyb.exprs.txt",sep='t',quote=F,row.names=F)           #补充缺失值(需要对Gastric.cancer.genesyb.exprs.txt进行处理，查找年号，删除无名基因)#最近邻居法（KNN，k-Nearest Neighbor）法：此方法是寻找和有缺失值的基因的表达谱相似的其他基因#通过这些基因的表达值（依照表达谱相似性加权）来填充缺失值setwd("C:\Users\sfksfj\Desktop\差异分析源文件\GSE123924_RAW\C")#调用函数imputelibrary(impute)#读取表达值gene_exp_matrix<-read.table("Gastric.cancer.genesyb.exprs.txt",header=T,sep="t",row.names=1)    gene_exp_matrix<-as.matrix(gene_exp_matrix)#KNN法计算缺失值imputed_gene_exp<-impute.knn(gene_exp_matrix,k=10,rowmax=0.5,colmax=0.8,maxp=3000,rng.seed=362436069)#读出经过缺失值处理的数据GeneExp<-imputed_gene_exp$data#写入表格genesymbol<-rownames(GeneExp)GeneExp<-cbind(genesymbol,GeneExp)write.table(GeneExp,file="Gastric.cancer.gene.exprs.txt",sep='t',quote=F,row.names=F)           #***********************（根据样本填写数量）setwd("C:\Users\sfksfj\Desktop\差异分析源文件\GSE123924_RAW\C")library(limma)rt<-read.table("Gastric.cancer.gene.exprs.txt",header=T,sep="t",row.names="genesymbol")class<-c(rep("normal",3),rep("tumor",3))#样本数量即为生成矩阵的CEL文件数量design<-model.matrix(~factor(class))colnames(design)<-c("normal","tumor")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)diffLab<-allDiff[with(allDiff, ((logFC>1 |logFC<(-1)) & adj.P.Val<0.05)),]write.table(diffLab,file="diffEXp.xls",sep="t",quote=F)#共表达diffExpLevel<-rt[rownames(diffLab),]qvalue=allDiff[rownames(diffLab),]$adj.P.ValdiffExpQvalue=cbind(qvalue,diffExpLevel)write.table(diffExpQvalue,file="diffExpLevel.xls",sep="t",quote=F)               #热图hmExp=log10(diffExpLevel+0.00001)library('gplots')hmMat=as.matrix(hmExp)pdf(file="heatmap.pdf",height=120,width=90)par(oma=c(3,3,3,5))heatmap.2(hmMat,col='greenred',trace="none",cexCol=1)        dev.off()

#火山图pdf(file="vol.pdf")xMax=max(-log10(allDiff$adj.P.Val))yMax=max(abs(allDiff$logFC))plot(-log10(allDiff$adj.P.Val),allDiff$logFC,xlab="adj.P.Val",ylab="logFC",main="Volcano",xlim=c(0,xMax),ylim=c(-yMax,yMax),pch=20,cex=0.4)diffSub=subset(allDiff,allDiff$adj.P.Val<0.05 & abs(allDiff$logFC)>1)#可修改，最低为0，logFC值越大差异倍数越高points(-log10(diffSub$adj.P.Val),diffSub$logFC,pch=20,col="red",cex=0.4)abline(h=0,lty=2,lwd=3)           dev.off()

(三)文章小结

今天，小果为大家详尽地展示了高通量基因表达数据（即microarray数据）从预处理到差异表达分析，再到结果可视化的整个流程。

首先，我们对normal样本和tumor样本数据进行了预处理，确保了数据的准确性和可靠性。接着，我们将normal和tumor样本数据进行了合并，以便在统一的框架下进行比较和分析，进而揭示不同样本间的基因表达差异。为了方便理解和解释，我们将Probe ID转换为了Gene symbol，使得每个基因的表达变化都能够直观地对应到具体的基因名称上。此外，针对数据中可能出现的缺失值，我们也进行了相应的处理，确保分析的完整性和准确性。紧接着，我们进行了差异表达分析，通过一系列统计和计算方法，识别出在normal和tumor样本间表达水平显著不同的基因。最后，为了更直观地展示这些差异表达基因的表达模式和显著性，我们生成了热图和火山图。热图以颜色的深浅展示了基因在不同样本间的表达差异，而火山图则清晰地标出了表达差异显著的基因，使我们迅速捕捉到关键信息。

整个流程不仅有助于我们深入理解基因表达的变化规律，还为后续的相关生物学实验和疾病研究提供了重要的参考依据。如果各位觉得自己运行代码太麻烦，欢迎用我们的云生信小工具，只要输入合适的数据就可以直接出想要的图呢，附云生信链接（http://www.biocloudservice.com/home.html）。希望本期内容能够对小伙伴们有所启发和帮助，让我们一起在生物信息学的海洋中探索更多未知的奥秘吧！

小果还提供思路设计、定制生信分析、文献思路复现；有需要的小伙伴欢迎直接扫码咨询小果，竭诚为您的科研助力！