TCGA代码分析流程 - 3.1 生存分析前的数据整理

生存分析只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；
clinical信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。
由于不同癌症的临床信息表格组织形式不同，这里的代码需要根据实际情况修改。

rm(list = ls())
options(stringsAsFactors = F)
library(stringr)
setwd("D:/BioInformationAnalyze/TCGAdata/CRC")  # 设置工作目录（需修改）
load("TCGA-COAD_01.2.lncRNA_Expdata_tpm.Rdata") # 表达矩阵 tpm形式
load("TCGA-COAD_02.lncRNA_Clinical.Rdata") # 临床信息
load("TCGA-COAD_03.2.lncRNA_cg.Rdata") # 差异基因

临床数据整理

tmp <- data.frame(colnames(clinical))
clinical <-  clinical[, c(
  "bcr_patient_barcode", # 病人id
  "vital_status", # 生存状态：存活/死亡
  "days_to_death", # 从初次病理诊断日期到死亡日期的天数
  "days_to_last_followup", # 从初次病理诊断日期到最后一次随访日期的天数
  "race_list", # 人种
  "age_at_initial_pathologic_diagnosis", # 初次病理诊断年龄，没有这一列的话可以用days_to_birth/365计算
  "gender" , # 性别
  "stage_event" # 分级，6th和7th代表第6/7代分级标准
)]
# rownames(clinical) <- NULL
rownames(clinical) <- clinical$bcr_patient_barcode
meta <- clinical[clinical$days_to_death != 0,] # 去除诊断时已经死亡的病人
# 简化meta的列名
colnames(meta) <- c("ID", "event", "death", "last_followup", "race", "age", "gender", "stage")

表达矩阵整理

exprSet <- exp[cg, group_list == "tumor"] # 从表达矩阵中提取癌症样本的差异基因表达量。
# 因为cg是由counts数据得来的，而生存分析用的exp是tpm形式的，所以取交集就会有NA，下一步先去掉全是NA的行。
exprSet <- exprSet[apply(exprSet, 1, function(y) any(!is.na(y))), ] # 去除全是NA的行
table(str_sub(colnames(exprSet),1,12) %in% meta$ID) # 表达矩阵和临床信息的匹配情况
exprSet <- exprSet[, str_sub(colnames(exprSet),1,12) %in% meta$ID] # 从表达矩阵中挑选有临床信息的样本
exprSet <- exprSet[, str_sort(colnames(exprSet))] # 表达矩阵按列名排序
exprSet <- exprSet[, !duplicated(str_sub(colnames(exprSet),1,12) )] # 表达矩阵去重，每个病人只保留一个肿瘤样本

#调整meta的ID顺序与exprSet列名一致
meta <- meta[match(str_sub(colnames(exprSet),1,12), meta$ID), ] # 调整排序并去除没有肿瘤样本的病人
all(meta$ID == str_sub(colnames(exprSet),1,12)) # 确保顺序完全一致

整理生存分析的输入数据

# 1.由随访时间和死亡时间计算生存时间(天)
is.empty.chr <- function(x){
  ifelse(stringr::str_length(x) == 0, T, F)
} # 建立函数来判断字符串是否为空
is.empty.chr(meta[2,3]) # 判断函数可行性
meta[,3][is.empty.chr(meta[,3])] = 0 # 将“death”列的空字符串改为0
meta[,4][is.empty.chr(meta[,4])] = 0 # 将“last_followup”列的空字符串改为0
meta$time <- (as.numeric(meta[,3]) + as.numeric(meta[,4])) # 由随访时间和死亡时间计算生存时间(天)

# 2.根据生死定义event，活着是0，死的是1
meta$event <- ifelse(meta$event == "Alive", 0, 1)

# 3.年龄和年龄分组
meta$age <- as.integer(meta$age)
meta$age_group <- ifelse(meta$age > median(meta$age), "older", "younger")

# 4.stage 
# 关键在于定位T的位置
tmp <- str_split(meta$stage, " ", simplify = T)[,2] # 按空格分开再取后者，但有个别是空字符串
meta$stage[str_count(tmp, "T") == 0] # 查看空字符串的原因，是因为没有空格
tmp <- str_replace(meta$stage, "thT", "th T") # 将th和T分开
tmp <- str_split(tmp, " ", simplify = T)[,2] # 再分一次

table(str_count(tmp,"T")) # 看看还有没有没有T的
table(str_locate(tmp,"T")[,1]) # 看看T在哪个位置
tmp <- str_sub(tmp, 1, str_locate(tmp,"T")[,1]-1)
table(tmp)
table(is.na(tmp))
meta$stage <- tmp

# 5.race
table(meta$race)

save(meta, exprSet, cancer_type, file = paste0(cancer_type, "_04.lncRNA_SurvData.Rdata"))  # 得到做生存分析的数据，包括K-M曲线和Cox回归。

想要TNM分期的话，代码同stage

 

可以先做Lasso回归对自变量进行降维：

# lasso回归数据
datExpr <- as.data.frame(with(meta, cbind(event, time, t(exprSet))))
rownames(datExpr) <- str_sub(rownames(datExpr),1,12)
save(datExpr, cancer_type, file = paste0(cancer_type, "_04.lncRNA_SurvModel_LassoData.Rdata"))