R tips：使用TCGAbiolinks包下载TCGA数据

library(TCGAbiolinks)
library(tidyverse)

# 获取全部的TCGA项目信息
AllProj <- getGDCprojects() 

# 每个项目的摘要
proj_READ <- TCGAbiolinks:::getProjectSummary("TCGA-READ")
proj_COAD <- TCGAbiolinks:::getProjectSummary("TCGA-COAD")

# query
query_READ <- GDCquery(
  project       = "TCGA-READ",
  data.category = "Transcriptome Profiling", 
  data.type     = "Gene Expression Quantification",
  workflow.type = "HTSeq - Counts", # "STAR - Counts"
)

query_COAD <- GDCquery(
  project       = "TCGA-COAD",
  data.category = "Transcriptome Profiling", 
  data.type     = "Gene Expression Quantification",
  workflow.type = "HTSeq - Counts" # "STAR - Counts"
)

# download
GDCdownload(
  query           = query_READ,
  method          = "api",
  files.per.chunk = 5
)
GDCdownload(
  query           = query_COAD,
  method          = "api",
  files.per.chunk = 5
)

# prepare
dat_READ <- GDCprepare(
  query         = query_READ,
  save          = TRUE,
  save.filename = "data_READ.rda",
  remove.files.prepared = FALSE
)
dat_COAD <- GDCprepare(
  query         = query_COAD,
  save          = TRUE,
  save.filename = "data_COAD.rda",
  remove.files.prepared = FALSE
)

library(DESeq2)

# 合并数据
dat_couts <- 
  list(dat_counts_READ,
     dat_counts_COAD) %>%
  map(~rownames_to_column(.x, "Gene")) %>%
  do.call(left_join, .) %>%
  column_to_rownames("Gene")

dds <- DESeq2::DESeqDataSetFromMatrix(
  dat_couts,
  colData = data.frame(row.names = colnames(dat_couts)),
  design = ~1
)

dds <- DESeq2::estimateSizeFactors(dds)

dat_couts_norm <- counts(dds, normalized = TRUE)

dat_surv_READ <- 
  data_READ@colData %>% 
  as.data.frame() %>%
  dplyr::select(
    barcode,
    sample_type, 
    vital_status, 
    days_to_death,
    days_to_last_follow_up
  ) %>%
  dplyr::filter(sample_type == "Primary Tumor") %>%
  mutate(status = ifelse(vital_status == 'Dead', "1", "0")) %>%
  mutate(
    OS = case_when(
      status == 1 ~ days_to_death,
      status == 0 ~ days_to_last_follow_up,
      TRUE ~ NA_integer_
    )
  ) %>%
  dplyr::filter(!is.na(OS))

dat_surv_COAD <- 
  data_COAD@colData %>% 
  as.data.frame() %>%
  dplyr::select(
    barcode,
    sample_type, 
    vital_status, 
    days_to_death,
    days_to_last_follow_up
  ) %>%
  dplyr::filter(sample_type == "Primary Tumor") %>%
  mutate(status = ifelse(vital_status == 'Dead', "1", "0")) %>%
  mutate(
    OS = case_when(
      status == 1 ~ days_to_death,
      status == 0 ~ days_to_last_follow_up,
      TRUE ~ NA_integer_
    )
  ) %>%
  dplyr::filter(!is.na(OS))
# 合并READ COAD
dat_surv <- rbind(dat_surv_READ, dat_surv_COAD)

# barcode的前15位是病人ID，根据barcode去重
dat_surv <- dat_surv %>% mutate(Patient_ID = str_sub(barcode, 15))
dat_surv_uniq <- dat_surv %>% distinct(Patient_ID, .keep_all = TRUE)

### 随机选100个基因
filter_gene <- sample(rownames(dat_couts_norm), 100)
filter_dat_norm <- 
  dat_couts_norm[filter_gene, ] %>% t %>%
  as.data.frame %>% 
  rownames_to_column("barcode")

filter_dat_norm_full <- left_join(dat_surv_uniq, filter_dat_norm, by = "barcode")

library(survival)
library(survminer)

res_surv_full <- 
  map(filter_gene,
      function(x){
        # 由于是随机选的基因，不少基因几乎没有表达量
        # 这里做一下判断，如果一个基因的均值低于5，就不做生存分析了
        if(mean(filter_dat_norm_full[[x]]) < 5){
          return(NULL)
        }

        # 根据基因表达的中位数标记为High和Low
        dat <- filter_dat_norm_full %>% 
          mutate('{x}' := ifelse(.data[[x]] > median(.data[[x]]), "High", "Low"))

        # 生存信息status保证是数值类型
        dat$status <- as.numeric(dat$status)

        # 手动构造生存分析的拟合公式
        pattern <- str_glue("Surv(OS, status) ~ {x}") 
        formu <- pattern %>% as.formula()

        # 拟合
        fit <- survfit(formu, data=dat)
        fit$call$formula <- formu

        # 差异分析
        diff <- survdiff(formu, data=dat)
        diff$call$formula <- formu

        p.value <- 1 - pchisq(diff$chisq, length(diff$n) -1)

        # 返回结果
        list(fit = fit,
             diff = diff,
             p.value = p.value,
             dat = dat)
      })
# 未做生存分析的基因直接筛除
res_surv_full <- res_surv_full %>% .[!map_lgl(., is.null)]

length(res_surv_full)
# [1] 42

# 此时的res_surv_full的每个元素都是一个四元素列表
res_surv_full[[1]] %>% str(max.level = 1)
#List of 4
# $ fit    :List of 17
#   ..- attr(*, "class")= chr "survfit"
# $ diff   :List of 6
#  ..- attr(*, "class")= chr "survdiff"
# $ p.value: num 0.294
# $ dat    :'data.frame':    94 obs. of  108 variables:

res_surv_plots <- 
  map(res_surv_full, 
      function(x){
        ggsurvplot(x$fit,
                   data = x$dat,
                   pval = TRUE, 
                   risk.table = TRUE)
      })

# 查看第一个结果
res_surv_plots[[1]]