ceRNAs network analysis using GDCRNATools
GDCRNATools: integrative analysis of protein coding genes, long non-coding genes, and microRNAs in GDC
github
library(GDCRNATools)
####### ceRNA network analysis #######
deLNC <- c('ENSG00000260920','ENSG00000242125','ENSG00000261211')
dePC <- c('ENSG00000043355','ENSG00000109586','ENSG00000144355')
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')
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)
rownames(rnaExpr) <- genes
colnames(rnaExpr) <- samples
mirExpr <- data.frame(matrix(c(7.7,7.4,7.9,8.9,8.6,9.5,
5.1,4.4,5.5,8.5,4.4,3.5,
4.9,5.5,6.9,6.1,5.5,4.1,
12.4,13.5,15.1,15.4,13.0,12.8,
2.5,2.2,5.3,4.4,4.4,2.9,
2.4,2.7,6.2,1.5,4.4,4.2),6,6),
stringsAsFactors=FALSE)
colnames(mirExpr) <- samples
rownames(mirExpr) <- c('hsa-miR-340-5p','hsa-miR-181b-5p',
'hsa-miR-181a-5p', 'hsa-miR-181c-5p',
'hsa-miR-199b-5p','hsa-miR-182-5p')
rnaExpr[1:3,1:3]
mirExpr[1:3,1:3]
图片alt
ceOutput <- gdcCEAnalysis(lnc = deLNC,
pc = dePC,
lnc.targets = 'starBase',
pc.targets = 'starBase',
rna.expr = rnaExpr,
mir.expr = mirExpr)
ceOutput
图片alt
# ceOutput2 <- ceOutput[ceOutput$hyperPValue<0.01 &
# ceOutput$corPValue<0.01 & ceOutput$regSim != 0,]
edges <- gdcExportNetwork(ceNetwork = ceOutput, net = 'edges')
nodes <- gdcExportNetwork(ceNetwork = ceOutput, net = 'nodes')
edges
nodes
图片alt
library(ggplot2)
library(ggridges)
head(diamonds)
ggplot(diamonds, aes(x = price, y = cut,fill=cut,
group=cut, height=..density..)) +
geom_density_ridges(scale = 4) +
scale_y_discrete(expand = c(0, 0)) + # will generally have to set the `expand` option
scale_x_continuous(expand = c(0, 0)) + # for both axes to remove unneeded padding
coord_cartesian(clip = "off") + # to avoid clipping of the very top of the top ridgeline
theme_ridges()
Venn
png(
file = "figure/ggvenn.png",
width = 15, height = 15, units = "cm", bg = "white", res = 300
)
ggvenn(
list("TCGA Mutation" = TCGA_mutation, "Sequencing Mutation" = my_mutation),
stroke_size = 0.5, set_name_size = 5, text_size = 5
)
dev.off()
################################################################################
#
# 食管癌GEO ceRNA
#
################################################################################
{
library("optparse")
library(tidyverse)
library(tidyverse)
library(ggrepel)
library(edgeR)
library(pheatmap)
library('SpidermiR')
library("ggvenn")
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggpubr)
library(lattice)
library(survivalROC)
option_list = list(make_option(c("-d","--directory"), type = "character", help = "test"))
args <- parse_args(OptionParser(option_list=option_list))
if(!is.null(args$directory)){
if(!file.exists(args$directory)){
message("请确保 ",args$directory," 存在!")
q()
}
setwd(args$directory)
message("切换工作目录到:",args$directory)
}
if(!file.exists("figures")){
dir.create(file.path("figures"))
}
if(!file.exists("figures/GEO_Volcano_heatmap")){
dir.create(file.path("figures/GEO_Volcano_heatmap"))
}
if(!file.exists("figures/venn")){
dir.create(file.path("figures/venn"))
}
if(!file.exists("figures/survival")){
dir.create(file.path("figures/survival"))
}
if(!file.exists("figures/cox")){
dir.create(file.path("figures/cox"))
}
if(!file.exists("result/enrichment")){
dir.create(file.path("result/enrichment"))
}
}
{
IS_DEBUG=T
LOG_FILENAME= "report_download.txt"
IS_EMPTY=T
LOG_FILENAME= source("/home/wangyang/workspace/bioinfo_analysis/Rscript/tools/tools.R")
}
####################################################
## 公共函数
###################################################
{
getDeg <- function(deg){
res <- deg %>% rownames_to_column("symbol") %>%
dplyr::select(symbol,logFC=logFC,Pvalue=adj.P.Val)
return(res)
}
Pheatmap <- function(deg,expr,filename){
deg_gene <- getDeg(deg)%>%
filter(abs(logFC)>2,Pvalue<0.01)%>%
arrange(desc(abs(logFC)))%>%
pull("symbol")
deg_expr <- expr[deg_gene,]%>%
{log10(cpm(.)+0.01)}
#deg_expr <- t(scale(t(deg_expr)))
#deg_expr[deg_expr > 2] = 2
#deg_expr[deg_expr < -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,show_rownames = F,
filename = paste0("figures/",filename,".png"))
return(paste0("figures/",filename,".png"))
}
Volcano <-function(deg,filename,title){
# 火山图----------------------
getDeg(deg)%>%
mutate(direction = factor(ifelse(Pvalue < 0.01 & abs(logFC)>2,
ifelse(logFC>1,"Up","Down"),
"NS"),levels = c("Up","Down","NS"))) %>%
na.omit()%>%
ggplot(aes(x=logFC,y=-log10(Pvalue),colour=direction))+
geom_point(alpha=0.6) +
scale_color_manual(values = c("red","green","#000000"))+
xlab(expression(log[2]("Fold Change")))+
ylab(expression(-log[10]("P Value")))+
theme_bw()+
theme(legend.title = element_blank(),
plot.title = element_text(hjust = 0.5,size = 24, face= "bold"),
axis.title = element_text(size=20),
legend.text= element_text(size=20),
aspect.ratio=1)+
geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6) +
labs(title = title)+
geom_text_repel(data=. %>%filter(Pvalue<0.001,abs(logFC)>2),
aes(label=symbol),
size=3,
segment.color="black",
show.legend = F)+
ggsave(filename = paste0("figures/",filename,".png"),width = 20,height = 15,units="cm")
return(paste0("figures/",filename,".png"))
}
get_mircode <- function(){
if(!file.exists("result/mircode.rda")){
mircode <- read.table("/home/wangyang/workspace/bioinfo_analysis/Rscript/data/mircode_highconsfamilies.txt",sep = "\t",header = T)
table(mircode$gene_class)
lncmiRcode = mircode[mircode$gene_class %in%c("lncRNA, intergenic","lncRNA, overlapping"),1:4]
head(lncmiRcode$microrna)
library(stringr)
p1 = str_starts(lncmiRcode$microrna,"miR-")
table(p1)
p2 = str_starts(lncmiRcode$microrna,"let-")
table(p2)
lncmiRcode_let <- lncmiRcode[p2,]%>%
dplyr::select(c("gene_symbol","microrna"))%>%
mutate(miRNA=str_replace_all(microrna,"let-",""))%>%
separate_rows(miRNA,sep = "/",convert = T)%>%
mutate_at(vars(contains("miRNA")), ~ paste0("hsa-let-",.))%>%
mutate(miRNA = str_extract(miRNA,"hsa-let-[0-9]+[a-z]?"))%>%
dplyr::select(c("gene_symbol","miRNA"))
lncmiRcode_miR <- lncmiRcode[p1,]%>%
dplyr::select(c("gene_symbol","microrna"))%>%
mutate(miRNA=str_replace_all(microrna,"miR-",""))%>%
separate_rows(miRNA,sep = "/",convert = T)%>%
mutate_at(vars(contains("miRNA")), ~ paste0("hsa-mir-",.))%>%
mutate(miRNA = str_extract(miRNA,"hsa-mir-[0-9]+[a-z]?"))%>%
dplyr::select(c("gene_symbol","miRNA"))
lncmiRcode_f <- rbind(lncmiRcode_let,lncmiRcode_miR)
saveRDS(lncmiRcode_f,file = "result/mircode.rda")
return(lncmiRcode_f)
}
return(readRDS("result/mircode.rda"))
}
Volcano_sig <-function(deg,title=NULL,filename=NULL,gene=NULL){
# 火山图----------------------
res <- getDeg(deg)%>%
mutate(direction = factor(ifelse(Pvalue < 0.01 & abs(logFC)>2,
ifelse(logFC>1,"Up","Down"),
"NS"),levels = c("Up","Down","NS"))) %>%
na.omit()%>%
ggplot(aes(x=logFC,y=-log10(Pvalue),colour=direction))+
geom_point(alpha=0.6) +
scale_color_manual(values = c("red","green","#000000"))+
xlab(expression(log[2]("Fold Change")))+
ylab(expression(-log[10]("P Value")))+
theme_bw()+
theme(legend.title = element_blank(),
plot.title = element_text(hjust = 0.5,size = 24, face= "bold"),
axis.title = element_text(size=20),
legend.text= element_text(size=20),
aspect.ratio=1)+
geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6)
if(!is.null(title)){
res <- res+labs(title = title)
}
if(!is.null(gene)){
res <- res + geom_text_repel(data=. %>%filter(symbol %in% gene),
aes(label=symbol),
size=3,
segment.color="black",
show.legend = F)
}
if(!is.null(filename)){
res <- res + ggsave(filename = paste0("figures/",filename,".png"),width = 20,height = 15,units="cm")
}
return(res)
}
basesave <- function(expr, filename, print=T) {
#extension
exten = stringr::str_match(filename, "\\.(\\w+)$")[, 2]
switch(exten,
png = {
png(filename)
eval(expr, envir = parent.frame())
dev.off()
},
{stop("filetype not recognized")})
#print?
if (print) eval(expr, envir = parent.frame())
invisible(NULL)
}
}
#####################################################
## 热图火山图
#####################################################
(function(){
deg <- readRDS("result/limma_mRNA_deg.rds")
Volcano(deg,"GEO_Volcano_heatmap/1_mRNA_Volcano",title = "mRNA Volcano")
expr <- readRDS("result/gset_expr_annot_mRNA_GSE89102.rda")
Pheatmap(deg,expr,filename="GEO_Volcano_heatmap/1_mRNA_pheatmap")
deg <- readRDS("result/limma_lnRNA_deg.rds")
Volcano(deg,"GEO_Volcano_heatmap/1_lnRNA_Volcano",title = "lnRNA Volcano")
expr <- readRDS("result/gset_expr_annot_lnRNA_GSE89102.rda")
Pheatmap(deg,expr,filename="GEO_Volcano_heatmap/1_lnRNA_pheatmap")
deg <- readRDS("result/limma_miRNA_deg_2.rds")
Volcano(deg,"GEO_Volcano_heatmap/1_miRNA_Volcano",title = "miRNA Volcano")
if(F){
res <- getDeg(deg)%>%
mutate(direction = factor(ifelse(Pvalue < 0.05 & abs(logFC)>2,
ifelse(logFC>1,"Up","Down"),
"NS"),levels = c("Up","Down","NS"))) %>%
na.omit()%>%
ggplot(aes(x=logFC,y=-log10(Pvalue),colour=direction))+
geom_point(alpha=0.6) +
scale_color_manual(values = c("red","green","#000000"))+
xlab(expression(log[2]("Fold Change")))+
ylab(expression(-log[10]("P Value")))+
theme_bw()+
theme(legend.title = element_blank(),
plot.title = element_text(hjust = 0.5,size = 24, face= "bold"),
axis.title = element_text(size=20),
legend.text= element_text(size=20),
aspect.ratio=1)+
geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)
print(res)
}
expr <- readRDS("result/gset_expr_annot_miRNA_GSE59973.rda")
Pheatmap(deg,expr,filename="GEO_Volcano_heatmap/1_miRNA_pheatmap")
if(F){
deg_gene <- getDeg(deg)%>%
filter(abs(logFC)>2,Pvalue<0.05)%>%
arrange(desc(abs(logFC)))%>%
pull("symbol")
deg_expr <- expr[deg_gene,]
#deg_expr <- t(scale(t(deg_expr)))
#deg_expr[deg_expr > 2] = 2
#deg_expr[deg_expr < -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,show_rownames = F)
# return(paste0("figures/",filename,".png"))
}
})()
if(F){
deg <- readRDS("result/miRNA_iso_expr_DESeq.rda")
Volcano(deg,"1_miRNA_iso_Volcano",title = "miRNA Volcano")
expr <- readRDS("result/miRNA_iso_expr.rda")
Pheatmap(deg,expr,filename="1_miRNA_iso_pheatmap")
deg_gene <- getDeg(deg)%>%
filter(abs(logFC)>1,Pvalue<0.05)%>%
arrange(desc(abs(logFC)))%>%
pull("symbol")
deg_expr <- expr[deg_gene,]%>%
{log10(cpm(.)+0.01)}
#deg_expr <- t(scale(t(deg_expr)))
#deg_expr[deg_expr > 2] = 2
#deg_expr[deg_expr < -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,show_rownames = F,
filename = "figures//1_miRNA_iso_pheatmap.png")
}
#############################################
## 差异基因提取 变量
#############################################
miRNA_deg_sig <- getDeg(readRDS("result/limma_miRNA_deg_2.rda"))%>%
filter(abs(logFC)>2,Pvalue<0.01)%>%
mutate(symbol=sub("R","r",symbol))%>%
mutate(symbol= case_when(grepl("hsa-mir",symbol)~str_extract(symbol,"hsa-mir-[0-9]+[a-z]?"),
grepl("hsa-let",symbol)~str_extract(symbol,"hsa-let-[0-9]+[a-z]?")))%>%
arrange(Pvalue)
report("以abs(logFC)>2 & Pvalue<0.01 为条件筛选差异miRNA共 ",length(unique(miRNA_deg_sig$symbol))," 个")
lncRNA_deg_sig <- getDeg(readRDS("result/limma_lnRNA_deg.rda"))%>%
filter(abs(logFC)>2,Pvalue<0.01)
cat("以abs(logFC)>2 & Pvalue<0.01 为条件筛选差异lncRNA共 ",length(unique(lncRNA_deg_sig$symbol))," 个")
mRNA_deg_sig <- getDeg(readRDS("result/limma_mRNA_deg.rda"))%>%
filter(abs(logFC)>2,Pvalue<0.01)
cat("以abs(logFC)>2 & Pvalue<0.01 为条件筛选差异mRNA共 ",length(unique(mRNA_deg_sig$symbol))," 个")
lnRNA_expr <- readRDS("result/lncRNA_expr.rda")%>%
dplyr::select(!matches("TCGA-([A-Z 0-9]*?)-([A-Z 0-9]*?)-(11.*?)-(.*?)"))%>%
rownames_to_column("symbol")%>%
mutate(symbol=str_extract(symbol,"[A-Z 0-9 -]+"))
lnRNA_expr <- lnRNA_expr[!duplicated(lnRNA_expr$symbol),]%>%
rownames_to_column("id")%>%
select(-id)%>%
column_to_rownames("symbol")
sum(rownames(lnRNA_expr) %in% lncRNA_deg_sig$symbol)
lnRNA_expr <- lnRNA_expr[rownames(lnRNA_expr) %in% lncRNA_deg_sig$symbol,]
mRNA_expr <- readRDS("result/mRNA_expr.rda")%>%
dplyr::select(!matches("TCGA-([A-Z 0-9]*?)-([A-Z 0-9]*?)-(11.*?)-(.*?)"))
miRNA_expr <- readRDS("result/miRNA_expr.rda")%>%
dplyr::select(!matches("TCGA-([A-Z 0-9]*?)-([A-Z 0-9]*?)-(11.*?)-(.*?)"))%>%
rownames_to_column("symbol")%>%
mutate(symbol= case_when(grepl("hsa-mir",symbol)~str_extract(symbol,"hsa-mir-[0-9]+[a-z]?"),
grepl("hsa-let",symbol)~str_extract(symbol,"hsa-let-[0-9]+[a-z]?")))%>%
aggregate(.~symbol,.,mean)%>%
column_to_rownames("symbol")
#############################################
## lnRNA ralated miRNA
#############################################
## 获得miRNA与lnRNA对应关系
(function(){
cat("\n***************************lnRNA ralated miRNA***************************")
if(T){
mircode <- read_tsv("/home/wangyang/workspace/pipeline/CANCERjupyter/Common/lncRNA_miRCODE_miRNA.csv")%>%
plyr::rename(c(lncRNA="gene_symbol"))
}else{
mircode <- get_mircode()
}
cat("从miRcode数据库获得lnRNA-miRNA的对应关系：共有 ",dim(mircode)[1]," 个")
lnRNAPrediction <- mircode%>%
filter(gene_symbol %in% rownames(lnRNA_expr))
cat("在miRcode中找到的lnRNA共有 ",length(unique(lnRNAPrediction$gene_symbol))," 个，","对应的miRNA有 ",length(unique(lnRNAPrediction$miRNA))," 个")
miRNA_lnRNA <- lnRNAPrediction%>%
filter(miRNA %in% miRNA_deg_sig$symbol)%>%
plyr::rename(c(gene_symbol="lnRNA"))%>%
dplyr::select(-gene_id)
cat("预测得到的 ",length(unique(lnRNAPrediction$miRNA))," 个miRNA与差异的miRNA取交集, 得到 ",length(unique(miRNA_lnRNA$miRNA)), "个miRNA")
miRNA_lnRNA <- miRNA_lnRNA[!duplicated(miRNA_lnRNA),]
saveRDS(miRNA_lnRNA,file = "result/miRNA_lnRNA.rda")
miRNA_intersect <- unique(miRNA_lnRNA$miRNA)
lnRNA_intersect <- unique(miRNA_lnRNA$lnRNA)
#length(miRNA_intersect)
#length(lnRNA_intersect)
cat("最后得到 ",length(lnRNA_intersect)," 个lncRNA可以与 ",length(miRNA_intersect)," 个miRNA相互作用")
#png(filename="plot_%02d.png")
png(file="figures/venn/ggvenn_lnRNA_prediction_miRNA.png",
width = 15,height = 15,units = "cm",bg = "white",res=300)
res <- ggvenn(
list("lncRNA-ralated miRNA"=lnRNAPrediction$miRNA,"DE miRNAs"=miRNA_deg_sig$symbol),
stroke_size = 0.5, set_name_size = 5,text_size = 5
)
print(res)
dev.off()
cat("lnRNA预测miRNA与差异miRNA Venn 图见figures/ggvenn_lnRNA_prediction_miRNA.png")
cat("***************************lnRNA ralated miRNA***************************")
})()
#############################################
## miRNA ralated mRNA
#############################################
(function(){
cat("\n***************************miRNA ralated mRNA***************************")
miRNA_lnRNA <- readRDS("result/miRNA_lnRNA.rda")
miRNA_intersect <- unique(miRNA_lnRNA$miRNA)
## test start hsa-let-7a-2 预测时去掉 -2
if(F){
miRNAprediction <- SpidermiRdownload_miRNAprediction(mirna_list=c("hsa-let-7a"))
}
## test end
if(!file.exists("result/miRNAprediction.rda")){
miRNAprediction <- SpidermiRdownload_miRNAprediction(mirna_list=miRNA_intersect)
save(miRNAprediction,file = "result/miRNAprediction.rda")
}else{
load("result/miRNAprediction.rda")
}
cat("用 ",length(miRNA_intersect)," 个miRMA使用[R package SpidermiR]预测得到 ",length(miRNAprediction$V2)," 个mRNA")
miRNA_mRNA <- miRNAprediction%>%
plyr::rename(c(V1="miRNA",V2="mRNA"))%>%
filter(mRNA %in% mRNA_deg_sig$symbol)
miRNA_mRNA <- miRNA_mRNA[!duplicated(miRNA_mRNA),]
saveRDS(miRNA_mRNA,file = "result/miRNA_mRNA.rda")
cat("用预测得到mRNA与差异的miRNA去交集, 得到 ",length(unique(miRNA_mRNA$mRNA)), "个mRNA")
miRNA_intersect <- unique(miRNA_mRNA$miRNA)
mRNA_intersect <- unique(miRNA_mRNA$mRNA)
#length(miRNA_intersect)
#length(mRNA_intersect)
cat("最后得到 ",length(miRNA_intersect)," 个miRNA可以靶向 ",length(mRNA_intersect)," 个mRNA相互作用")
library("ggvenn")
png(file="figures/venn/ggvenn_miRNA_predication_mRNA.png",
width = 15,height = 15,units = "cm",bg = "white",res=300)
res <- ggvenn(
list("miRNA-ralated mRNA"=miRNAprediction$V2,"DE mRNAs"=mRNA_deg_sig$symbol),
stroke_size = 0.5, set_name_size = 5,text_size = 5
)
print(res)
dev.off()
cat("miRNA预测mRNA与差异miRNA Venn 图见figures/ggvenn_miRNA_predication_mRNA.png")
cat("***************************lnRNA ralated miRNA***************************")
})()
if(F){
list<-SpidermiRdownload_miRNAvalidate(validated)
list%>%
filter(V1 %in% miRNA_intersect)
org<-SpidermiRquery_species(species)
net_type<-SpidermiRquery_networks_type(organismID=org[6,])
net_shar_prot<-SpidermiRquery_spec_networks(organismID = org[6,]
,network = "pred")
}
#############################################
### 构建ceRNA网络
#############################################
if(F){
getDeg(readRDS("result/miRNA_expr_DESeq.rda"))%>%
mutate(direction = factor(ifelse(Pvalue < 0.01 & abs(logFC)>2,
ifelse(logFC>1,"Up","Down"),
"NS"),levels = c("Up","Down","NS"))) %>%
na.omit()%>%
ggplot(aes(x=logFC,y=-log10(Pvalue),colour=direction))+
geom_point(alpha=0.6) +
scale_color_manual(values = c("red","green","#000000"))+
xlab(expression(log[2]("Fold Change")))+
ylab(expression(-log[10]("P Value")))+
theme_bw()+
theme(legend.title = element_blank(),
plot.title = element_text(hjust = 0.5,size = 24, face= "bold"),
axis.title = element_text(size=20),
legend.text= element_text(size=20),
aspect.ratio=1)+
geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6) +
geom_text_repel(data=. %>%filter(symbol %in% c("hsa-mir-301b")),
aes(label=symbol),
size=3,
segment.color="black",
show.legend = F)
}
report("\n***************************构建ceRNA网络***************************")
miRNA_lnRNA <- readRDS("result/miRNA_lnRNA.rda")
miRNA_mRNA <- readRDS("result/miRNA_mRNA.rda")
lnRNA_miRNA_mRNA <- merge(miRNA_lnRNA,miRNA_mRNA,by="miRNA")
if(F){
miRNA <- unique(lnRNA_miRNA_mRNA$miRNA)
length(miRNA)
#Volcano_sig(readRDS("result/miRNA_expr_DESeq.rda"),gene=miRNA)
lnRNA <- unique(lnRNA_miRNA_mRNA$lnRNA)
length(lnRNA)
#Volcano_sig(readRDS("result/lnRNA_expr_DESeq.rda"),gene=lnRNA)
mRNA <-unique(lnRNA_miRNA_mRNA$mRNA)
length(mRNA)
}
if(T){
lncRNA_deg_ce <- lncRNA_deg_sig%>%
plyr::rename(c(symbol="lnRNA"))%>%
mutate(lnRNA_direction=ifelse(logFC>0,"up","down"))%>%
dplyr::select(lnRNA,lnRNA_direction)
miRNA_deg_ce <- miRNA_deg_sig%>%
plyr::rename(c(symbol="miRNA"))%>%
mutate(miRNA_direction=ifelse(logFC>0,"up","down"))%>%
dplyr::select(miRNA,miRNA_direction)
mRNA_deg_ce <- mRNA_deg_sig%>%
plyr::rename(c(symbol="mRNA"))%>%
mutate(mRNA_direction=ifelse(logFC>0,"up","down"))%>%
dplyr::select(mRNA,mRNA_direction)
lnRNA_miRNA_mRNA_ <- lnRNA_miRNA_mRNA
lnRNA_miRNA_mRNA <- lnRNA_miRNA_mRNA_ %>%
inner_join(lncRNA_deg_ce,by="lnRNA")%>%
inner_join(miRNA_deg_ce,by="miRNA")%>%
inner_join(mRNA_deg_ce,by="mRNA")%>%
filter(lnRNA_direction!=miRNA_direction)%>%
filter(miRNA_direction!=mRNA_direction)
}
#saveRDS(lnRNA_miRNA_mRNA,file = "result/lnRNA_miRNA_mRNA.rda")
lnRNA_miRNA_mRNA <- readRDS("result/lnRNA_miRNA_mRNA.rda")
unique(lnRNA_miRNA_mRNA$miRNA)
length(unique(lnRNA_miRNA_mRNA$lnRNA))
ce_lnRNA <- unique(lnRNA_miRNA_mRNA$lnRNA)
ce_lnRNA <- ce_lnRNA[ce_lnRNA %in% rownames(lnRNA_expr)]
ce_miRNA <- unique(lnRNA_miRNA_mRNA$miRNA)
ce_miRNA <- ce_lnRNA[ce_miRNA %in% rownames(mRNA_expr)]
ce_mRNA <- unique(lnRNA_miRNA_mRNA$mRNA)
ce_mRNA <- ce_lnRNA[ce_miRNA %in% rownames(mRNA_expr)]
cytoscape <- function(lnRNA_miRNA_mRNA,filename){
ce_lnRNA <- unique(lnRNA_miRNA_mRNA$lnRNA)
ce_miRNA <- unique(lnRNA_miRNA_mRNA$miRNA)
ce_mRNA <- unique(lnRNA_miRNA_mRNA$mRNA)
report("ceRNA网络中有: ", length(ce_lnRNA),
" 个lnRNA对应 ",length(ce_miRNA)," 个miRNA, ",
length(ce_miRNA)," 的miRNA对应 ",
length(ce_mRNA), "的mRNA")
miRNA_lnRNA <<- unique(lnRNA_miRNA_mRNA%>%dplyr::select(lnRNA,miRNA))
miRNA_mRNA <<- unique(lnRNA_miRNA_mRNA%>%dplyr::select(miRNA,mRNA))
lncRNA_deg_ce <- lncRNA_deg_sig%>%
plyr::rename(c(symbol="lnRNA"))%>%
filter(lnRNA %in% ce_lnRNA)%>%
mutate(direction=ifelse(logFC>0,"up","down"),type="lncRNA")%>%
dplyr::select(name=lnRNA,direction,type)
miRNA_deg_ce <- miRNA_deg_sig%>%
plyr::rename(c(symbol="miRNA"))%>%
filter(miRNA %in% ce_miRNA)%>%
mutate(direction=ifelse(logFC>0,"up","down"),type="miRNA")%>%
dplyr::select(name=miRNA,direction,type)
mRNA_deg_ce <- mRNA_deg_sig%>%
plyr::rename(c(symbol="mRNA"))%>%
filter(mRNA %in% ce_mRNA)%>%
mutate(direction=ifelse(logFC>0,"up","down"),type="mRNA")%>%
dplyr::select(name=mRNA,direction,type)
cytoscape_type <- bind_rows(mRNA_deg_ce,miRNA_deg_ce,lncRNA_deg_ce)%>%
mutate(type = str_c(direction,type,sep="_"))%>%
dplyr::select(-2)
write.csv(cytoscape_type,file = paste0("result/",filename,"_type.csv"),row.names = F,quote = F)
report("写入网络节点类型到 ",paste0("result/",filename,"_type.csv"), " 共有 ",dim(cytoscape_type)[1]," 个")
ceRAN_pair1 <- miRNA_lnRNA%>%
dplyr::select(miRNA,name=lnRNA)
ceRAN_pair2 <- miRNA_mRNA%>%
dplyr::select(miRNA,name=mRNA)
cytoscape_input <- bind_rows(ceRAN_pair1,ceRAN_pair2)
write.csv(cytoscape_input,file = paste0("result/",filename,"_input.csv"),row.names = F,quote = F)
report("写入网络节点关系:",paste0("result/",filename,"_input.csv")," 共有 ",dim(cytoscape_input)[1]," 个")
}
cytoscape(lnRNA_miRNA_mRNA,filename = "geo__")
a <- unique(lnRNA_miRNA_mRNA$mRNA)
length(a)
write.csv(a,file = "result/gene_all.csv",quote = F,row.names = F)
genes <- "MDM2,CCN1,HIF1A,TGFB1,MYCN,MXD1,NSD3,BTG1,PPP3R1,PPP6C"
genes <- unlist(str_split(genes,pattern = ","))
lnRNA_miRNA_mRNA_top10 <- lnRNA_miRNA_mRNA%>%
filter(mRNA %in% genes)
cytoscape(lnRNA_miRNA_mRNA_top10,filename = "geo_top10")
genes %in% unique(lnRNA_miRNA_mRNA_top10$mRNA)
unique(lnRNA_miRNA_mRNA_top10$mRNA)
unique(lnRNA_miRNA_mRNA_top10$lnRNA)
#Volcano_sig(readRDS("result/mRNA_expr_DESeq.rda"),gene=mRNA)
report("***************************构建ceRNA网络***************************")
#############################################
## lncRNA miRNA mRNA 生存分析
#############################################
report("\n***************************lncRNA miRNA mRNA 生存分析***************************")
library(survival)
library(survminer)
survival_p <- function(data,gene,prefix,pcutoff=0.05){
low_high <- ifelse(data[,gene]<=median(data[,gene]),"Low","Hight")
diff <- survdiff(Surv(futime_year,fustat)~low_high,data = data)
pVal <- 1 -pchisq(diff$chisq,df =1)
pValue <- signif(pVal,4)
if(pValue<pcutoff){
fit <- surv_fit(Surv(futime_year,fustat) ~ low_high, data = data)
ggsurvplot(fit, pval=T, risk.table=F,
legend.title=gene,
risk.table.height = 0.3,data = data)+
ggsave(filename = paste0("figures/survival/",prefix,"-",gene,".png"))
}
return(pValue)
}
Survival <- function(expr,gene,clinical,prefix,pcutoff=0.05){
expr_survival <- expr[gene,]%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(clinical,by="patient_id")
mRNA_survival <- sapply(gene,
function(x){ survival_p(expr_survival,x,prefix,pcutoff)})
return(mRNA_survival)
}
mRNA_lnRNA_clinical <- readRDS("result/mRNA_lnRNA_clinical.rda")
res <- Survival(mRNA_expr,ce_mRNA,mRNA_lnRNA_clinical,prefix="mRNA")
report("使用[R pacakge Survival]得到 ",length(ce_mRNA)," 个mRNA中(以p<0.05筛选)有 ",sum(res<0.05)," 个mRNA与总生存率相关")
report(paste(names(res[res<0.05]),collapse = ","))
## test start ceRNAw网络mRNA的表达与患者总生存率之间关系
if(F){
mRNA_expr_survival <- mRNA_expr[mRNA,]%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(mRNA_lnRNA_clinical,by="patient_id")
survival_p(mRNA_expr_survival,"PHOX2B")
mRNA_survival <- sapply(mRNA,
function(x){ survival_p(mRNA_expr_survival,x)})
table(mRNA_survival<0.05)
mRNA_expr_survival <- mRNA_expr_survival%>% filter(futime_year!=0)
low_high <- ifelse(mRNA_expr_survival[,"PHOX2B"] <= median(mRNA_expr_survival[,"PHOX2B"]),"Low","Hight")
mRNA_expr_survival$futime_year <- as.numeric(mRNA_expr_survival$futime_year)
mRNA_expr_survival$fustat <- as.numeric(mRNA_expr_survival$fustat)
fit <- surv_fit(Surv(futime_year,fustat) ~ low_high, data = mRNA_expr_survival)
# pdf(file = paste0("result/PHOX2B.pdf"),width = 5.5,height = 5)
ggsurvplot(fit, pval=T, risk.table=F,
legend.title='PHOX2B',
risk.table.height = 0.3,data = mRNA_expr_survival)+
ggsave(filename = "result/aa.png")
# dev.off()
# ggsurvplot(fit, data = mRNA_expr_survival,
# surv.median.line = "hv", # 增加中位生存时间
# conf.int = TRUE) # 增加置信区间
}
## test End
mRNA_lnRNA_clinical <- readRDS("result/mRNA_lnRNA_clinical.rda")
res <- Survival(lnRNA_expr,ce_lnRNA,mRNA_lnRNA_clinical,prefix="lnRNA")
report("使用[R pacakge Survival]得到 ",length(ce_lnRNA)," 个lnRNA中(以p<0.05筛选)有 ",sum(res<0.05)," 个lnRNA与总生存率相关")
report(paste(names(res[res<0.05]),collapse = ","))
## 由于miRNA存在hsa-let-7a-2 这种形式，预测mRNA时是去掉-2 ，以hsa-let-7a预测
miRNA_clinical <- readRDS("result/miRNA_clinical.rda")
res <- Survival(miRNA_expr,ce_miRNA,miRNA_clinical,prefix="miRNA",pcutoff = 0.09)
report("使用[R pacakge Survival]得到 ",length(ce_miRNA)," 个lnRNA中(以p<0.05筛选)有 ",sum(res<0.09)," 个lnRNA与总生存率相关")
report(paste(names(res[res<0.09]),collapse = ","))
report("***************************lncRNA miRNA mRNA 生存分析***************************")
#############################################
## cox 回归分析
#############################################
library(survivalROC)
library(survival)
library(survminer)
Cox_single <- function(data,gene,filename=NULL){
outTab <- data.frame()
for(i in gene){
cox <- coxph(Surv(futime_year,fustat) ~ data[,i],data=data)
coxSummary <- summary(cox)
outTab <- rbind(outTab,cbind(id=i,
HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
}
if(!is.null(filename)){
write.csv(outTab,file = paste0("result/",filename,".csv"),quote = F,row.names = F)
}
return(outTab)
}
Cox_multi <- function(data,gene,filename=NULL){
if(sum(grepl("-",gene))>0){
gene <- sapply(gene,function(x){sub("-","",x)})
data<- data%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(gene,collapse = "+")))
cox <- coxph(formula,data=data)
#cox <- step(cox,direction = "both")
coxSummary <- summary(cox)
outTab <- data.frame()
outTab <- rbind(outTab,cbind(HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
riskScore <- predict(cox,type = "risk",newdata = data)
expr_survival_riskScore <- data %>%
mutate(riskScore = as.vector(ifelse(riskScore>median(riskScore),"high","low")))
fit <- surv_fit(Surv(futime_year,fustat) ~ riskScore, data = expr_survival_riskScore)
res <- ggsurvplot(fit, pval=T, risk.table=F,
risk.table.height = 0.3,data = expr_survival_riskScore)
if(!is.null(filename)){
res+ ggsave(filename = paste0("figures/cox/",filename,".png"))
write.csv(outTab,file = paste0("result/",filename,".csv"),quote = F,row.names = F)
}
return(list(outTab=outTab,image=res))
}
riskScore_pheatmap <- function(data,gene,filename=NULL){
if(sum(grepl("-",gene))>0){
gene <- sapply(gene,function(x){sub("-","",x)})
data<- lnRNA_expr_survival%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(gene,collapse = "+")))
cox <- coxph(formula,data=data)
#cox <- step(cox,direction = "both")
riskScore <- predict(cox,type = "risk",newdata = data)
expr_survival_riskScore <- data %>%
mutate(riskScore = as.vector(ifelse(riskScore>median(riskScore),"high","low")))
deg_expr <- expr[deg_gene,]%>%
{log10(cpm(.)+0.01)}
#deg_expr <- t(scale(t(deg_expr)))
#deg_expr[deg_expr > 2] = 2
#deg_expr[deg_expr < -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,show_rownames = F,
filename = paste0("figures/cox/",filename,".png"))
return(paste0("figures/",filename,".png"))
}
Cox_multi_ROC <- function(data,gene,filename=NULL,predict.time=5){
if(sum(grepl("-",gene))>0){
gene <- sapply(gene,function(x){sub("-","",x)})
data<- data%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(gene,collapse = "+")))
cox <- coxph(formula,data=data)
#message(summary(cox))
#cox <- step(cox,direction = "both")
riskScore <- predict(cox,type = "risk",newdata = data)
roc <- survivalROC(Stime=data$futime_year ,status=data$fustat,marker=riskScore,
predict.time=predict.time,method="KM")
res <- xyplot(TP~FP,roc,main=paste("ROC curve(","AUC=",round(roc$AUC,3),")"),
xlab="False postive rate",ylab="True positive rate",panel = function(x,y){
panel.xyplot(x,y,type="l",lwd=2,cex.main=1.3,cex.lab=1.2,cex.axis=1.2,font=1.2,col="red")
panel.abline(c(0,1))
})
if(!is.null(filename)){
png(filename = paste0("figures/cox/",filename,".png"))
print(res)
dev.off()
}
return(res)
}
lnRNA_expr <- readRDS("result/lncRNA_expr.rda")
mRNA_lnRNA_clinical <- readRDS("result/mRNA_lnRNA_clinical.rda")
lnRNA_expr_survival <- lnRNA_expr[lnRNA_intersect,]%>%
filter()%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(mRNA_lnRNA_clinical,by="patient_id")
if(F){
lnRNA_expr_ <- lnRNA_expr[ce_lnRNA,]
#lnRNA_expr_ <- log2(lnRNA_expr_+0.01)
par(mfcol=c(2,1))
boxplot(lnRNA_expr_)
boxplot(log2(lnRNA_expr_+1))
boxplot( log2(lnRNA_expr_+0.01))
boxplot( log2(lnRNA_expr_+1))
boxplot(scale(cpm(lnRNA_expr_),center = F, scale = TRUE))
lnRNA_expr_survival <- lnRNA_expr_%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(mRNA_lnRNA_clinical,by="patient_id")
Cox_single_res <- Cox_single(lnRNA_expr_survival,ce_lnRNA,filename = "singleCox_lnRNA")
genes <- Cox_single_res%>%filter(pvalue<0.09)%>%pull("id")
Cox_multi(lnRNA_expr_survival,genes)$image
Cox_multi_ROC(lnRNA_expr_survival,genes,predict.time=3)
}
## 22个lnRNA 先单因素cox再多因素cox LINC00486,LINC00365
Cox_single_res <- Cox_single(lnRNA_expr_survival,lnRNA_intersect,filename = "singleCox_lnRNA")
genes <- Cox_single_res%>%filter(pvalue<0.09)%>%pull("id")
saveRDS(genes,file = "result/lnRNA_cox_mutil.rda")
res <- Cox_multi(lnRNA_expr_survival,genes,filename = "muitiCox_lnRNA")
res$outTab
print(res$image)
Cox_multi_ROC(lnRNA_expr_survival,genes,predict.time=3,filename = "muitiCox_lnRNA_lnRNA_ROC")
if(F){
cox <- coxph(Surv(futime_year,fustat) ~ lnRNA_expr_survival[,"HOTAIR"],data=lnRNA_expr_survival)
coxSummary <- summary(cox)
outTab <- data.frame()
outTab <- rbind(outTab,cbind(HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
wald.test <- signif(coxSummary$wald["pvalue"],digits = 3)
beta <- signif(coxSummary$wald["test"],digits = 3)
HR <- signif(coxSummary$coef[1],digits = 3)
HR.lower = coxSummary$conf.int[,"lower .95"]
HR.upper = coxSummary$conf.int[,"upper .95"]
}
if(F){
if(sum(grepl("-",genes))>0){
genes <- sapply(genes,function(x){sub("-","",x)})
lnRNA_expr_survival<- lnRNA_expr_survival%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(genes,collapse = "+")))
cox <- coxph(formula,data=lnRNA_expr_survival)
cox <- step(cox,direction = "both")
coxSummary <- summary(cox)
outTab <- data.frame()
outTab <- rbind(outTab,cbind(HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
riskScore <- predict(cox,type = "risk",newdata = lnRNA_expr_survival)
lnRNA_expr_survival_riskScore <- lnRNA_expr_survival %>%
mutate(riskScore = as.vector(ifelse(riskScore< median(riskScore),"high","low")))
fit <- surv_fit(Surv(futime_year,fustat) ~ riskScore, data = lnRNA_expr_survival_riskScore)
ggsurvplot(fit, pval=T, risk.table=F,
risk.table.height = 0.3,data = lnRNA_expr_survival_riskScore,xlab="Time(Years)")
roc <- survivalROC(Stime=lnRNA_expr_survival_riskScore$futime_year ,status=lnRNA_expr_survival_riskScore$fustat,marker=riskScore,
predict.time=5,method="KM")
plot(roc$FP,roc$TP,type="l",xlim=c(0,1),ylim=c(0,1),col="red",main=paste("ROC curve(","AUC=",round(roc$AUC,3),")"),
xlab="False postive rate",ylab="True positive rate",
lwd=2,cex.main=1.3,cex.lab=1.2,cex.axis=1.2,font=1.2)
abline(0,1)
lnRNA_expr_survival <- lnRNA_expr_survival%>%
select(c(genes,futime_year,fustat))
library(ggrisk)
ggrisk(lnRNA_expr_survival,time="futime_year",event="fustat",heatmap.genes = genes)
}
if(F){
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(genes,collapse = "+")))
cox <- coxph(formula,data=lnRNA_expr_survival)
riskScore <- predict(cox,type = "risk",newdata = lnRNA_expr_survival)
expr_survival_riskScore <- lnRNA_expr_survival %>%
mutate(riskScore = as.vector(ifelse(riskScore>median(riskScore),"high","low")))
lnRNA_expr_group <- data.frame(row.names = expr_survival_riskScore$TCGA_full_id,group=expr_survival_riskScore$riskScore)%>%
arrange(desc(group))
lnRNA_expr <- lnRNA_expr[,rownames(lnRNA_expr_group)]
identical(rownames(lnRNA_expr_group),colnames(lnRNA_expr))
deg_expr <- lnRNA_expr[genes,]%>%
{log2(cpm(.)+ 1)}
deg_expr <- t(scale(t(deg_expr)))
deg_expr[deg_expr>2] = 2
deg_expr[deg_expr< -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,scale = "row",annotation_col=lnRNA_expr_group,cluster_cols=F)
}
if(F){
## 单基因
my_comparisons <- list( c("i", "ii"))
lnRNA_expr_survival%>%
dplyr::select(c(tumor_stage,LINC00365))%>%
filter(tumor_stage!="unknow") %>%
ggplot(aes(x=tumor_stage,y=LINC00365))+
stat_boxplot(geom="errorbar",width=0.15,aes(color=tumor_stage))+
geom_boxplot(aes(fill=tumor_stage))+
ylim(0, 1000)+
stat_compare_means(comparisons =my_comparisons)+
stat_compare_means(method = "kruskal.test",label = "p.format")
}
mRNA_lnRNA_clinical <- readRDS("result/mRNA_lnRNA_clinical.rda")
mRNA_expr_survival <- mRNA_expr[ce_mRNA,]%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
#filter(substr(TCGA_full_id,14,15)!=11)%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(mRNA_lnRNA_clinical,by="patient_id")
Cox_single_res <- Cox_single(mRNA_expr_survival,ce_mRNA,filename = "singleCox_mRNA")
genes <- Cox_single_res%>%filter(pvalue<0.05)%>%pull("id")
res <- Cox_multi(mRNA_expr_survival,genes,filename = "muitiCox_mRNA")
res$outTab
#print(res$image)
Cox_multi_ROC(mRNA_expr_survival,genes,predict.time=5,filename = "muitiCox_mRNA_ROC")
if(F){
if(sum(grepl("-",genes))>0){
genes <- sapply(genes,function(x){sub("-","",x)})
data<- mRNA_expr_survival%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(genes,collapse = "+")))
cox <- coxph(formula,data=data)
riskScore <- predict(cox,type = "risk",newdata = data)
expr_survival_riskScore <- data %>%
mutate(riskScore = as.vector(ifelse(riskScore>median(riskScore),"high","low")))
lnRNA_expr_group <- data.frame(row.names = expr_survival_riskScore$TCGA_full_id,group=expr_survival_riskScore$riskScore)%>%
arrange(desc(group))
lnRNA_expr <- lnRNA_expr[,rownames(lnRNA_expr_group)]
identical(rownames(lnRNA_expr_group),colnames(lnRNA_expr))
deg_expr <- lnRNA_expr[genes,]%>%
{log2(cpm(.)+ 1)}
deg_expr <- t(scale(t(deg_expr)))
deg_expr[deg_expr>2] = 2
deg_expr[deg_expr< -2] = -2
col <- colorRampPalette(c("green","white","red"))(256)
pheatmap(deg_expr,col=col,fontsize=8,show_colnames =F,scale = "row",annotation_col=lnRNA_expr_group,cluster_cols=F)
}
if(F){
if(sum(grepl("-",genes))>0){
genes <- sapply(genes,function(x){sub("-","",x)})
lnRNA_expr_survival<- lnRNA_expr_survival%>%
rename_all(~sub("-","",.))
}
formula <-as.formula(paste0("Surv(futime_year,fustat)~",paste0(genes,collapse = "+")))
cox <- coxph(formula,data=lnRNA_expr_survival)
cox <- step(cox,direction = "both")
coxSummary <- summary(cox)
outTab <- data.frame()
outTab <- rbind(outTab,cbind(HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
riskScore <- predict(cox,type = "risk",newdata = lnRNA_expr_survival)
lnRNA_expr_survival_riskScore <- lnRNA_expr_survival %>%
mutate(riskScore = as.vector(ifelse(riskScore< median(riskScore),"high","low")))
fit <- surv_fit(Surv(futime_year,fustat) ~ riskScore, data = lnRNA_expr_survival_riskScore)
ggsurvplot(fit, pval=T, risk.table=F,
risk.table.height = 0.3,data = lnRNA_expr_survival_riskScore,xlab="Time(Years)")
roc <- survivalROC(Stime=lnRNA_expr_survival_riskScore$futime_year ,status=lnRNA_expr_survival_riskScore$fustat,marker=riskScore,
predict.time=5,method="KM")
xyplot(TP~FP,roc,main=paste("ROC curve(","AUC=",round(roc$AUC,3),")"),
xlab="False postive rate",ylab="True positive rate",panel = function(x,y){
panel.xyplot(x,y,type="l",lwd=2,cex.main=1.3,cex.lab=1.2,cex.axis=1.2,font=1.2,col="red")
panel.abline(c(0,1))
})
lnRNA_expr_survival <- lnRNA_expr_survival%>%
select(c(genes,futime_year,fustat))
library(ggrisk)
ggrisk(lnRNA_expr_survival,time="futime_year",event="fustat",heatmap.genes = genes)
}
if(F){
my_comparisons <- list( c("i", "ii"))
mRNA_expr_survival%>%
dplyr::select(c(tumor_stage,CNTN3))%>%
filter(tumor_stage!="unknow") %>%
ggplot(aes(x=tumor_stage,y=CNTN3))+
stat_boxplot(geom="errorbar",width=0.15,aes(color=tumor_stage))+
geom_boxplot(aes(fill=tumor_stage))+
stat_compare_means(comparisons =my_comparisons)+
stat_compare_means(method = "kruskal.test",label = "p.format")
}
#############################################
## ceRNA 子网络构建
#############################################
lnRNA_cox_mutil <- readRDS("result/lnRNA_cox_mutil.rda")
lnRNA_miRNA_mRNA <- readRDS("result/lnRNA_miRNA_mRNA.rda")
lnRNA_miRNA_mRNA_sub <- lnRNA_miRNA_mRNA%>%
filter(lnRNA %in% lnRNA_cox_mutil)
cytoscape(lnRNA_miRNA_mRNA_sub,filename = "small")
saveRDS(lnRNA_miRNA_mRNA_sub,file = "result/sub_netwrok_gene.rda")
if(F){
miRNA_lnRNA_clinical <- readRDS("result/mRNA_lnRNA_clinical.rda")
miRNA_expr <- readRDS("result/miRNA_expr.rda")%>%
rownames_to_column("symbol")%>%
mutate(miRNA_= case_when(grepl("hsa-mir",symbol)~str_extract(symbol,"hsa-mir-[0-9]+[a-z]?"),
grepl("hsa-let",symbol)~str_extract(symbol,"hsa-let-[0-9]+[a-z]?")))%>%
filter(miRNA_ %in% miRNA)%>%
dplyr::select(-miRNA_)%>%
column_to_rownames("symbol")
miRNA_expr_survival <- miRNA_expr%>%
t()%>%
as.data.frame()%>%
rownames_to_column("TCGA_full_id")%>%
mutate(patient_id =substr(TCGA_full_id,1,12))%>%
inner_join(mRNA_lnRNA_clinical,by="patient_id")
Cox_single_res <- Cox_single(miRNA_expr_survival, rownames(miRNA_expr))
genes <- Cox_single_res%>%filter(pvalue<0.05)%>%pull("id")
Cox_multi(mRNA_expr_survival,genes,filename = "ce-lnRNA-Cox_multi")
Cox_multi_ROC(miRNA_expr_survival,predict.time=3,filename = "ce-lnRNA-Cox_multi")
}
##############################################
## GO KEGG
##############################################
(function(){
if(!file.exists("enrichment")){
dir.create(file.path("enrichment"))
}
sub_netwrok_gene <- readRDS("result/sub_netwrok_gene.rda")
if(F){
gene <- unique(sub_netwrok_gene[sub_netwrok_gene$miRNA=="hsa-mir-34b",]$mRNA)
write.csv(gene,file = "result/GO_KEKK_GENE.csv",quote = F,row.names = F)
}
gene <- unique(sub_netwrok_gene$mRNA)
entrez_id <- bitr(gene, fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db)%>%
pull("ENTREZID")%>%unique()
egobp_BP<-enrichGO(gene = entrez_id,OrgDb = org.Hs.eg.db,keyType = 'ENTREZID',ont = 'BP',pAdjustMethod="none") # BP MF CC ALL
#dotplot(egobp_BP,showCategory=20)
egobp_BP@result%>%
arrange(p.adjust)%>%
write.csv(file = "result/enrichment/egobp_BP.csv",quote = F,row.names = F)
write.csv(egobp_BP@result)
egobp_MF<-enrichGO(gene = entrez_id,OrgDb = org.Hs.eg.db,keyType = 'ENTREZID',ont = 'MF',pAdjustMethod="none") # BP MF CC ALL
#dotplot(egobp_MF,showCategory=20)
egobp_MF@result%>%
arrange(p.adjust)%>%
write.csv(file = "result/enrichment/egobp_MF.csv",quote = F,row.names = F)
egobp_CC<-enrichGO(gene = entrez_id,OrgDb = org.Hs.eg.db,keyType = 'ENTREZID',ont = 'CC',pAdjustMethod="none") # BP MF CC ALL
#dotplot(egobp_CC,showCategory=20)
egobp_CC@result%>%
arrange(p.adjust)%>%
write.csv(file = "result/enrichment/egobp_CC.csv",quote = F,row.names = F)
kk <- clusterProfiler::enrichKEGG(gene =entrez_id,organism = 'hsa',pvalueCutoff = 0.4,qvalueCutoff = 0.3,pAdjustMethod="none")
#dotplot(kk,showCategory=20)
kk@result%>%
arrange(p.adjust)%>%
write.csv(file = "result/enrichment/kegg.csv",quote = F,row.names = F)
})()