################某基因与多基因间相关性 表达
rm(list = ls())
setwd("/data3/ff/data/xena_TCGA/TPM/gene.lnc")
#####读入表达文件
library(data.table)
PAAD<-fread("TCGA-PAAD.htseq-TPM.Gene.lnc.txt",data.table = F)
row.names(PAAD)<-PAAD$genenames
###########提取所需基因
gene<-read.table("/data3/ff/data/TF/AL5.TF.txt",header = T)
gene<-gene$TF
#####在表达文件中提取所需基因
expr<-PAAD[gene,]
#转置
expr<-t(expr)
expr<-as.data.frame(expr)
#删除第一行
expr<-expr[-1,]
#使用Hmisc 包，计算矩阵相关系数及其对应的显著性水平
library(Hmisc)
#写函数来提取矩阵相关及其P值
CorMatrix <- function(cor,p) {
ut <- upper.tri(cor)
data.frame(row = rownames(cor)[row(cor)[ut]] ,
column = rownames(cor)[col(cor)[ut]],
cor =(cor)[ut],
p = p[ut] )
}
res <- rcorr(as.matrix(expr), type=c("pearson"))
?rcorr #查看函数说明 rcorr(x, y, type=c("pearson","spearman"))
library(dplyr)
results <- CorMatrix (res$r, res$P) %>% as.data.frame() %>%
filter(row == "AL5" & cor <0)
write_csv(results,"AL5.negative.TF.csv",col_names = T)
#####################两基因间相关性图
rm(list = ls())
setwd("/data3/ff/data/xena_TCGA/TPM/gene.lnc")
PAAD<-fread("TCGA-PAAD.htseq-TPM.Gene.lnc.txt",data.table = F)
row.names(PAAD)<-PAAD$genenames
PAAD<-PAAD[,-1]
##############################################################---LA.HE1相关性
LA.HE1<-PAAD[c("LA","HE1"),]
#转置
LA.HE1<-t(LA.HE1)
LA.HE1<-as.data.frame(LA.HE1)
#################################画图
library(ggplot2)
library(ggpmisc)
library(ggpubr)
#使用ggpubr包stat_cor，将相关系数和显著性水平加入
my.formula <- y ~ x
ggplot(data=LA.HE1,aes(x=LA,y=HE1)) +
scale_x_log10() + scale_y_log10() +
labs(title="TCGA PAAD", x="LA expression", y="HE1 expression")+
theme_bw()+#移除背景
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+#移除网格线
theme(plot.title = element_text(hjust = 0.5))+ #设置标题居中
geom_point(size=2,alpha=1,color="black") +
#stat_poly_eq(formula = my.formula,
#aes(label = paste(#..eq.label..,#展示方程
#..rr.label..,#展示r值
# sep = "~~~")),
#parse = TRUE) +
stat_cor(method = "pearson",#将相关系数和显著性水平加入
#label.x =-2, label.y = -2,
color='black',p.accuracy = 0.001)+
geom_smooth(method="lm",color="red",alpha=.7,size=1.2,se=FALSE, formula = my.formula)+
scale_y_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10))+
scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10));
ggsave("LA.HE1.cor.pdf",width = 10, height =8,dpi = 300) #保存成pdf
