1、目的
如果同学们筛选的是lncRNA的话可以用共表达为指标预测这个lncRNA对应的miRNA，最后构建出ceRNA网络。
我们的单基因也可以利用共表达找到与目标基因共调控或负调控的基因（其实这个也是一个发现新指标的路数）。
2、代码
remove(list = ls()) ##清空当前环境
setwd("F:/生信代码复现/相关性分析WQS") ##💚💚💚💚💚设置路径
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("limma")
#install.packages("ggplot2")
#install.packages("ggpubr")
#install.packages("ggExtra")
#加载所需的包
library(limma)
library(ggplot2)
library(ggpubr)
library(ggExtra)
library(future)
library(future.apply)
gene="ARL15" #💚💚💚💚💚目标基因的名称
corFilter=0.3 #💚💚💚💚💚筛选相关系数的阈值
pFilter=0.05 #💚💚💚💚💚筛选显著性检验的pvalue的阈值
expFile="F:/生信代码复现/相关性分析WQS/TPM.xlsx" #💚💚💚💚💚表达量文件的名称💙💙💙💙这个是最核心的文件，需要通过这个就是用来分析用的
#读取表达量文件，并将其转换为矩阵形式
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[rowMeans(data)>1,]
#删除样本分组信息
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=log2(data+1)
#获取目标基因的表达量
x=as.numeric(data[gene,])
#对每个基因进行循环，计算与目标基因的相关性和显著性检验
res3 <- data.frame()
genes <- colnames(rt)[-c(1:2)]
system.time({
outTab <- data.frame() # 定义一个空的数据框用于存储结果
res3 <- future_lapply(rownames(data), function(j) {
if(gene != j) { # 只有当基因不等于目标基因时才执行以下代码
y=as.numeric(data[j,])
corT=cor.test(x, y, method = 'pearson') #使用皮尔逊相关系数
cor=corT$estimate
pvalue=corT$p.value
outTab=rbind(outTab, cbind(Query=gene, Gene=j, cor, pvalue)) #将结果保存到一个表格中
#根据阈值筛选出显著相关的基因
if((abs(cor)>corFilter) & (pvalue<pFilter)){
#绘制散点图
df1=as.data.frame(cbind(x,y))
p1=ggplot(df1, aes(x, y)) +
xlab(paste0(gene, " expression"))+ ylab(paste0(j, " expression"))+
geom_point()+ geom_smooth(method="lm", formula=y~x) + theme_bw()+
stat_cor(method = 'pearson', aes(x =x, y =y))
pdf(file=paste0("cor.", j, ".pdf"), width=5, height=4.6) #保存为pdf文件
print(p1)
dev.off()
}
}
})
})
res3 <- data.frame(do.call(rbind,res3))
for(j in rownames(data)){
if(gene==j){next} #跳过目标基因本身
y=as.numeric(data[j,])
corT=cor.test(x, y, method = 'pearson') #使用皮尔逊相关系数
cor=corT$estimate
pvalue=corT$p.value
outTab=rbind(outTab, cbind(Query=gene, Gene=j, cor, pvalue)) #将结果保存到一个表格中
#根据阈值筛选出显著相关的基因
if((abs(cor)>corFilter) & (pvalue<pFilter)){
#绘制散点图
df1=as.data.frame(cbind(x,y))
p1=ggplot(df1, aes(x, y)) +
xlab(paste0(gene, " expression"))+ ylab(paste0(j, " expression"))+
geom_point()+ geom_smooth(method="lm", formula=y~x) + theme_bw()+
stat_cor(method = 'pearson', aes(x =x, y =y))
pdf(file=paste0("cor.", j, ".pdf"), width=5, height=4.6) #保存为pdf文件
print(p1)
dev.off()
}
}
#输出相关性结果的文件
write.table(file="corResult.txt", outTab, sep="\t", quote=F, row.names=F)
outTab=outTab[abs(as.numeric(outTab$cor))>corFilter & as.numeric(outTab$pvalue)<pFilter,]
write.table(file="corSig.txt", outTab, sep="\t", quote=F, row.names=F)