4.1.Survival Analysis
1) Background
1958年，Edward L. Kaplan 和Paul Meier也首次在临床研究中提出了生存曲线的概念，又被称作Kaplan-Meier曲线，主要用来对各组患者的生存状况进行描述。绘制生存曲线最主要的目的是进行生存分析，即通过将终点事件和出现这一终点所经历的时间结合起来进行统计分析，从而对两组患者的预后进行比较。
2) Pipeline
Survival analysis is a branch of statistics for analyzing the expected duration of time until one or more events happen, such as death in biological organisms and failure in mechanical systems. (https://en.wikipedia.org/wiki/Survival_analysis)
生存分析是研究生存时间和相关因素有无关系以及样本生存时间的分布规律的一种统计分析方法。
生存分析使用的方法：
描述生存过程：Kaplan-Meier plots to visualize survival curves（根据生存时间分布，估计生存率及其标准误，绘制生存曲线。常用Kaplan-Meier法，还有寿命法）
比较生存过程：Log-rank test to compare the survival curves of two or more groups（通过比较两组或者多组之间的的生存曲线，一般是生存率及其标准误，从而研究之间的差异，一般用log rank检验)
影响生存时间的因素分析：Cox proportional hazards regression to describe the effect of variables on survival（用Cox风险比例模型来分析变量对生存的影响，可以两个及两个以上的因素）
Reference: http://www.sthda.com/english/wiki/cox-proportional-hazards-model
3) Data structure
File name | Description |
rna.rds | RSEM normalized counts value matrix |
clinical_info.rds, LIHC.merged_only_clinical_clin_format.txt, all_clin.rds | Clinical information for TCGA samples |
3a) Input data
cd /home/test/clinical_analysisImport data
setwd("/home/test/clinical_analysis")
rna = readRDS(file="/home/test/clinical_analysis/rna.rds")
#if you don't want to normalize and scale expression data, just load data:
#exp = readRDS(file="/home/test/clinical_analysis/exp.rds")
clinical_info = readRDS(file="/home/test/clinical_analysis/clinical_info.rds")Data character
rna[1:4,1:4]
# TCGA.2V.A95S.01A.11R.A37K.07 TCGA.2Y.A9GS.01A.12R.A38B.07 TCGA.2Y.A9GT.01A.11R.A38B.07 TCGA.2Y.A9GU.01A.11R.A38B.07
#?|100133144 1.5051 26.4120 0.0000 5.7222
#?|100134869 3.7074 2.6663 4.4833 5.1216
#?|10357 90.1124 71.0054 95.5122 61.6679
#?|10431 1017.1038 639.2311 742.4344 1186.9807
#
dim(rna)
#[1] 16910 423
#Data rna: TCGA rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.Level3 data
#Each row represents a the RSEM normalized counts value of a gene, each column represents a sample
head(clinical_info)
# new_tumor_days death_days followUp_days new_time new_death death_event
#TCGA-2Y-A9GU NA NA 1939 1939 1939 0
#TCGA-2Y-A9GX 2133 NA 2442 2133 2442 0
#TCGA-2Y-A9H0 NA NA 3675 3675 3675 0
#TCGA-2Y-A9H2 NA NA 1731 1731 1731 0
#TCGA-2Y-A9H3 NA NA 1516 1516 1516 0
#TCGA-2Y-A9H4 NA NA 1452 1452 1452 0
#new_tumor_days = days_to_new_tumor_event_after_initial_treatment
#death_days = days_to_death
#followUp_days = days_to_last_followup
#new_time = time to new tumor containing data to censor for new_tumor
#new_death = time to death containing values to censor for death
#death_event = death censoring
dim(clinical_info)
#[1] 313 6TCGA barcode information: https://docs.gdc.cancer.gov/Encyclopedia/pages/images/TCGA-TCGAbarcode-080518-1750-4378.pdf
3b) Data preprocessing
#get the index of the normal/control samples
n_index = which(substr(colnames(rna),14,14) == '1')
t_index = which(substr(colnames(rna),14,14) == '0')
#apply voom function from limma package to normalize the data
library("limma")
vm <- function(x){
cond <- factor(ifelse(seq(1,dim(x)[2],1) %in% t_index, 1, 0))
d <- model.matrix(~1+cond)
x <- t(apply(x,1,as.numeric))
ex <- voom(x,d,plot=F)
return(ex$E)
}
rna_vm = vm(rna)
colnames(rna_vm) = gsub('\\.','-',substr(colnames(rna),1,12))
#check how data look
rna_vm[1:4,1:4]
# TCGA-2V-A95S TCGA-2Y-A9GS TCGA-2Y-A9GT TCGA-2Y-A9GU
#?|100133144 -3.787793 0.2055691 -6.143200 -2.380251
#?|100134869 -2.718538 -2.8818107 -2.826099 -2.526695
#?|10357 1.710169 1.6153717 1.441946 0.940422
#?|10431 5.199493 4.7767132 4.393891 5.196012
#Counts are transformed to log2 counts per million reads (CPM)
#z = [(value gene X in tumor Y)-(mean gene X in normal)]/(standard deviation X in normal)
#We transform data to z-score so that per each patient for each gene we will have a measure of
#how many SD away from the mean that is and we will consider those with Z > +/- 1.96 (roughly p=0.05 or 2 SD away)
#to be differentially expressed.
#calculate z-scores
scal <- function(x,y){
mean_n <- rowMeans(y) # mean of normal
sd_n <- apply(y,1,sd) # SD of normal
# z score as (value - mean normal)/SD normal
res <- matrix(nrow=nrow(x), ncol=ncol(x))
colnames(res) <- colnames(x)
rownames(res) <- rownames(x)
for(i in 1:dim(x)[1]){
for(j in 1:dim(x)[2]){
res[i,j] <- (x[i,j]-mean_n[i])/sd_n[i]
}
}
return(res)
}
z_rna = scal(rna_vm[,t_index],rna_vm[,n_index])
#set the rownames keeping only gene name
rownames(z_rna) = sapply(rownames(z_rna), function(x) unlist(strsplit(x,'\\|'))[[1]])
#check how data look
z_rna[1:4,1:4]
# TCGA-2V-A95S TCGA-2Y-A9GS TCGA-2Y-A9GT TCGA-2Y-A9GU
#? 0.4781493 3.2921169 -1.1816147 1.4699895
#? 1.4040622 1.2730221 1.3177359 1.5580328
#? 0.9726639 0.7364384 0.3042798 -0.9454642
#? 1.2534915 -0.0275773 -1.1875702 1.2429450
#Since we need the same number of patients in both clinical and RNASeq data, we take the indices for the matching samples
ind_tum = which(unique(colnames(z_rna)) %in% rownames(clinical_info))
expression = z_rna[,ind_tum]
#remove gene name=="?"
exp = expression[which(rownames(expression) != '?'),]
#remove death_event==NA, new_death==NA
clinical_info = clinical_info[!is.na(clinical_info$new_death),]
clinical_info = clinical_info[!is.na(clinical_info$death_event),]
#uniq same samples
ind_exp = which(unique(colnames(exp)) %in% rownames(clinical_info))
ind_clin = which(rownames(clinical_info) %in% colnames(exp))
exp = exp[,ind_exp]
clinical_info = clinical_info[ind_clin,]
#check how data look
exp[1:4,1:4]
# TCGA-2Y-A9GU TCGA-2Y-A9GX TCGA-2Y-A9H0 TCGA-2Y-A9H2
#A1BG -9.683236 -3.4303191 -5.5250004 -8.7919389
#A1CF -1.436172 -1.1263191 0.3910154 -7.4946049
#A2LD1 5.466810 -0.8570834 -0.5426471 -0.7720493
#A2M 1.665745 -0.8879948 -5.1409194 -5.5017694
dim(exp)
#[1] 16897 313
#16897 genes, 313 samples
#You could read this exp matrix from:
#exp = readRDS(file="/home/test/clinical_analysis//exp.rds")4) Running steps
4a) Install packages
# install.packages(c("survival", "survminer", "limma"))4b) Library package
library("survival")
library("survminer")
library("limma")4c) Create event vector for RNASeq data
rna_event = t(apply(exp, 1, function(x) ifelse(abs(x) > 1.96, 1, 0)))
#rna_event = 1 : differentially expressed; = 0 : not differentially expressed
#Or another method to classify the gene expression:
rna_event2 = t(apply(exp, 1, function(x) ifelse(x > 1.96, 2, ifelse(x < -1.96, 1, 0 ))))
#rna_event = 2: up-regulated differentially expressed 1 : down-regulated differentially expressed; = 0 : not differentially expressed
#Or you could use Foldchange of a gene's expression as the event vector.
#Selected interested gene(s)
ind_gene = which(rownames(exp) == 'CCDC58')
#check how data look
table(rna_event2[ind_gene,])
# 0 2
#168 145
#In the total 313 samples, CCDC58 gene are not differentially expressed in 168 samples, up-regulated differentially expressed in 145 samples.4d) Fit survival curves
survplotdata = cbind(as.numeric(as.character(clinical_info$new_death)), clinical_info$death_event, rna_event[ind_gene,])
colnames(survplotdata) = c('new_death', 'death_event', 'CCDC58')
survplotdata = as.data.frame(survplotdata)
#check how data look
head(survplotdata)
# new_death death_event CCDC58
#TCGA-2Y-A9GU 1939 0 0
#TCGA-2Y-A9GX 2442 0 0
#TCGA-2Y-A9H0 3675 0 0
#TCGA-2Y-A9H2 1731 0 0
#TCGA-2Y-A9H3 1516 0 1
#TCGA-2Y-A9H4 1452 0 1fit = survfit(Surv(new_death, death_event) ~ CCDC58, data = survplotdata)
#if your interested genes contain more than one genes, you could do the following steps:
#ind_gene1 = which(rownames(exp) == 'CCDC58')
#ind_gene2 = which(rownames(exp) == 'TP53')
#survplotdata2 = cbind(as.numeric(as.character(clinical_info$new_death)), clinical_info$death_event, rna_event[ind_gene1,], rna_event[ind_gene2,])
#colnames(survplotdata2) = c('new_death', 'death_event', 'CCDC58', 'TP53')
#survplotdata2 = as.data.frame(survplotdata2)
#fit2 = survfit(Surv(new_death, death_event) ~ CCDC58+TP53, data = survplotdata2)4e) Draw survival curves
#draw survival curves
pdf(file="/home/test/clinical_analysis/CCDC58_LIHC.pdf")
ggsurvplot(fit,
risk.table = TRUE,
pval = TRUE,
break.time.by = 500,
ggtheme = theme_bw(),
legend.title = "CCDC58 expression",
legend.labs = c("Up-regulated ", "NotAltered"),
risk.table.y.text.col = TRUE,
risk.table.y.text = FALSE,
risk.table.col = "strata",
palette = c("blue", "red"))
dev.off()*: each '+' represent a censored sample.
5) Appendix
5a) Download TCGA RNAseq data and clinical data
We could download data of TCGA liver cancer (LIHC) following:
go to FireBrowse (http://gdac.broadinstitute.org/), select "LIHC" -> "Browse"
from "mRNASeq" select "illuminahiseq_rnaseqv2-RSEM_genes_normalized" and save it
from "Clinical" select "Merge_Clinical" and download it
unzip the files
rename the folders as "RNA" and "Clinical"
5b) Data preprocessing
#read RNA file
rna = read.table('RNA/LIHC.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data.txt',nrows=20533, header=T,row.names=1,sep='\t')
#and take off first row cause we don't need it
rna = rna[-1,]
#first we remove genes whose expression is == 0 in more than 50% of the samples:
rem <- function(x){
x <- as.matrix(x)
x <- t(apply(x,1,as.numeric))
r <- as.numeric(apply(x,1,function(i) sum(i == 0)))
remove <- which(r > dim(x)[2]*0.5)
return(remove)
}
remove = rem(rna)
rna = rna[-remove,]
#The two digits at position 14-15 of the barcode will indicate the sample type:
#"Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29."
#see the values
table(substr(colnames(rna),14,14))
#So we have 373 tumor and 50 normal
#OR you could read this data from:
#rna = readRDS(file="/home/test/clinical_analysis/rna.rds")
####The following data processing steps are in the "2.3 Data preprocessing" part####z_rna = readRDS(file="/home/test/clinical_analysis/z_rna.rds")
#read the Clinical file, in this case i transposed it to keep the clinical feature title as column name
clinical <- t(read.table('Clinical/LIHC.merged_only_clinical_clin_format.txt2',header=T, row.names=1, sep='\t'))
clinical = as.data.frame(clinical)
clinical$IDs = toupper(clinical$patient.bcr_patient_barcode)
#match the patient ID in clinical data with the colnames of z_rna
sum(clinical$IDs %in% colnames(z_rna))
#we have 371 patients that we could use
#get the columns that contain data we can use: days to death, new tumor event, last day contact to....
ind_keep = grep('days_to_new_tumor_event_after_initial_treatment',colnames(clinical))
#this is a bit tedious, since there are numerous follow ups, let's collapse them together and keep the first value (the higher one) if more than one is available
new_tum = as.matrix(clinical[,ind_keep])
new_tum_collapsed <- c()
for (i in 1:dim(new_tum)[1]){
if ( sum ( is.na(new_tum[i,])) < dim(new_tum)[2]){
m <- min(new_tum[i,],na.rm=T)
new_tum_collapsed <- c(new_tum_collapsed,m)
} else {
new_tum_collapsed <- c(new_tum_collapsed,'NA')
}
}
#do the same to death
ind_keep <- grep('days_to_death',colnames(clinical))
death <- as.matrix(clinical[,ind_keep])
death_collapsed <- c()
for (i in 1:dim(death)[1]){
if ( sum ( is.na(death[i,])) < dim(death)[2]){
m <- max(death[i,],na.rm=T)
death_collapsed <- c(death_collapsed,m)
} else {
death_collapsed <- c(death_collapsed,'NA')
}
}
#and days last follow up here we take the most recent which is the max number
ind_keep = grep('days_to_last_followup',colnames(clinical))
fl = as.matrix(clinical[,ind_keep])
fl_collapsed <- c()
for (i in 1:dim(fl)[1]){
if ( sum ( is.na(fl[i,])) < dim(fl)[2]){
m <- max(fl[i,],na.rm=T)
fl_collapsed <- c(fl_collapsed,m)
} else {
fl_collapsed <- c(fl_collapsed,'NA')
}
}
#and put everything together
all_clin = data.frame(new_tum_collapsed, death_collapsed, fl_collapsed)
colnames(all_clin) = c('new_tumor_days', 'death_days', 'followUp_days')
#create vector with time to new tumor containing data to censor for new_tumor
all_clin$new_time <- c()
for (i in 1:length(as.numeric(as.character(all_clin$new_tumor_days)))){
all_clin$new_time[i] <- ifelse ( is.na(as.numeric(as.character(all_clin$new_tumor_days))[i]),
as.numeric(as.character(all_clin$followUp_days))[i],as.numeric(as.character(all_clin$new_tumor_days))[i])
}
#create vector time to death containing values to censor for death
all_clin$new_death <- c()
for (i in 1:length(as.numeric(as.character(all_clin$death_days)))){
all_clin$new_death[i] <- ifelse ( is.na(as.numeric(as.character(all_clin$death_days))[i]),
as.numeric(as.character(all_clin$followUp_days))[i],as.numeric(as.character(all_clin$death_days))[i])
}
#create vector for death censoringtable(clinical$patient.follow_ups.follow_up.vital_status)
#alive dead
#226 94 = 320
all_clin$death_event = ifelse(clinical$patient.follow_ups.follow_up.vital_status == 'alive', 0, 1)
#finally add row.names to clinical
rownames(all_clin) = clinical$IDs
#new_tumor_days = days_to_new_tumor_event_after_initial_treatment
#death_days = days_to_death
#followUp_days = days_to_last_followup
#new_time = time to new tumor containing data to censor for new_tumor
#new_death = time to death containing values to censor for death
#death_event = death censoring
# since we need the same number of patients in both clinical and RNASeq data take the indices for the matching samples
ind_tum = which(unique(colnames(z_rna)) %in% rownames(all_clin))
ind_clin = which(rownames(all_clin) %in% colnames(z_rna))
clinical_info = all_clin[ind_clin,]
#OR you could read this data from:
#clinical_info = readRDS(file="/home/test/clinical_analysis/clinical_info.rds")
####The following data processing steps are in the "2.3 Data preprocessing" part####6) Homework
Please plot the survival curves about the patients with up-regulated differentially expressed and not altered expressed AFP gene in TCGA LIHC data.
input files:
exp = readRDS(file="/home/test/clinical_analysis/exp.rds") #if you don't want to normalize and scale expression data, just load the exp.rds
clinical_info = readRDS(file="/home/test/clinical_analysis/clinical_info.rds")
or
rna = readRDS(file="/home/test/clinical_analysis/rna.rds")
clinical_info = readRDS(file="/home/test/clinical_analysis/clinical_info.rds")7) References
Last updated