################某基因與多基因間相關(guān)性 表達(dá)
rm(list = ls())
setwd("/data3/ff/data/xena_TCGA/TPM/gene.lnc")
#####讀入表達(dá)文件? ? ? ? ? ? ? ? ? ?
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
#####在表達(dá)文件中提取所需基因
expr<-PAAD[gene,]
#轉(zhuǎn)置
expr<-t(expr)
expr<-as.data.frame(expr)
#刪除第一行
expr<-expr[-1,]
#使用Hmisc 包，計(jì)算矩陣相關(guān)系數(shù)及其對(duì)應(yīng)的顯著性水平
library(Hmisc)
#寫(xiě)函數(shù)來(lái)提取矩陣相關(guān)及其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 #查看函數(shù)說(shuō)明? ? 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)
#####################兩基因間相關(guān)性圖
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相關(guān)性
LA.HE1<-PAAD[c("LA","HE1"),]
#轉(zhuǎn)置
LA.HE1<-t(LA.HE1)
LA.HE1<-as.data.frame(LA.HE1)
#################################畫(huà)圖
library(ggplot2)
library(ggpmisc)
library(ggpubr)
#使用ggpubr包stat_cor，將相關(guān)系數(shù)和顯著性水平加入
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())+#移除網(wǎng)格線
? theme(plot.title = element_text(hjust = 0.5))+ #設(shè)置標(biāo)題居中
? 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",#將相關(guān)系數(shù)和顯著性水平加入
? ? ? ? ? #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
