小云带你做单基因相关miRNA的分析绘图
基因-微小RNA(miRNA)的相关性分析,可以探索基因和miRNA之间的潜在相互作用和调控关系。意义在于帮助研究人员分析基因和miRNA之间的关系,并找到与目标基因相关的miRNA,以及它们之间的调控机制。这对于深入理解基因调控网络、发现潜在的生物标志物研究疾病的发生机制等方面都具有重要意义。
代码中选择了一个目标基因(gene)并对其与miRNA的表达量进行相关性分析。通过计算相关系数和p值,可以确定与目标基因表达量显著相关的miRNA。下面是原始下载的基因表达量数据和mirna的原始数据。
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("limma")
#install.packages("ggpubr")
#install.packages("ggExtra")
# 首先,检查是否已经安装了名为
"BiocManager"的包,如果没有则进行安装。 接下来,安#装"limma"、"ggpubr"和"ggExtra"包。 #引用包 library(limma) library(reshape2) library(ggpubr) library(ggExtra) corFilter=0.2 #相关系数过滤标准 pvalueFilter=0.001 #pvalue过滤标准 gene="VCAN" #示例展示的基因名称
expFile="diffGeneExp.txt" #基因表达文件
miExpFile="miRNAmatrix.txt" #miRNA表达文件
mirnaFile="miRNA.txt" #miRNA列表文件,之前rownames提取出的行名
#初始化了一些变量,包括相关系数过滤阈值(corFilter)、p值过滤阈值(pvalueFilter)、
#目标基因名称(gene)、基因表达文件名(expFile)、miRNA表达文件名(miExpFile)、#miRNA列表文件名(mirnaFile)
#读取基因表达文件,提取目标基因的表达量
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[gene,,drop=F]
#从指定的文件(expFile)中读取基因表达数据。基因名称作为行名,提取目标基因的表达
#量数据。将表达量数据转换为矩阵形式,存储在变量"data"中
#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
conData=data[,group==1,drop=F]
treatData=data[,group==0,drop=F]
logFC=mean(treatData)-mean(conData)
geneExp=data[,group==0,drop=F]
# 从基因表达数据中分离出正常样本和肿瘤样本。计算肿瘤样本和正常样本的平均表达
#量之差(logFC)。将肿瘤组的基因表达数据存储在变量"geneExp"中。
#读取miRNA表达数据
rt1=read.table(miExpFile, header=T, sep="\t", check.names=F)
rt1=as.matrix(rt1)
rownames(rt1)=rt1[,1]
exp1=rt1[,2:ncol(rt1)]
dimnames1=list(rownames(exp1),colnames(exp1))
miRNA=matrix(as.numeric(as.matrix(exp1)), nrow=nrow(exp1), dimnames=dimnames1)
miRNA=avereps(miRNA)
miRNA=miRNA[rowMeans(miRNA)>0,]
#从指定的文件(miExpFile)中读取miRNA表达数据。将miRNA名称作为行名,提取miRNA
#的表达量数据。将表达量数据转换为矩阵形式,存储在变量"miRNA"中
#提取目标miRNA表达量
miList=read.table(mirnaFile, header=F, sep="\t", check.names=F)
sameMirna=intersect(row.names(miRNA), as.vector(miList[,1]))
miRNA=miRNA[sameMirna,]
miRNA=log2(miRNA+1)
#正常和肿瘤样品数目
group=sapply(strsplit(colnames(miRNA),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
miRNAtumor=miRNA[,group==0]
conNum=length(group[group==1]) #正常组样品数目
treatNum=length(group[group==0]) #肿瘤组样品数目
sampleType=c(rep("Normal",conNum), rep("Tumor",treatNum))
colnames(miRNAtumor)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3-\\4",colnames(miRNAtumor))
#样品取交集
sameSample=intersect(colnames(miRNAtumor), colnames(geneExp))
miRNAtumor=miRNAtumor[,sameSample,drop=F]
geneExp=geneExp[,sameSample,drop=F]
y=as.numeric(geneExp[gene,])
# 从miRNA和基因表达数据中获取样本的类型(正常/肿瘤)信息。获取样本的类型和
#样本数目。确定miRNA和基因表达数据中样本的交集,即在两个数据中都存在的样本。
#相关性检验
outTab=data.frame()
for(i in row.names(miRNAtumor)){
if(sd(miRNAtumor[i,])>0.01){
#miRNA差异分析
miLogFC=mean(miRNA[i,((conNum+1):ncol(miRNA))])-mean(miRNA[i,1:conNum])
test=wilcox.test(miRNA[i,] ~ sampleType)
diffPvalue=test$p.value
#相关性分析
x=as.numeric(miRNAtumor[i,])
corT=cor.test(x, y, method="spearma")
cor=corT$estimate
pvalue=corT$p.value
outTab=rbind(outTab,cbind(Gene=gene, miRNA=i, cor, pvalue, logFC=miLogFC, diffPval=diffPvalue))
#对满足相关性过滤条件的miRNA进行可视化
if((cor< -corFilter) & (pvalue 对miRNA表达数据中的每个miRNA进行循环遍历:
如果该miRNA的表达量的标准差大于0.01:
计算miRNA在肿瘤组和正常组之间的平均表达量差异(miLogFC)。
进行miRNA表达量与目标基因表达量之间的相关性分析,使用的是Spearman相关系数和对应的p值。
将相关性结果存储在数据框"outTab"中。
对满足相关性过滤条件的miRNA进行可视化,包括绘制相关性散点图和差异的箱线图,并将结果保存为PDF文件。
#输出相关性结果
outTab=outTab[order(as.numeric(as.vector(outTab[,"cor"]))),]
write.table(file="net.network.txt",outTab,sep="\t",quote=F,row.names=F)
#输出相关性节点属性
miNode=data.frame(Node=unique(as.vector(outTab[,"miRNA"])), Type="miRNA")
geneNode=data.frame(Node=unique(as.vector(outTab[,"Gene"])), Type="Gene")
nodeOut=rbind(miNode, geneNode)
write.table(nodeOut, file="net.node.txt", sep="\t", quote=F, row.names=F)
#输出相关性显著的结果
outTab=outTab[((as.numeric(as.vector(outTab[,"cor"]))<-corFilter)&(as.numeric(as.vector(outTab[,"pvalue"])) < pvalueFilter)),]
write.table(file="cor.sig.txt",outTab,sep="\t",quote=F,row.names=F)
#输出相关性显著miRNA的表达量
corMirna=unique(as.vector(outTab[,"miRNA"]))
corMirnaExp=miRNA[corMirna,,drop=F]
corMirnaExp=rbind(ID=colnames(corMirnaExp), corMirnaExp)
write.table(corMirnaExp, file="cor.MirnaExp.txt", sep="\t", quote=F, col.names=F)
通过基因和miRNA之间的相关性分析,可以找到基因和miRNA之间的调控关系。通过分析它们之间的相关性,可以识别出与目标基因表达密切相关的miRNA。这有助于理解基因调控网络中的调节模式和机制。发现潜在的miRNA调控靶点:相关性分析可以帮助确定miRNA的潜在靶点基因。如果一个miRNA与某个基因呈现显著的相关性,那么该基因很可能是该miRNA的调控靶点。
下面我们来看看结果输出文件列表和结果图示范。
好了小云今天的讲解就到这里了,下期再见~