GDCRNATools: integrative analysis of protein coding genes, long non-coding genes, and microRNAs in GDCgithub
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
# 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
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()
https://wilkelab.org/ggridges/
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) })()