学徒带你7步3251行代码+300行注释完成TCGA数据库挖掘实战全文复现

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
 #0 准备R包,使用的是RTCGA包下载数据，我觉得最简单
#若已安装，则可以忽略
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(version = "3.10")
BiocManager::install("RTCGA")
BiocManager::install("RTCGA.miRNASeq")
BiocManager::install("RTCGA.clinical")
#加载包
#我们要的是乳腺癌的数据(TCGA中的BRCA)
library(RTCGA.miRNASeq)
?miRNASeq #使用帮助文档，来查看RTCGA.miRNASeq包中有哪些癌症的miRNA表达信息
#我们要探索BRCA.miRNASeq
 #1 下载数据
#1.1简单探索一下数据的结构
BRCA.miRNASeq[1:10,1:6]
#1.2获取表达矩阵
expr <- expressionsTCGA(BRCA.miRNASeq)
dim(expr) #查看一下维度
expr[1:10,1:6]
#探索结果发现
#1)从第3列开始才是miRNA的表达量
#2)每3行是一个样本，依次为这一样本的read_count,RPM,cross-mapped，我们要的是read_count
#3)获取的expr表达矩阵的行名没有了，且表达量都是字符串型
 #1.3整理表达矩阵，取出每个样本的read_count，把character转换为numeric
#从第3列开始取，取到最后一列,也就是3:ncol(expr)
expr1<-expr[seq(1,nrow(expr),by=3),3:ncol(expr)]
#[seq(1,nrow(expr),by=3),]表示间隔3行取，也就是取1、4、7、10...行，
dim(expr1)
expr1[1:6,1:6]
str(expr1) #查看结构也可以发现表达量都是character,要批量转换为numeric数值型
 expr1<-apply(expr1,2,as.numeric) #产生了NA
expr1[1:6,1:6]
#补上行名,不要忘记每三行是一个样本
rownames(expr1)=rownames(BRCA.miRNASeq)[seq(1,nrow(BRCA.miRNASeq),by=3)]
 expr1[1:6,1:6] #行名有了
 #2.数据过滤，过滤低表达量的miRNA
#有的miRNA在几百个样本中的表达量都为0，需要去除
 #2.1准备阶段：行列转换，转换成列名为样本名，行名为gene名，为后续进行批量处理做准备
expr1<-t(expr1)
dim(expr1)
#看一下有没有NA
table(apply(expr1,1,is.na)) #有1190个NA
Expr<-na.omit(expr1) #剔除有NA的行
table(apply(Expr,1,is.na)) #现在没有NA了
 #2.2过滤低表达量的miRNA
dim(Expr)
Expr<-Expr[apply(Expr, 1, function(x){
            sum(x>1)>10}),]
#sum()函数统计某一miRNA有多少个样本表达量>1,sum的结果>10，则符合
#apply()函数，对Expr这个数据框的行进行批量,然后返回的是每个miRNA在10个样本中的表达量是否>1的逻辑向量，TRUE还是FALSE
#只取那些逻辑值为True的miRNA的表达量，过滤掉FALSE的
dim(Expr)
#过滤掉369个miRNA
 #3 处理临床信息
#3.1获取临床信息
library(RTCGA.clinical)
BRCA_clinical<-BRCA.clinical
dim(BRCA_clinical)
View(BRCA_clinical)
#有3703列的临床信息，我们只需要根据自己所需来取就可以
#这篇文章其实只要生存状态和生存时间这2个主要临床信息
#患者的编码(样本号)是必须要的
write.csv(BRCA_clinical,quote = F,"BRCA_clinical_Data.csv")
#输出一下数据
 #3.2获取整理我们所需要的临床信息
#根据文章只需要patient.bcr_patient_barcode',
#'patient.vital_status',
#'patient.days_to_death',
#'patient.days_to_last_followup',这几列，我还选取了其他几列
#这里说明一下，TCGA中如果患者的生存状态为dead，其生存时间就是days_to_death计算
#如果为alive，生存时间就是days_to_last_followup
#我们只需要将days_to_last_followup和days_to_last_followup，就可以得到某一患者的生存时间
#但我发现有时候days_to_last_followup会存在负数和0，这时候我们要注意一下
#TCGA这里的时间是以天计算
BRCA_clinicaldata<-BRCA_clinical[c('patient.bcr_patient_barcode',
                                   'patient.vital_status',
                                   'patient.days_to_death',
                                   'patient.days_to_last_followup',
                                   'patient.race',
                                   'patient.age_at_initial_pathologic_diagnosis',
                                   'patient.gender',
                                   'patient.stage_event.pathologic_stage'
)]
str(BRCA_clinicaldata) #查看数据结构
head(BRCA_clinicaldata)
BRCA_clinicaldata[1:4,1:4] #现在行名不太对，要把patient.bcr_patient_barcode（患者的样本号）这一列作为行名
rownames(BRCA_clinicaldata) <- NULL #不设行名
head(BRCA_clinicaldata)
 #把patient.bcr_patient_barcode（患者的样本号）这一列作为行名
#tibble::column_to_rownames(data,var="")用到这个函数，把数据框中的指定列转换为行名
BRCA_clinicaldata<-tibble::column_to_rownames(BRCA_clinicaldata,var="patient.bcr_patient_barcode")
 head(BRCA_clinicaldata) #行名有了
head(Expr)#但需要大写，要和表达矩阵匹配
rownames(BRCA_clinicaldata)<-stringr::str_to_upper(rownames(BRCA_clinicaldata))
#stringr::str_to_upper()函数，把行名大写，和Expr一致
 head(Expr)
head(BRCA_clinicaldata)#但可以发现临床信息中的样本名只显示12个字符，而Expr表达矩阵里的样本名不止12个字符
dim(BRCA_clinicaldata)
dim(Expr)
Expr[1:4,1:4]
#同时，可以看到临床信息里是1098个样本；Expr表达矩阵里是1190个样本，677个miRNA信息（过滤后）
#后续需要对此进行处理
#保存数据
save(Expr,BRCA_clinicaldata,file="step1_TCGA_BRCA_miRNAexpr_clinical.Rdata")

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
 load("step1_TCGA_BRCA_miRNAexpr_clinical.Rdata")
#准备R包
#文章用的是edgeR包进行差异分析
BiocManager::install('edgeR')
#加载
library(edgeR)
#1 准备数据
#1.1设置样本分组
#TCGA样本编号的14-15位显示了该样本是癌症组还是正常组，01-09是Tumor,10-29是Normal
install.packages("stringr") #要用到这个包提取信息，没有的话安装一下，
library(stringr)
group_list<-ifelse(as.numeric(substr(colnames(Expr),14,15))<10, 'Tumor','Normal')
#substr(x,start,stop)函数，取指定字符位数区间，这里取14-15位
#使用ifelse(Condition,A,B)函数,如果满足这一Condition，则输出A，否则为B
#这里，如果满足14-15位小于10，则定义为'Tumor'，大于10，则为'Normal'
 group_list <- factor(group_list,levels = c("Normal","Tumor"))
group_list
table(group_list)
#1086个Tumor,104个Normal
#文章里是1103个Tumor,104个Normal,差别不大，大概是因为RTCGA包的数据不是最新的
 #1.2构建DGEList对象
dge<-DGEList(counts=Expr,group=group_list)
#DGEList 对象的主要组成是一个 counts 矩阵：每行对应一个基因，每列对应一个样本的数据，包含样本分组信息的数据
 #1.2过滤低表达的miRNA
#根据经验，如果 miRNA在所有条件下的所有样本都没有表达，那么就应该过滤掉。
#通常一个基因的 counts 为 5——10 才认为它表达了。
#edgeR包使用 count-per-million（CPM）指标过滤，考虑了文库大小
keep <- rowSums(cpm(dge)>1) >= 3
#CPM 为 1 代表在最小的样本中，count 为 6——7。
#这里所求至少在 3 个库中的CPM表达量>1（>=3）
table(keep) #有12个miRNA不符合
dge <- dge[keep, , keep.lib.sizes=FALSE] #进行过滤，重置文库大小
 #edgeR包建议过滤后重新计算库的大小，尽管它的影响非常小
dge$samples$lib.size <- colSums(dge$counts)
 #1.3对数据进行标准化
#去除文库大小差异，解决测序深度不同的问题导致read count差异的问题
#edgeR 不需要考虑基因长度的影响，它认为每个基因对每个样本的 read counts 的影响是一样的。
#edgeR考虑的是文库大小， 通过不同库的大小修正测序深度对差异表达分析的影响
 #edgeR默认采用TMM归一化法，计算每个样本的 sizefactor
dge <- calcNormFactors(dge, method = 'TMM')
dge$samples
 #后续需要对筛选出的差异表达的miRNA的表达量值进行绘制热图和PCA图。
#edgeR建议使用logCPM这个指标来画图
#To draw a heatmap of individual RNA-seq samples, we suggest using moderated log-counts-per-million.
logCPM <- cpm(dge, log=TRUE, prior.count=2)
#prior.count=2 也就是logCPM=log2(CPM+2/L),避免得到的logCPM为0
 #1.4 构建实验矩阵，为估计离散值做准备
design<-model.matrix(~0+group_list)
colnames(design)<-levels(group_list)
rownames(design)<-colnames(dge)
 #2开始差异分析
#2.1估计离散值
#对于多因素实验，edgeR 方法计算离散度，它通过一个设计矩阵拟合 GLM 模型来考虑多因素影响。
#一次性计算 common dispersion、trended dispersion 和 tagwise dispersion 的命令：
dge <-estimateDisp(dge,design) #使用这个代替下面三步
#dge <- estimateGLMCommonDisp(dge,design)
#dge <- estimateGLMTrendedDisp(dge, design)
#dge <- estimateGLMTagwiseDisp(dge, design)

#2.2 进行差异基因分析
###  使用likelihood ratio test (似然比检验)
#对应的是glmFit()和glmLRT()；
#说明书推荐用 quasi-likelihood (QL) F-test
#对应的是glmQLFit()和glmQLFTest()
fit <- glmFit(dge, design) #拟合模型
fit2 <- glmLRT(fit, contrast=c(-1,1))  #统计检验
#我这里只有Tumor和Normal组，所以只需要contrast=c(-1,1)；如果进行多组，可以查看说明书

#2.3展示差异分析结果
allDEM<-topTags(fit2, n=nrow(dge))
allDEM<-as.data.frame(allDEM)
head(allDEM)
write.csv(allDEM,file="allDEM.csv") #输出结果
 #3 设定差异基因的阈值
#文章是FDR<0.05,fold change(FC)>2，也就是logFC>1
FDR_cutoff<-0.05
logFC_cutoff<-1
allDEM$change<-as.factor(ifelse(allDEM$FDR>FDR_cutoff,'No', #如果FDR>0.05，则定义为'No'，否则进行下一步的ifelese
                                ifelse(allDEM$logFC>logFC_cutoff,'Up', #在FDR<=0.05前提下，logFC>0.05则定义'Up',否则进行下一步的ifelse
                                       ifelse(allDEM$logFC< -logFC_cutoff,'Down','No'))))
                                    #如果FDR<0.05,且logFC<1情况下，logFC<-1则定义为'Down',否则定义为'No'
 table(allDEM$change)
#文章共筛选出370个差异表达 miRNAs，其中有 108 个下调的miRNAs，262 个 上调的miRNAs
#我们这里共筛选出269个差异表达miRNAs,其中64个下调，204个上调
 #4 把过滤好的DEG提取出来
#4.1提取出所有的上调和下调miRNA
SigDEMname<-rownames(allDEM)[allDEM$change !="No"]
#这里只需要选择change不等于No的就可以了，选出的就是Up和Down
SigDEM<-allDEM[SigDEMname,]
write.csv(SigDEM,file="SigDEM.csv") #输出数据

#4,2只提取出所有的上调miRNA,这里后续要用到
Up_DEMname<-rownames(allDEM)[allDEM$change=="Up"]
Up_DEM<-allDEM[Up_DEMname,]
write.csv(Up_DEM,file="Up_DEM.csv") #输出数据
 #4,3只提取出所有的下调miRNA
Down_DEMname<-rownames(allDEM)[allDEM$change=="Down"]
Down_DEM<-allDEM[Down_DEMname,]
write.csv(Down_DEM,file="Down_DEM.csv") #输出数据
 #5 volcano plot 绘miRNA差异分析结果的的火山图
#文章没有绘制火山图，这里绘制一下
#5.1检查数据
head(allDEM)
#发现hsa-mir-133a-1的logFC<-6，应该在Tumor中是低表达的
#我们检查一下
Expr['hsa-mir-133a-1',]
#验证了这一结果，说明分组正确
install.packages("ggplot2")
library(ggplot2)
volcano_plot <- ggplot(data = allDEM,
            aes(x = logFC,
                y = -log10(FDR))) +
  geom_point(alpha=0.4, size=3.5,
             aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="grey",lwd=0.8) +
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.8) +
  theme_bw()+ labs(title="Volcano Plot",x=expression(log[2](FC)), y=expression(-log[10](FDR)))
#这里根据筛选差异miRNA的筛选阈值，绘制了渐近线
#labs(x=expression(log[2](FC)), y=expression(-log[10](FDR)))使用这个命令，可以让log的书写更加美丽
volcano_plot
ggsave("volcano plot.tiff") #保存图片
#5.2 绘制DEMs for heatmap
#5.2.1这里用的是gplots包中的heatmap.2，在生信20题见到过有印象，一般性我都用pheatmap
install.packages("gplots")
library(gplots)
?heatmap.2
#准备数据
head(SigDEM)
#这里用到前面准备好的logCPM来绘制热图（说明书建议使用logCPM来绘制热图和PCA图）
logCPM[1:4,1:4]
logCPM_SigDEM_expr<-logCPM[rownames(SigDEM),] #制作只含SigDEM的表达矩阵
head(logCPM_SigDEM_expr)
logCPM_SigDEM_expr[1:4,1:4]
 #对绘制的热图设定上下限值
n<-logCPM_SigDEM_expr
n[1:4,1:4]
n[n>2]=2 #设置上限，logCPM>2的都等于2
n[n< -2]= -2 #设置上限，logCPM<-2的都等于-2
 #发现绘制heatmap.2的数据格式必须是数值矩阵
class(n) #是matrix
?heatmap.2
heatmap.2(n,col='greenred',  #颜色可以自己选，文章用的红绿颜色
          labCol = NA, #不显示列名
          trace = "none" #一般性都是trace = "none"
)
#这里我选择手动Export导出，可以自己挪动尺寸
 #保存数据
save(Expr,SigDEM,logCPM_SigDEM_expr,Up_DEM,Down_DEM,group_list,BRCA_clinicaldata,file="step2_DEMdata.Rdata")

#载入数据
rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
load("step2_DEMdata.Rdata")
 #1安装R包
install.packages("survival")
install.packages("survminer")
install.packages("stringr")
 ##2 整理表达矩阵和临床信息
#文章使用筛选出的上调的DEMs做Cox生存分析，不需要Normal组了，只需要Tumor组
#2.1 提取上调的miRNA的表达量数据框,如果是对所有的上调和下调做Cox。那就改一下命令
#这里可以使用edgeR包中TMM标准化后的logCPM进行后续的临床分析
#也可以使用log2(x+1)对表达量进行归一化，文章似乎用的这种方法（文章没有具体说明）
#也可以使用RPKM,FPKM....
#这里我展示的是使用log2(x+1),和文章结果更相近
#但logCPM绘制的PCA效果区分明显（毕竟edgeR包推荐用logCPM画PCA），log2(X+1)的PCA不好
 #这里使用log2(x+1)进行COX回归分析
Up_DEM_expr<-Expr[rownames(Up_DEM),]  #制作只含Up_DEM的表达矩阵
Up_DEM_expr[1:4,1:4]
Up_DEM_expr<-as.data.frame(Up_DEM_expr)
Expr[1:4,1:4]
Up_DEM_expr<-log2(Up_DEM_expr+1)
 #如果使用edgeR得到的logCPM做COX回归分析
#Up_DEM_expr<-logCPM_SigDEM_expr[rownames(Up_DEM),] #制作只含Up_DEM的表达矩阵
#检查一下数据
Up_DEM_expr[1:4,1:4]
dim(Up_DEM_expr)
dim(BRCA_clinicaldata)
table(group_list)
 ###插曲，这里如果你要画PCA的画,这里推荐用logCPM绘制热图，也是edgeR包推荐的
install.packages("FactoMineR")
install.packages("factoextra")
library("FactoMineR")#进行PCA分析
library("factoextra")#PCA可视化
 Exprset_PCA=t(Up_DEM_expr)#画PCA图时要求是行名是样本名，列名是探针名，因此此时需要行列转换
Exprset_PCA<-as.data.frame(Exprset_PCA)#PCA图要求是data.frame，因此将matrix转换为data.frame
Exprset_PCA.pca <- PCA(Exprset_PCA, graph = F)
Exprset_PCA=cbind(Exprset_PCA,group_list) #将样本分组信息追加到最后一列
 #PCA可视化，fviz_pca_ind按样本  fviz_pca_var按基因
fviz_pca_ind(Exprset_PCA.pca,
             title = "Principal Component Analysis",
             geom.ind = "point", ###  show points only (nbut not "text"),这里c("point", "text)2选1
             col.ind = Exprset_PCA$group_list, ###  color by groups
             ###  palette = c("#00AFBB", "#E7B800"),自定义颜色
             addEllipses = TRUE, ###  Concentration ellipses加圆圈
             legend.title = "Groups")

#2.2取出Tumor组的表达矩阵(做生存分析，不需要Normal样本，对差异的上调的miRNA的表达矩阵只取Tumor的）
Up_DEM_expr<-Up_DEM_expr[,group_list=='Tumor']
dim(Up_DEM_expr)
dim(BRCA_clinicaldata)  #可以发现临床信息中有1098个样本，表达矩阵中只有1086个，后续要进行匹配一致
 #2.3 整理临床信息
#2.3.1 为临床信息重新命名
str(BRCA_clinicaldata) #查看一下临床信息的结构，发现有些数据的类型需要进行转换
colnames(BRCA_clinicaldata)<-c("event","days_to_death","days_to_last_followup","race","age","gender","stage")
 #2.3.2 使得表达矩阵和临床信息的样本一致
head(BRCA_clinicaldata)
head(Up_DEM_expr)
#让临床信息和表达矩阵的样本名的1-12个字符匹配
library(stringr)
BRCA_clinicaldata<-BRCA_clinicaldata[match(substr(colnames(Up_DEM_expr),1,12),rownames(BRCA_clinicaldata)),]
#使用到match函数，match[A,B]，A是被匹配的，B是要匹配的，返回的是Up_DEM_expr的列名（样本名）在BRCA_clinicaldata行名中的位置
 dim(Up_DEM_expr)
dim(BRCA_clinicaldata)
head(colnames(Up_DEM_expr))
head(BRCA_clinicaldata)
#现在数目对了
#检查一下
all(substr(rownames(BRCA_clinicaldata),1,12)==substr(colnames(Up_DEM_expr),1,12)) #现在正确了
rownames(BRCA_clinicaldata)<-colnames(Up_DEM_expr)
 #2.4整理生存分析Input数据，进行生存分析需要生存状态和生存时间2个数据
#2.4.1计算生存时间(Dead按照days_to_death计算，Alive按照days_to_last_followup计算，要注意如果最后得到的生存时间可能有0和负数，要删去)
str(BRCA_clinicaldata)
#把NA都赋值为0
BRCA_clinicaldata[,2][is.na(BRCA_clinicaldata[,2])]<-0
BRCA_clinicaldata[,3][is.na(BRCA_clinicaldata[,3])]<-0
str(BRCA_clinicaldata)
#2.4.1.1计算生存时间,按照天
BRCA_clinicaldata$time_days<-as.numeric(BRCA_clinicaldata[,2])+as.numeric(BRCA_clinicaldata[,3])
table(BRCA_clinicaldata$time_days<= 0) #有53个数据的生存时间<=0要删去
BRCA_clinicaldata<-BRCA_clinicaldata[BRCA_clinicaldata$time_days>0,]   #删去生存时间为<=0的数据
table(BRCA_clinicaldata$time_days<= 0)
#2.4.1.2生存时间以年记录
BRCA_clinicaldata$time_year<-BRCA_clinicaldata$time_days/365
 #2.4.2根据dead和alive定义生存状态event，alive=0，dead=1
BRCA_clinicaldata$event<-ifelse(BRCA_clinicaldata$event=='alive',0,1)
str(BRCA_clinicaldata)
table(BRCA_clinicaldata$event) #剔除<=0的生存时间后，一共有928个alive,105个dead
 #2.4.3要再把表达矩阵和临床信息匹配，因为删去了生存<=0的数据
dim(BRCA_clinicaldata)
dim(Up_DEM_expr)
Up_DEM_expr<-Up_DEM_expr[,rownames(BRCA_clinicaldata)]
dim(Up_DEM_expr)
all(colnames(Up_DEM_expr)==rownames(BRCA_clinicaldata))
 #3进行单因素COX回归分析
#3.1首先做一个基因的生存分析
#COX回归的教程：http://www.sthda.com/english/wiki/cox-proportional-hazards-model
#COX回归分析的格式：
#coxph(Surv(time, status) ~ sex, data = lung)
###  coxph(formula, data, method)
 #3.1.1构建生存对象
#Surv(time, event)
library(survival)
Surv1<-with(BRCA_clinicaldata,Surv(time_year,event))
#我们的生存时间和生存状态储存在BRCA_clinicaldata中
 Up_DEM_expr<-t(as.data.frame(Up_DEM_expr)) #将表达矩阵行列转换，变成行为样本名，列为基因
dim(Up_DEM_expr)
Up_DEM_expr<-as.data.frame(Up_DEM_expr)
res.cox<-coxph(Surv1~Up_DEM_expr[,1],data=Up_DEM_expr) #对第一个miRNA做COX回归
res.cox #查看分析结果
res.cox_summary<-summary(res.cox) #使用summary()函数可以查看更完整关于COX模型的报告
res.cox_summary
 #3.2对所有miRNA进行COX回归分析，开始批量分析啦
#3.2.1 创建生存对象
head(BRCA_clinicaldata)
?Surv
mySurv<-with(BRCA_clinicaldata,Surv(time_year,event))
class(Up_DEM_expr)
head(Up_DEM_expr)
#3.2.2批量分析每个基因的表达量与生存的关系
univ_models <-apply(Up_DEM_expr,2,function(gene){
  res_cox<-coxph(mySurv ~  gene, data =  Up_DEM_expr)
})

#3.2.3把得到的COX结果用lappy()批量整理成文章需要的表格形式
options(scipen=200) #在200个数字以内都不使用科学计数法，如果不设置这步，可能会对后续按照P排序造成Bug
univcox_results<-lapply(univ_models,function(x){
  x<-summary(x) #后续从中提取信息
  P.value<-round(x$wald["pvalue"], digits=4) #P-value #round(x,digits=4),将x保留4位小数
  wald.test<-round(x$wald["test"], digits=4) #wald.test
  beta<-round(x$coef[1],digits=4)#coef beta
  HR<-round(x$coef[2],digits=4) #exp(coef)，也就是HR值
  HR_CI_lower<-round(x$conf.int[,"lower .95"],4)
  HR_CI_upper<-round(x$conf.int[,"upper .95"],4)
  HR<-paste0(HR,"(",HR_CI_lower,"-",HR_CI_upper,")")
  z<-round(x$coef[4],digits=4) #z值，z = coef/se(coef)
  res<-c(beta, HR, z,wald.test, P.value)
  names(res)<-c("beta","HR (95% CI)","z","wald.test","P-value")
  return(res)
})
 results <- t(as.data.frame(univcox_results, check.names = FALSE)) #整理表格形式
class(results)
results<-as.data.frame(results)
View(results)
head(results) #整理成我们需要放在文章的表格形式了
 #把P-value按照升序排序，为后期按照P筛选做准备
results<-results[order(results$`P-value`,decreasing = F),]
View(results)
head(rownames(results),21) #查看一下大致上和文章的结果一样吗，好像差不多
#输出结果
save(results,Expr,SigDEM,Up_DEM_expr,BRCA_clinicaldata,file="step3-univariate2_log2(x+1) COX_analysis.Rdata")
write.csv(results,file="univariate COX analysis_log2(x+1).csv") #保存表格

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
#多因素COX回归分析（逐步回归）
#加载数据
load("step3-univariate2_log2(x+1) COX_analysis.Rdata")
#1准备R包
library(survival)
 #2准备数据
#2.1筛选P<0.05的上调miRNA，去做多因素COX回归（文章筛选的是P<0.05，我这里也筛选P<0.05的）
Up_univariate.res<-results[results$`P-value`< 0.05,] #P<0.05为阈值
dim(Up_univariate.res) #一共30个
write.csv(Up_univariate.res,file="Up_univariate.res.csv") #输出数据
#2.2准备多因素COX的Input数据
Up_DEM_expr[1:4,1:4]
Up_univariate.res[1:4,1:4]
multi_input.expr<-Up_DEM_expr[,rownames(Up_univariate.res)] #提取出30个miRNA的表达量数据
dim(multi_input.expr)
multi_input.expr[1:4,1:4]
#3进行多因素COX回归分析(逐步回归)
str(BRCA_clinicaldata)
#3.1首先产生生存对象
mySurv<-with(BRCA_clinicaldata,Surv(time_year,event))
#3.2开始多因素COX分析
multi_COX<-coxph(mySurv ~ ., data=multi_input.expr)  
###  "~ ." 这个的"."表示把所有变量都放进去进行多因素分析
#如果你只想分析hsa-mir-181c和hsa-mir-466，可以这样写
#multi_COX<-coxph(mySurv ~ `hsa-mir-181c` + `hsa-mir-466`, data=multi_input.expr)
###  ~后面用"+"把你要分析的变量放进去
 summary(multi_COX)#看一下结果COX回归分析的结果
 #文中使用step法（逐步回归）建立COX模型
#逐步回归可以进一步筛选自变量，剔除对因变量不重要或者不显著的自变量，得到最优的回归模型
?step #查看一下说明书
step.multi_COX=step(multi_COX,direction = "both")  #方法可以选择"both"（两边，这是默认的方法）,"backward"（向后逐步回归）,"forward"（向前逐步回归）
step.multi_COX #查看结果,最后筛选出10个
 #根据公式，h(t)=h0(t)×exp(b1x1+b2x2+...+bpxp)
#t代表生存时间
#h(t)代表t时刻死亡的风险，由x1、x2、x3这些变量组成来预测这个风险
#b1 、b2等就是自变量的x1、x2的coef(beta值)，即回归系数
#h0代表基线风险，也就是x都为0的生活，ho(t)就是1，因为exp(0)=1
 step.multi_COX
#coef,即回归系数b(beta)
#exp(coef)，即HR，风险比(HR-hazard ratio)
#如果这里P<0.05则代表在纳入了其他自变量的情况下，仍然显著，是独立预测因子
#HR = 1: No effect
#HR < 1: Reduction in the hazard
#HR > 1: Increase in Hazard
#对于癌症研究，HR>1: bad prognostic factor; HR<1:good prognostic factor
#这里拿'hsa-mir-301a'举例，P>0.05，HR=1.2
#说明，保持其他自变量不变的情况下，随着hsa-mir-301a的表达量的升高，患者死亡的风险会增高，
#但是它在多因素分析中P>0.05，原本在单因素中P<0.05，说明在纳入其他变量后，受到其他变量影响，贡献不是那么大,故不是独立预测因子
 #这里写成COX风险比例的方程就是：
#SRS(Survival Risk Score)=(-0.33743hsa-mir-181c)+(0.44684hsa-mir-466)+(0.32558hsa-mir-548f-1)+
#(-0.45569hsa-mir-106a)+(0.19587hsa-mir-3200)+(-0.34022hsa-mir-449c)+(-0.15158hsa-mir-618)+
#(0.31046hsa-mir-549)+(0.38831hsa-mir-3927)+(0.18810hsa-mir-301a)
#hsa-mir-106a这里自变量都是指的是miRNA的表达量
 #4风险评分计算
#4.1使用Survival程序包的Predict函数，计算出每位患者的风险评
RiskScore<-predict(step.multi_COX,type = "risk",newdata =multi_input.expr)
RiskScore #可以看到每个患者都有一个对应的风险得分
 #4.2，以风险评分的中位值将患者分为两组，大于中位值的 患者为高风险组，小于或等于中位值的患者为低风 险组
risk_group<-ifelse(RiskScore>=median(RiskScore),'high','low')
risk_group
#5保存数据
save(BRCA_clinicaldata,Expr,Up_DEM_expr,step.multi_COX,RiskScore,risk_group,file="step4-Multivariate Cox regression analysis.Rdata")

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
 #miRNA heatmap 文章的热图是low都在一边，high都在一边
 #载入数据
load("step4-Multivariate Cox regression analysis.Rdata")
 #1准备R包
install.packages("dplyr")
install.packages("pheatmap")
#2准备输入数据
#2.1提取出这10个miRNA的名字
step.multi_COX
ten_miRNA_names<-c("hsa-mir-181c","hsa-mir-466","hsa-mir-548f-1","hsa-mir-106a",
                   "hsa-mir-3200","hsa-mir-449c","hsa-mir-618",
                   "hsa-mir-549","hsa-mir-3927","hsa-mir-301a")
#2.2取出这10个miRNA的log2(x+1)表达量和准备分组注释
Up_DEM_expr[1:4,1:4]
miRNA_heatmap_input<-Up_DEM_expr[,ten_miRNA_names]
miRNA_heatmap_input[1:4,1:4]
risk_group<-as.data.frame(risk_group) #对样本的分组注释
#检查样本名字是否对上
all(rownames(risk_group)==rownames(miRNA_heatmap_input))#正确
heatmap.dat<-cbind(id=rownames(miRNA_heatmap_input),risk_group,miRNA_heatmap_input)
heatmap.dat$risk_group<-factor(heatmap.dat$risk_group,levels=c("low","high"),ordered = TRUE)
#这一步很重要，把风险组按照，low和high排序，ordered=TRUE很重要
heatmap.dat$risk_group
library(dplyr)
heatmap.dat<-arrange(heatmap.dat,risk_group) #使用arrange()函数，按照risk_group排序，生信星球学到的
heatmap.dat[1:4,1:4]#发现没有行名了
rownames(heatmap.dat)<-heatmap.dat[,1] #补上行名
##分组注释信息
risk_group<-heatmap.dat[,2] #获取第2列的risk_group
risk_group<-as.data.frame(risk_group)
rownames(risk_group)<-rownames(heatmap.dat)
colnames(risk_group)<-c("risk group")
#表达量矩阵
miRNA_heatmap_input<-heatmap.dat[,3:ncol(heatmap.dat)] #第3列开始是表达量
 #3绘制热图，这里用的是pheatmap绘制的热图，第一次的热图使用的heatmap.2
library(pheatmap)
miRNA_heatmap_input[1:4,1:4]
#需要行名位基因名，列名位样本的数据框
miRNA_heatmap_input<-t(miRNA_heatmap_input)
 ten_miRNA_heatmap<-pheatmap(miRNA_heatmap_input, #准备的数据
              color = colorRampPalette(c("green", "black", "red"))(50), #这里用这个颜色和文章颜色差不多
            annotation_col = risk_group, #分组注释
          show_colnames =FALSE, #不显示列名
         show_rownames = TRUE, #显示行名
         cluster_cols=FALSE,  #这里不对列聚类，因为我们已经自己分好组了
         cluster_rows=TRUE)   #对行聚类
ten_miRNA_heatmap
 save(ten_miRNA_heatmap,file="step5_Ten-miRNA_heatmap.Rdata")

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
 #KM生存曲线
#http://www.sthda.com/english/wiki/survival-analysis-basics
#这里很详细的教程
#载入数据
load("step4-Multivariate Cox regression analysis.Rdata")
 #载入R包
install.packages("stringr")
install.packages("survival") #用于生存分析
install.packages("survminer") #用于画图
 #1.准备KM的Input数据
#整理成含有time,event,risk_group(高低风险组)的数据
library(stringr)
str(BRCA_clinicaldata)
dim(BRCA_clinicaldata)
class(names(RiskScore))
#检查一下样本名
all(substr(names(RiskScore),1,12)==substr(rownames(BRCA_clinicaldata),1,12))
all(substr(names(risk_group),1,12)==substr(rownames(BRCA_clinicaldata),1,12))
KM.input<-cbind(BRCA_clinicaldata[,c("event","time_year")],RiskScore,risk_group) #用到cbind()函数
#3.进行KM生存分析
library(survival)
library(survminer)
#3.1计算生存曲线：survfit()
str(KM.input)
fit<-survfit(Surv(time_year,event) ~ risk_group, data=KM.input)
###  ~risk_group 表示通过高低风险组来计算患者的生存率 如果是按照性别那就 ~sex（sex是你的变量名）
 print(fit) #查看一下结果
summary(fit) #展示更详细的结果
summary(fit)$table #这是以table形式展示
 #3.2进行可视化  我使用的是ggsurvplot()这个函数 [in Survminer R package]
KMsurvival_plot<-ggsurvplot(fit,pval = TRUE, #show p-value of log-rank test，显示log-rank分析得到的P值
           conf.int = TRUE, #添加置信区间
           conf.int.style = "step",  ###  customize style of confidence intervals,改变置信区间的样子
           risk.table = "abs_pct",  ###  absolute number and percentage at risk，这里以n(%)的形式展示risk table
           risk.table.y.text.col = T,###  colour risk table text annotations.
           risk.table.y.text = FALSE,###  show bars instead of names in text annotations in legend of risk table.不显示注释名字
           xlab = "Time in years",   ###  customize X axis label.自定义x的标签为time in years
           surv.median.line = "hv", #添加中位生存时间的线
           ncensor.plot = FALSE, #我这里不显示删失的图，TRUE就显示
           legend.labs =
             c("high risk", "low risk"),    ###  对legend的标签重新命名
           palette = c("#E7B800", "#2E9FDF"), ###  自定义颜色
           ggtheme = theme_light()#绘图主题
               )
KMsurvival_plot
 #3.3 生存曲线的总结，Kaplan-Meier life table: summary of survival curves
#这个更为详细
KMres.sum  <- surv_summary(fit)
head(KMres.sum)
attr(KMres.sum , "table") #获取表格形式
 #3.3查看统计学结果  Log-Rank test comparing survival curves: survdiff()
#The log-rank test is the most widely used method of comparing two or more survival curves.
#前面我们在画图中也可以直接看到P值
surv_diff <- survdiff(Surv(time_year, event) ~ risk_group, data = KM.input)
surv_diff
#p-values<0.05 说明高低风险组的生存概率有显著差异
 #保存数据
save(KM.input,KMres.sum,BRCA_clinicaldata,KMsurvival_plot,step.multi_COX,Up_DEM_expr,risk_group,RiskScore,file="step6_Kaplan-Meier Survival Analysis.Rdata")

rm(list = ls())  ####  魔幻操作，一键清空~
options(stringsAsFactors = F)#在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式
 #时间依赖的ROC曲线
#载入数据
load("step6_Kaplan-Meier Survival Analysis.Rdata")
 #1.载入R包
install.packages("survivalROC")
install.packages("timeROC") #2个包都可以绘制生存时间依赖的ROC曲线
#1.使用SurvivalROC （不能显示置信区间和SD）
library(survivalROC)
?survivalROC #查看这个函数的格式
#输入数据
Survival_ROC_input<-KM.input
#这里我预测五年的生存率
survival_ROC<-survivalROC(Stime=Survival_ROC_input$time_year, #生存时间，Event time or censoring time for subjects
                          status=Survival_ROC_input$event, #生存状态,dead or alive
                          marker=Survival_ROC_input$RiskScore, #风险得分，Predictor or marker value
                          predict.time=5, #预测5年的生存时间
                          method="KM" #使用KM法进行拟合，默认的方法是method="NNE"
                          )
survival_ROC_AUC<-round(survival_ROC$AUC,3)#ROC曲线的AUC保留3位小数（文章保留了3位）
 #画图
#x轴为False positive rate，y轴为True positive rate
plot(survival_ROC$FP,survival_ROC$TP,type="l",xlim=c(0,1),ylim=c(0,1),
     xlab="False positive rate",  
     ylab="True positive rate",
     main=paste0("ROC Curve", " (", "AUC = ",survival_ROC_AUC," )"),  #标题
     cex.main=1.5,#放大标题
     cex.lab=1.3,#坐标轴标签（名称）的缩放倍数
     cex.axis=1.3, font=1.2, #放大轴上的数字
     lwd=1.5, #指定线条宽度
     col="red"  #红色
)
abline(a=0,b=1) #添加一条y=x的线
 #2.使用timeROC(可以计算置信区间和SD）
#Time ROC可以时计算多个时间的AUC
#文章没有说明计算的几年的，我这里计算3，5，10年
library(timeROC)
?timeROC #看一下说明书
#输入数据
time_ROC_input<-KM.input
time_ROC<-timeROC(T=time_ROC_input$time_year, #生存时间(dead和alive的生存时间).
                  delta=time_ROC_input$event, #生存结局，Censored的样本必须用0表示
                  marker=time_ROC_input$RiskScore, #预测的变量，这里是风险评分，在没有特殊情况下，值越大，越容易发生危险事件
                  cause=1, #阳性结局的赋值（必须是1或2），也就是dead的赋值，这里dead是1表示的
                  weighting = "marginal", #权重计算方法，这是默认方法，选择KM计算删失分布，weighting="aalen" [选用COX]，weighting="aalen" [选用Aalen]
                  times = c(3,5,10), #计算3、5、10年的ROC曲线
                  ROC=TRUE,
                  iid=TRUE #计算AUC
                 )
time_ROC #查看结果，可以看到，还包括了SE
#绘制ROC曲线啦
summary(time_ROC) #返回12个参数
time_ROC$AUC
#绘制3年的ROC
plot(time_ROC,time=3,col="red",title=FALSE,lwd=2) #将生成一条两倍于 默认宽度的线条
#在此基础上添加5年的ROC
plot(time_ROC,time=5,col="blue",add=TRUE,title=FALSE,lwd=2)
#add 10年的ROC
plot(time_ROC,time=10,col="black",add=TRUE,title=FALSE,lwd=2)
#添加图例
?legend
legend("bottomright", #图例设置在右下角
       c(paste0("AUC at 3 years = ", round(time_ROC$AUC[1],3)),
         paste0("AUC at 5 years = ", round(time_ROC$AUC[2],3)),
         paste0("AUC at 10 years = ", round(time_ROC$AUC[3],3))),
       col=c("red","blue","black"),lwd=2,bty="n" #或者bty+“o"
  )
 #使用ggplot2来画图
library(ggplot2)
time_ROC$TP
summary(time_ROC)
#整理成数据框，gpplot需要数据框
time_ROC.res<-data.frame(TP_3year=time_ROC$TP[,1], #获取3年的ROC的TP
                         FP_3year=time_ROC$FP[,1],  #获取3年的ROC的FP
                         TP_5year=time_ROC$TP[,2],  #获取5年的ROC的TP
                         FP_5year=time_ROC$FP[,2], #获取5年的ROC的FP
                         TP_10year=time_ROC$TP[,3], #获取10年的ROC的TP
                         FP_10year=time_ROC$FP[,3]) #获取10年的ROC的FP
time_ROC$AUC #这里放了3，5，10年的AUC，后续分别取子集time_ROC$AUC[[1]]，time_ROC$AUC[[2]]，time_ROC$AUC[[3]]
TimeROC_plot<-ggplot()+
  geom_line(data=time_ROC.res,aes(x=FP_3year,y=TP_3year),size=1,color="red")+
  geom_line(data=time_ROC.res,aes(x=FP_5year,y=TP_5year),size=1,color="blue")+
  geom_line(data=time_ROC.res,aes(x=FP_10year,y=TP_10year),size=1,color="black")+
  geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey",size=1, linetype = 2 #添加虚线
            )+
  theme_bw()+
  annotate("text",x=0.75,y=0.25,size=4.5,
           label=paste0("AUC at 3 years = ",round(time_ROC$AUC[[1]],3)),color="red")+
  annotate("text",x=0.75,y=0.15,size=4.5,
           label=paste0("AUC at 5 years = ",round(time_ROC$AUC[[2]],3)),color="blue")+
  annotate("text",x=0.75,y=0.05,size=4.5,
           label=paste0("AUC at 10 years = ",round(time_ROC$AUC[[3]],3)),color="black")+
  labs(x="False positive rate",y="True positive rate")+
  theme(axis.text=element_text(face="bold", size=11,  color="black"),#加粗刻度标签
        axis.title=element_text(face="bold", size=14, color="black")) #加粗xy轴标签名字
TimeROC_plot  #非常美丽