R计算mRNA和lncRNA之间的相关性+散点图

#差异表达的lncRNA
deLNC <- c('ENSG00000260920','ENSG00000242125','ENSG00000261211')
#差异表达的mRNA
dePC <- c('ENSG00000043355','ENSG00000109586','ENSG00000144355')
#所有的RNA名字
genes <- c(deLNC, dePC)
#所有的样本名字
samples <- c('TCGA-2F-A9KO-01', 'TCGA-2F-A9KP-01',
            'TCGA-2F-A9KQ-01', 'TCGA-2F-A9KR-01',
            'TCGA-2F-A9KT-01', 'TCGA-2F-A9KW-01')
#样本类型
type=factor(c("normal","disease","disease","normal","normal","disease"))

#表达谱数据
rnaExpr <- data.frame(matrix(c(2.7,7.0,4.9,6.9,4.6,2.5,
                    0.5,2.5,5.7,6.5,4.9,3.8,
                    2.1,2.9,5.9,5.7,4.5,3.5,
                    2.7,5.9,4.5,5.8,5.2,3.0,
                    2.5,2.2,5.3,4.4,4.4,2.9,
                    2.4,3.8,6.2,3.8,3.8,4.2),6,6),
                    stringsAsFactors=FALSE)
#设置表达谱的行名为RNA的名字
rownames(rnaExpr) <- genes
#设置表达谱的列名为样本的名字
colnames(rnaExpr) <- samples

#构建所有mRNA和lncRNA的组合
combination <- expand.grid(deLNC, dePC)
#第一列为lncRNA，第二列为mRNA
names(combination)=c("lnc","pc")

#通过循环来计算所有lncRNA和mRNA之间表达的相关性以及p值
cor_result=apply(combination,1,function(x){
  lnc=as.character(x[1])
  pc=as.character(x[2])
  result=cor.test(as.numeric(rnaExpr[lnc,]), as.numeric(rnaExpr[pc,]))
  score=c(pval=result$p.value,result$estimate)
  return(score)
})

#将lncRNA，mRNA的名字和相关性检验的结果合并起来
result=cbind(combination,t(cor_result))

#挑选p值小于0.05并且相关性大于0的lncRNA-mRNA对
index=which(result$pval<0.05 & result$cor>0)

#创建一个文件夹corplot来存放相关性图
dir.create("corplot")
#循环画出显著相关的mRNA和lncRNA的相关性散点图
for(i in index){
#获取p值
pval=round(result[i,3],4)
#获取相关系数
cor=round(result[i,4],2)
#获取lncRNA名字
lnc=as.character(result[i,1])
#获取lncRNA的表达量
lncExpr=as.numeric(rnaExpr[lnc,])
#获取mRNA名字
pc=as.character(result[i,2])
#获取mRNA的表达量 
pcExpr=as.numeric(rnaExpr[pc,])
#创建pdf文件
pdf(file=paste0("corplot/",lnc,"_",pc,".pdf"),width=10)
#绘制散点图
plot(lncExpr,pcExpr,
     col=c("red","blue")[type],
     pch=c(15,24)[type],
     xlab=paste(lnc," (lncRNA)"),
     ylab=paste(pc," (mRNA)"),     
     )
#添加拟合的直线
abline(lm(pcExpr~lncExpr),lty=2,col="green",lwd=2)
#添加图注
legend("bottomright",legend=levels(type),col=c("red","blue"),pch=c(15,24))
#在散点图上添加p值和相关系数
mtext(paste0("cor=",cor,"\npval=",pval), side=3,line= -2,adj = 0.1)
dev.off()
}