代码分享：基于TCGA甲基化数据挖掘验证分子表型代码分享（第三期）
2026-03-10 来源:本站 点击次数:94
DNA甲基化作为表观遗传调控的重要机制，在肿瘤发生、发展及预后中发挥关键作用。前两期已经完成了甲基化数据挖掘的基础分析，这一期将立足于解释具体的生物学问题来展开。
8.每个基因的甲基化信号值的比较 --------------------------------------------------------
结果如下：
结果如下：
结果如下：
结果如下：
结果如下：
结果如下：
结果如下：
最后出图结果如下：
Fig1样本质控
Fig2样本筛选
Fig3-标志物筛选
本文所使用的数据均已公开，欢迎各位老师获取数据下载链接。
8.每个基因的甲基化信号值的比较 --------------------------------------------------------
# 8.每个基因的甲基化信号值的比较 --------------------------------------------------------
##对膀胱癌中差异上调的cg进行绘制均值的折线图
load('./00.data//pangguangai_21paired_meLoad_data_412481.Rdata')
# exp_count_gene_name =myLoad$beta
sample= pd
normalized_counts =myNorm
# exp_count_gene_name = exp_count_gene_name[,c(1:10,13:18,21:24,27:30,33:42)]
sample=sample[c(1:10,13:18,21:24,27:30,33:42),]
normalized_counts=normalized_counts[,c(1:10,13:18,21:24,27:30,33:42)]
colnames(marker_plot)
marker_1 = marker[,c(1,8,9,15,25)]
marker_plot <- reshape2::melt(marker_1, id.vars=c("rowname","gene",'change'), variable.name="group",value.name="AVG")
color <- c("#E63863","#0228A5")
gene <- "IFFO1"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "IFFO1", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_IFFO1",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_IFFO1",".png"), width = 6*1.5, height = 6*1.5)
##对膀胱癌中差异上调的cg进行绘制均值的折线图
load('./00.data//pangguangai_21paired_meLoad_data_412481.Rdata')
# exp_count_gene_name =myLoad$beta
sample= pd
normalized_counts =myNorm
# exp_count_gene_name = exp_count_gene_name[,c(1:10,13:18,21:24,27:30,33:42)]
sample=sample[c(1:10,13:18,21:24,27:30,33:42),]
normalized_counts=normalized_counts[,c(1:10,13:18,21:24,27:30,33:42)]
colnames(marker_plot)
marker_1 = marker[,c(1,8,9,15,25)]
marker_plot <- reshape2::melt(marker_1, id.vars=c("rowname","gene",'change'), variable.name="group",value.name="AVG")
color <- c("#E63863","#0228A5")
gene <- "IFFO1"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "IFFO1", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_IFFO1",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_IFFO1",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
gene <- "TMEM106A"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TMEM106A", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_TMEM106A",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_TMEM106A",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TMEM106A", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_TMEM106A",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_TMEM106A",".png"), width = 6*1.5, height = 6*1.5)
gene <- "SHANK2"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "SHANK2", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_SHANK2",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_SHANK2",".png"), width = 6*1.5, height = 6*1.5)
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "SHANK2", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_SHANK2",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_SHANK2",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
gene <- "IFFO1"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "IFFO1", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_IFFO1",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_IFFO1",".png"), width = 6*1.5, height = 6*1.5)
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "IFFO1", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_IFFO1",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_IFFO1",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
gene <- "TMEM106A"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TMEM106A", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_TMEM106A",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_TMEM106A",".png"), width = 6*1.5, height = 6*1.5)
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TMEM106A", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_TMEM106A",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_TMEM106A",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
gene <- "SHANK2"
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "SHANK2", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_SHANK2",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_SHANK2",".png"), width = 6*1.5, height = 6*1.5)
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "SHANK2", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_SHANK2",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_SHANK2",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
gene <- c("IFFO1","TMEM106A","SHANK2")
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
marker_plot_sub
marker_plot_sub$rowname <- factor(marker_plot_sub$rowname, levels = c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,
"cg18222083", "cg19548479" ,
"cg20170028"
))
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TCGA BLCA", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_all",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_all",".png"), width = 6*1.5, height = 6*1.5)
marker_plot_sub = marker_plot[marker_plot$gene %in% gene,]
marker_plot_sub
marker_plot_sub$rowname <- factor(marker_plot_sub$rowname, levels = c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,
"cg18222083", "cg19548479" ,
"cg20170028"
))
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_line(size=0.5)+
geom_point(size=3)+
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TCGA BLCA", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "line_βval_all",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "line_βval_all",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
color <- c("#EE9B7E","#D9D7BC")
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TCGA BLCA", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_all",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_all",".png"), width = 6*1.5, height = 6*1.5)
ggplot(marker_plot_sub, aes(x = rowname, y = AVG, group = group, colour=group, shape= group, fill = group)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_text(aes(label=AVG),position = position_dodge(.9),vjust =-1,hjust=1.3)+
scale_shape_manual(values = c(24,22)) +
scale_fill_manual(values =color)+
scale_color_manual(values = color)+
labs(x="",y="β value", title = "TCGA BLCA", subtitle = "Padj < 0.00001 & |logFC| >= 0.4") +
# ylim(0,3200) +
theme(
#legend.position="none",
axis.text.x = element_text(size = 12,colour="black"),
#axis.text.x = element_text(size = 12,colour="black", angle = 45,hjust = 1),
axis.text.y = element_text(size = 12,colour="black"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 12),
#axis.ticks = element_blank(),
plot.title = element_text(hjust = 0.5, vjust=1),
#panel.background = element_rect(fill="white",color = "black"),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA))
ggsave(paste0("../03.plot/", "barplot_βval_all",".pdf"), width = 6*1.5, height = 6*1.5)
ggsave(paste0("../03.plot/", "barplot_βval_all",".png"), width = 6*1.5, height = 6*1.5)
结果如下：
#上面的结果再加上对应的igv的图这里我选择hg19
9.补充生存分析以及AUC分析 ----------------------------------------------
# 输入文件
dd =read.delim('../00.data/clinical.tsv')
# clinical.tsv，包含生存数据和临床亚组信息。
# dd <- data.table::fread("easy_input.csv",data.table = F)
dd[1:3,]
#处理数据为生存分析需要的格式
meta = dd
# 生存时间需要为月份
# meta[meta==""]=NA
meta$time = ifelse(meta$vital_status == 'Alive',meta$days_to_last_follow_up,meta$days_to_death)
meta$time = as.numeric(meta$time)/30
# 1.2生死定义
meta$event = ifelse(meta$vital_status== 'Alive',0,1)
table(meta$event)
# 区分年龄
meta$age = ceiling(abs(as.numeric(meta$days_to_birth))/365)
meta$age_group=ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')
table(meta$age_group)
## 自定义颜色
orange <- "#EB736B"
blue <- "#57B5AB"
library(stringr)
# 1.4提取分期信息
b= meta$ajcc_pathologic_stage
a= str_extract_all(meta$ajcc_pathologic_stage,"I|V");head(a)#提取分期信息
b= sapply(a,paste,collapse = '')
head(b)
meta$stage = b
#去掉生存信息不全或者生存时间小于0.1(月)的病人，样本纳排标准不唯一，且差别很大
# k1 = meta$time>=0.1;table(k1)
# # 8
# k2 =!(is.na(meta$time)|is.na(meta$event));table(k2)
# 6
##只先看目前的数据
#注意生存分析只需要肿瘤的病人
head(normalized_counts)
# nm_tumor = normalized_counts[]
sample$sample_type_fix[sample$sample_type_fix=='Primary Tumor']='Tumor'
colnames(sample)
exp = normalized_counts[,sample$sample_type_fix %in% "Tumor"]
dim(exp)#412481 17 只存贮这17个
meta$case_submitter_id<- str_replace_all(meta$case_submitter_id ,"-", ".")
# meta$case_submitter_id =
# meta = meta[meta$case_submitter_id %in% str_sub(colnames(exp),1,12),]
meta = meta[meta$case_id %in% sample_tumor$case_id,]
sample_tumor = sample[sample$sample_type_fix %in% c("Tumor"),]
dim(sample_tumor) #34 215
dim(meta) #34 215
meta <- meta[!duplicated(meta$case_id), ]
dim(meta) #17 215
# exp = exp[,k1&k2]
# meta = meta[k1&k2,]
# rownames(meta)
# colnames(meta)
save(meta,exp,sample_tumor,file = './/surve_model_tumor.Rdata')
# load('./00.data//surve_model_tumor.Rdata')
## 提取数据
# rt <- meta
# colnames(rt)
## 生存分析
# rt$days_to_death =as.numeric(rt$days_to_death)
fit <- survfit(Surv(time,event) ~age_group, data = meta)
## 作图
ggsurvplot(fit,
pval = TRUE,
palette = c(orange, blue),
title="age_group",
xlab="Time")
orange <- "#EB736B"
blue <- "#57B5AB"
ggsurvplot(fit,
palette =c("#E7B800","#2E9FDF"),
risk.table =TRUE,pVal =TRUE,
conf.int =TRUE,xlab ="Time in months",
ggtheme =theme_light(),
ncensor.plot = TRUE)
ggsurvplot(fit,
palette =c("#EB736B","#57B5AB"),
risk.table =TRUE,pVal =TRUE,
conf.int =TRUE,xlab ="Time in months",
ggtheme =theme_light(),
ncensor.plot = TRUE)
#性别年龄
sfit1=survfit(Surv(time, event)~gender, data=meta)
sfit2=survfit(Surv(time, event)~age_group, data=meta)
splots<- list()
splots[[1]]<- ggsurvplot(sfit1,pval =TRUE,data = meta,risk.table = TRUE)
splots[[2]]<- ggsurvplot(sfit2,pval =TRUE,data = meta,risk.table = TRUE)
arrange_ggsurvplots(splots, print = TRUE,
ncol=2,nrow=1,risk.table.height=0.4)
dev.off()
#单个基因(选择对应的cg展示他们的风险)
g=rownames(exp)[1]
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
#多个基因
gs=rownames(exp)[1:4]
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=2,nrow=2,risk.table.height=0.4)
##log_rank_p
#加上KM生存曲线即可
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event == 1表示“1为结局时间发生“
log_rank_p <- pbapply(exp, 1, function(gene){
# gene <- exp[1,]
meta$group <- ifelse(gene > median(gene), "high", "low")
data.survdiff <- survdiff(mySurv ~ group, data = meta)
# survdiff()用于检验两条或多条生存曲线之间是否存在差异，
# 或者针对已知替代项检验单条曲线，也适用于分组资料以及多组间生存曲线的比较；
p.val <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) # 计算log-rank检验的p值
return(p.val)
})
log_rank_p <- sort(log_rank_p) # 取出每一个基因生存分析的P值，形成表
table(log_rank_p < 0.01) # 哪些是P<0.01的
table(log_rank_p < 0.05) # 哪些是P<0.05的
log_rank <- names(log_rank_p)[log_rank_p < 0.05] # log-rank检验筛选的差异基因 0.05
log_rank_p
# Cox比例风险回归（HR）
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event = 1表示“1为结局事件发生”
unicox_results <- pbapply(exp, 1, function(gene){ # 对表达矩阵exp中的每一行（基因）执行函数function
# gene = as.numeric(exp[1,]) # 举个例子
group <- ifelse(gene > median(gene), "high","low") # 标记不同样本中该基因的上下调，将表达量转换为二分类变量
survival_dat <- data.frame(group = group,
stage = meta$stage,
age = meta$age,
gender = meta$gender,
gene = gene, # 基因表达量
stringsAsFactors = F) # 选择用于生存分析的数据
# 拟合Cox比例风险回归模型 m ：
m <- coxph(mySurv ~ gender + age + stage + group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ gene, data = survival_dat) # 表达量也可直接使用连续型变量gene
# sm <- summary(m) # 看拟合出的Cox比例风险回归模型的参数
# n：总人数。
# number of events：发生结局事件的人数，即死亡人数
# coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
# exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
# se(coef)：coef的标准误
# z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
# Pr(>|z|)：p值
# low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
tmp <- round(cbind(coef = beta, # 模型的偏回归系数β
se = se, # 偏回归系数β的标准误
z = beta/se, # Wald检验的统计量
p = 1 - pchisq((beta/se)^2, 1), # p值，差异的统计学意义
HR = exp(beta), # 风险值，计算公式为exp(coef)
HRCILL = exp(beta - qnorm(.975, 0, 1) * se), # HR的95%置信区间低点
HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3) # HR的95%置信区间高点
return(tmp["grouplow",]) # 分类变量
# return(tmp["gene",]) # 连续型变量
})
unicox_results <- as.data.frame(t(unicox_results))
table(unicox_results[,4] < 0.01) # p值<0.01的有多少
table(unicox_results[,4] < 0.05) # p值<0.05的有多少
unicox <- rownames(unicox_results)[unicox_results[,4] < 0.05] # Cox单因素分析筛选的差异基因
save(log_rank_p,unicox_results,file = './/surve_cox_log_rank_p.Rdata')
#多个基因
# c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028"
gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
load('./00.data//surve_cox_log_rank_p.Rdata')
# 取交集
cd =intersect(log_rank,gs) #0
cd =intersect(unicox,gs) #0
gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> mean(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
#ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5.pdf"), width = 6*1.2, height = 6*1)
##log_rank_p
# KM生存曲线即可
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event == 1表示“1为结局时间发生“
log_rank_p <- pbapply(exp, 1, function(gene){
# gene <- exp[1,]
meta$group <- ifelse(gene > mean(gene), "high", "low")
data.survdiff <- survdiff(mySurv ~ group, data = meta)
# survdiff()用于检验两条或多条生存曲线之间是否存在差异，
# 或者针对已知替代项检验单条曲线，也适用于分组资料以及多组间生存曲线的比较；
p.val <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) # 计算log-rank检验的p值
return(p.val)
})
log_rank_p <- sort(log_rank_p) # 取出每一个基因生存分析的P值，形成表
table(log_rank_p < 0.01) # 哪些是P<0.01的
table(log_rank_p < 0.05) # 哪些是P<0.05的
log_rank <- names(log_rank_p)[log_rank_p < 0.05] # log-rank检验筛选的差异基因 0.05
log_rank_p
# Cox比例风险回归（HR）
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event = 1表示“1为结局事件发生”
unicox_results <- pbapply(exp, 1, function(gene){ # 对表达矩阵exp中的每一行（基因）执行函数function
# gene = as.numeric(exp[1,]) # 举个例子
group <- ifelse(gene > mean(gene), "high","low") # 标记不同样本中该基因的上下调，将表达量转换为二分类变量
survival_dat <- data.frame(group = group,
stage = meta$stage,
age = meta$age,
gender = meta$gender,
gene = gene, # 基因表达量
stringsAsFactors = F) # 选择用于生存分析的数据
# 拟合Cox比例风险回归模型 m ：
m <- coxph(mySurv ~ gender + age + stage + group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ gene, data = survival_dat) # 表达量也可直接使用连续型变量gene
# sm <- summary(m) # 看拟合出的Cox比例风险回归模型的参数
# n：总人数。
# number of events：发生结局事件的人数，即死亡人数
# coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
# exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
# se(coef)：coef的标准误
# z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
# Pr(>|z|)：p值
# low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
tmp <- round(cbind(coef = beta, # 模型的偏回归系数β
se = se, # 偏回归系数β的标准误
z = beta/se, # Wald检验的统计量
p = 1 - pchisq((beta/se)^2, 1), # p值，差异的统计学意义
HR = exp(beta), # 风险值，计算公式为exp(coef)
HRCILL = exp(beta - qnorm(.975, 0, 1) * se), # HR的95%置信区间低点
HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3) # HR的95%置信区间高点
return(tmp["grouplow",]) # 分类变量
# return(tmp["gene",]) # 连续型变量
})
unicox_results <- as.data.frame(t(unicox_results))
table(unicox_results[,4] < 0.01) # p值<0.01的有多少
table(unicox_results[,4] < 0.05) # p值<0.05的有多少
unicox <- rownames(unicox_results)[unicox_results[,4] < 0.05] # Cox单因素分析筛选的差异基因
save(log_rank_p,unicox_results,file = './00.data//surve_cox_log_rank_p_mean.Rdata')
cd =intersect(log_rank,gs) #0
cd =intersect(unicox,gs) #0
#输出均值的结果
# cg19548479
write.csv(log_rank_p,'./00.data//log_rank_p_mean.csv')
write.csv(unicox_results,'./00.data//unicox_results_mean.csv')
# cg19548479符合后面的几个检验
# gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
g='cg19548479'
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
dd =read.delim('../00.data/clinical.tsv')
# clinical.tsv，包含生存数据和临床亚组信息。
# dd <- data.table::fread("easy_input.csv",data.table = F)
dd[1:3,]
#处理数据为生存分析需要的格式
meta = dd
# 生存时间需要为月份
# meta[meta==""]=NA
meta$time = ifelse(meta$vital_status == 'Alive',meta$days_to_last_follow_up,meta$days_to_death)
meta$time = as.numeric(meta$time)/30
# 1.2生死定义
meta$event = ifelse(meta$vital_status== 'Alive',0,1)
table(meta$event)
# 区分年龄
meta$age = ceiling(abs(as.numeric(meta$days_to_birth))/365)
meta$age_group=ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')
table(meta$age_group)
## 自定义颜色
orange <- "#EB736B"
blue <- "#57B5AB"
library(stringr)
# 1.4提取分期信息
b= meta$ajcc_pathologic_stage
a= str_extract_all(meta$ajcc_pathologic_stage,"I|V");head(a)#提取分期信息
b= sapply(a,paste,collapse = '')
head(b)
meta$stage = b
#去掉生存信息不全或者生存时间小于0.1(月)的病人，样本纳排标准不唯一，且差别很大
# k1 = meta$time>=0.1;table(k1)
# # 8
# k2 =!(is.na(meta$time)|is.na(meta$event));table(k2)
# 6
##只先看目前的数据
#注意生存分析只需要肿瘤的病人
head(normalized_counts)
# nm_tumor = normalized_counts[]
sample$sample_type_fix[sample$sample_type_fix=='Primary Tumor']='Tumor'
colnames(sample)
exp = normalized_counts[,sample$sample_type_fix %in% "Tumor"]
dim(exp)#412481 17 只存贮这17个
meta$case_submitter_id<- str_replace_all(meta$case_submitter_id ,"-", ".")
# meta$case_submitter_id =
# meta = meta[meta$case_submitter_id %in% str_sub(colnames(exp),1,12),]
meta = meta[meta$case_id %in% sample_tumor$case_id,]
sample_tumor = sample[sample$sample_type_fix %in% c("Tumor"),]
dim(sample_tumor) #34 215
dim(meta) #34 215
meta <- meta[!duplicated(meta$case_id), ]
dim(meta) #17 215
# exp = exp[,k1&k2]
# meta = meta[k1&k2,]
# rownames(meta)
# colnames(meta)
save(meta,exp,sample_tumor,file = './/surve_model_tumor.Rdata')
# load('./00.data//surve_model_tumor.Rdata')
## 提取数据
# rt <- meta
# colnames(rt)
## 生存分析
# rt$days_to_death =as.numeric(rt$days_to_death)
fit <- survfit(Surv(time,event) ~age_group, data = meta)
## 作图
ggsurvplot(fit,
pval = TRUE,
palette = c(orange, blue),
title="age_group",
xlab="Time")
orange <- "#EB736B"
blue <- "#57B5AB"
ggsurvplot(fit,
palette =c("#E7B800","#2E9FDF"),
risk.table =TRUE,pVal =TRUE,
conf.int =TRUE,xlab ="Time in months",
ggtheme =theme_light(),
ncensor.plot = TRUE)
ggsurvplot(fit,
palette =c("#EB736B","#57B5AB"),
risk.table =TRUE,pVal =TRUE,
conf.int =TRUE,xlab ="Time in months",
ggtheme =theme_light(),
ncensor.plot = TRUE)
#性别年龄
sfit1=survfit(Surv(time, event)~gender, data=meta)
sfit2=survfit(Surv(time, event)~age_group, data=meta)
splots<- list()
splots[[1]]<- ggsurvplot(sfit1,pval =TRUE,data = meta,risk.table = TRUE)
splots[[2]]<- ggsurvplot(sfit2,pval =TRUE,data = meta,risk.table = TRUE)
arrange_ggsurvplots(splots, print = TRUE,
ncol=2,nrow=1,risk.table.height=0.4)
dev.off()
#单个基因(选择对应的cg展示他们的风险)
g=rownames(exp)[1]
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
#多个基因
gs=rownames(exp)[1:4]
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=2,nrow=2,risk.table.height=0.4)
##log_rank_p
#加上KM生存曲线即可
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event == 1表示“1为结局时间发生“
log_rank_p <- pbapply(exp, 1, function(gene){
# gene <- exp[1,]
meta$group <- ifelse(gene > median(gene), "high", "low")
data.survdiff <- survdiff(mySurv ~ group, data = meta)
# survdiff()用于检验两条或多条生存曲线之间是否存在差异，
# 或者针对已知替代项检验单条曲线，也适用于分组资料以及多组间生存曲线的比较；
p.val <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) # 计算log-rank检验的p值
return(p.val)
})
log_rank_p <- sort(log_rank_p) # 取出每一个基因生存分析的P值，形成表
table(log_rank_p < 0.01) # 哪些是P<0.01的
table(log_rank_p < 0.05) # 哪些是P<0.05的
log_rank <- names(log_rank_p)[log_rank_p < 0.05] # log-rank检验筛选的差异基因 0.05
log_rank_p
# Cox比例风险回归（HR）
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event = 1表示“1为结局事件发生”
unicox_results <- pbapply(exp, 1, function(gene){ # 对表达矩阵exp中的每一行（基因）执行函数function
# gene = as.numeric(exp[1,]) # 举个例子
group <- ifelse(gene > median(gene), "high","low") # 标记不同样本中该基因的上下调，将表达量转换为二分类变量
survival_dat <- data.frame(group = group,
stage = meta$stage,
age = meta$age,
gender = meta$gender,
gene = gene, # 基因表达量
stringsAsFactors = F) # 选择用于生存分析的数据
# 拟合Cox比例风险回归模型 m ：
m <- coxph(mySurv ~ gender + age + stage + group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ gene, data = survival_dat) # 表达量也可直接使用连续型变量gene
# sm <- summary(m) # 看拟合出的Cox比例风险回归模型的参数
# n：总人数。
# number of events：发生结局事件的人数，即死亡人数
# coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
# exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
# se(coef)：coef的标准误
# z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
# Pr(>|z|)：p值
# low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
tmp <- round(cbind(coef = beta, # 模型的偏回归系数β
se = se, # 偏回归系数β的标准误
z = beta/se, # Wald检验的统计量
p = 1 - pchisq((beta/se)^2, 1), # p值，差异的统计学意义
HR = exp(beta), # 风险值，计算公式为exp(coef)
HRCILL = exp(beta - qnorm(.975, 0, 1) * se), # HR的95%置信区间低点
HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3) # HR的95%置信区间高点
return(tmp["grouplow",]) # 分类变量
# return(tmp["gene",]) # 连续型变量
})
unicox_results <- as.data.frame(t(unicox_results))
table(unicox_results[,4] < 0.01) # p值<0.01的有多少
table(unicox_results[,4] < 0.05) # p值<0.05的有多少
unicox <- rownames(unicox_results)[unicox_results[,4] < 0.05] # Cox单因素分析筛选的差异基因
save(log_rank_p,unicox_results,file = './/surve_cox_log_rank_p.Rdata')
#多个基因
# c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028"
gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
load('./00.data//surve_cox_log_rank_p.Rdata')
# 取交集
cd =intersect(log_rank,gs) #0
cd =intersect(unicox,gs) #0
gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
splots <- lapply(gs,function(g){
meta$gene = ifelse(exp[g,]> mean(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene,data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
})
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
#ggsave(paste0("../03.plot/", "volcano","_logFC_0.4_P.Value_10_5.pdf"), width = 6*1.2, height = 6*1)
##log_rank_p
# KM生存曲线即可
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event == 1表示“1为结局时间发生“
log_rank_p <- pbapply(exp, 1, function(gene){
# gene <- exp[1,]
meta$group <- ifelse(gene > mean(gene), "high", "low")
data.survdiff <- survdiff(mySurv ~ group, data = meta)
# survdiff()用于检验两条或多条生存曲线之间是否存在差异，
# 或者针对已知替代项检验单条曲线，也适用于分组资料以及多组间生存曲线的比较；
p.val <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) # 计算log-rank检验的p值
return(p.val)
})
log_rank_p <- sort(log_rank_p) # 取出每一个基因生存分析的P值，形成表
table(log_rank_p < 0.01) # 哪些是P<0.01的
table(log_rank_p < 0.05) # 哪些是P<0.05的
log_rank <- names(log_rank_p)[log_rank_p < 0.05] # log-rank检验筛选的差异基因 0.05
log_rank_p
# Cox比例风险回归（HR）
mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象，event = 1表示“1为结局事件发生”
unicox_results <- pbapply(exp, 1, function(gene){ # 对表达矩阵exp中的每一行（基因）执行函数function
# gene = as.numeric(exp[1,]) # 举个例子
group <- ifelse(gene > mean(gene), "high","low") # 标记不同样本中该基因的上下调，将表达量转换为二分类变量
survival_dat <- data.frame(group = group,
stage = meta$stage,
age = meta$age,
gender = meta$gender,
gene = gene, # 基因表达量
stringsAsFactors = F) # 选择用于生存分析的数据
# 拟合Cox比例风险回归模型 m ：
m <- coxph(mySurv ~ gender + age + stage + group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ group, data = survival_dat) # 表达量用二分类变量group
# m <- coxph(mySurv ~ gene, data = survival_dat) # 表达量也可直接使用连续型变量gene
# sm <- summary(m) # 看拟合出的Cox比例风险回归模型的参数
# n：总人数。
# number of events：发生结局事件的人数，即死亡人数
# coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
# exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
# se(coef)：coef的标准误
# z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
# Pr(>|z|)：p值
# low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
tmp <- round(cbind(coef = beta, # 模型的偏回归系数β
se = se, # 偏回归系数β的标准误
z = beta/se, # Wald检验的统计量
p = 1 - pchisq((beta/se)^2, 1), # p值，差异的统计学意义
HR = exp(beta), # 风险值，计算公式为exp(coef)
HRCILL = exp(beta - qnorm(.975, 0, 1) * se), # HR的95%置信区间低点
HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3) # HR的95%置信区间高点
return(tmp["grouplow",]) # 分类变量
# return(tmp["gene",]) # 连续型变量
})
unicox_results <- as.data.frame(t(unicox_results))
table(unicox_results[,4] < 0.01) # p值<0.01的有多少
table(unicox_results[,4] < 0.05) # p值<0.05的有多少
unicox <- rownames(unicox_results)[unicox_results[,4] < 0.05] # Cox单因素分析筛选的差异基因
save(log_rank_p,unicox_results,file = './00.data//surve_cox_log_rank_p_mean.Rdata')
cd =intersect(log_rank,gs) #0
cd =intersect(unicox,gs) #0
#输出均值的结果
# cg19548479
write.csv(log_rank_p,'./00.data//log_rank_p_mean.csv')
write.csv(unicox_results,'./00.data//unicox_results_mean.csv')
# cg19548479符合后面的几个检验
# gs=c("cg08875705","cg01979888", "cg01493517", "cg00983904", "cg23737737" ,"cg18222083", "cg19548479" ,"cg20170028")
g='cg19548479'
meta$gene = ifelse(exp[g,]> median(exp[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
arrange_ggsurvplots(splots, print = TRUE,ncol=4,nrow=2,risk.table.height=0.4)
最后出图结果如下：
Fig1样本质控
相关文章
更多 >