TCGA数据的另一个重要的包：RTCGA
参考学习资料：
TCGA之miRNA预测临床表现
RTCGA（1） 数据概况与数据下载
RTCGA（2） 数据分析与可视化
https://www.bioconductor.org/packages/release/bioc/vignettes/RTCGA/inst/doc/RTCGA_Workflow.html
https://www.bioconductor.org/packages/release/bioc/manuals/RTCGA/man/RTCGA.pdf
TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析
RTCGA包安装
if(!require(BiocManager)){
install.packages("BiocManager")
}
BiocManager::install("RTCGA")
library('RTCGA')
library('magrittr')
infoTCGA函数用于查看所有癌症类型：
> infoTCGA()%>%
+ rownames()%>%
+ sub('-counts','',.) -> cohort
> cohort
[1] "ACC" "BLCA" "BRCA" "CESC" "CHOL" "COAD" "COADREAD"
[8] "DLBC" "ESCA" "FPPP" "GBM" "GBMLGG" "HNSC" "KICH"
[15] "KIPAN" "KIRC" "KIRP" "LAML" "LGG" "LIHC" "LUAD"
[22] "LUSC" "MESO" "OV" "PAAD" "PCPG" "PRAD" "READ"
[29] "SARC" "SKCM" "STAD" "STES" "TGCT" "THCA" "THYM"
[36] "UCEC" "UCS" "UVM"
共收录了38种癌症
checkTCGA用于查询数据集的名字及数据集的公开日期：
checkTCGA查询的信息可以用于downTCGA以下载特定的TCGA数据。
checkTCGA('DataSets', 'BRCA' ) #BRCA的TCGA数据集
checkTCGA('Dates') #所有的TCGA数据的公布日期，最新的为2016-01-28
数据下载与导入
有两种方式下载：downTCGA下载数据、R包的形式下载数据，推荐以R包的形式下载数据。
第一种
下载ACC的miRNA基因表达数据
dir.create( 'hre')
downloadTCGA( cancerTypes = 'ACC', dataSet = 'miR_gene_expression',
destDir = 'hre', date = tail( checkTCGA('Dates'), 2 )[1] )
下载的数据需要使用readTCGA的方式去读取：
dir.create('data')
# 下载临床数据
# dataset = "clinical" 是默认的参数，可以忽略
downloadTCGA( cancerTypes = c('BRCA', 'OV'),
destDir = 'data' )
# 读取数据集
sapply( c('BRCA', 'OV'), function( element ){
folder <- grep( paste0( '(_',element,'\\.', '|','_',element,'-FFPE)', '.*Clinical'),
list.files('data/'),value = TRUE)
path <- paste0( 'data/', folder, '/', element, '.clin.merged.txt')
assign( value = readTCGA( path, 'clinical' ),
x = paste0(element, '.clin.data'), envir = .GlobalEnv)
})
第一种方法比较繁琐，不是大众首推的方式。
第二种方式：通过R包的形式下载数据
RTCGA将TCGA分类汇总后制成了bioconductor包来管理，现在在用的有两个版本2015-11-01与2016-01-28，建议使用最新版，数据是最新的。
- RTCGA.mutations
- RTCGA.rnaseq
- RTCGA.clinical
- RTCGA.PANCAN12
- RTCGA.mRNA
- RTCGA.miRNASeq
- RTCGA.RPPA
- RTCGA.CNV
-
RTCGA.methylation
在这里面，数据包RTCGA.rnaseq,RTCGA.mRNA,RTCGA.RPPA,RTCGA.miRNASeq,RTCGA.methylation的内容通过expressionsTCGA函数获取，数据包RTCGA.clinical的内容通过survivalTCGA函数获取，数据包RTCGA.mutations的内容通过mutationsTCGA函数获取。
以RTCGA.rnaseq包为例，安装之后可以通过expressionsTCGA函数获取相应的TCGA数据，其包含的数据是以“癌症缩写+rnaseq”格式的数据框，如下图：
library("BiocManager")
BiocManager::install("RTCGA.rnaseq")
library("RTCGA.rnaseq")
#获取BRCA、OV及HNSC中的VENTX基因的表达量
expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq,
extract.cols = "VENTX|27287")
如下图，获取的结果中，第一列为病人barcode，不同的病人的不同样本有不同的barcode用于区分，第二列为dataset，代表来源于哪个数据框，本例中为BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq，第三列及以后列（如有）为自己选择的感兴趣的基因表达量。
> expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq,
+ extract.cols = "VENTX|27287")
# A tibble: 2,085 x 3
bcr_patient_barcode dataset `VENTX|27287`
<chr> <chr> <dbl>
1 TCGA-3C-AAAU-01A-11R-A41B-07 BRCA.rnaseq 6.55
2 TCGA-3C-AALI-01A-11R-A41B-07 BRCA.rnaseq 8.70
3 TCGA-3C-AALJ-01A-31R-A41B-07 BRCA.rnaseq 13.6
4 TCGA-3C-AALK-01A-11R-A41B-07 BRCA.rnaseq 6.21
5 TCGA-4H-AAAK-01A-12R-A41B-07 BRCA.rnaseq 34.5
6 TCGA-5L-AAT0-01A-12R-A41B-07 BRCA.rnaseq 26.4
7 TCGA-5L-AAT1-01A-12R-A41B-07 BRCA.rnaseq 39.9
8 TCGA-5T-A9QA-01A-11R-A41B-07 BRCA.rnaseq 6.28
9 TCGA-A1-A0SB-01A-11R-A144-07 BRCA.rnaseq 14.4
10 TCGA-A1-A0SD-01A-11R-A115-07 BRCA.rnaseq 12.9
# ... with 2,075 more rows
另外，这些bioconductor数据包还可以通过convertTCGA转换成相应的bioconductor对象，如RTCGA.rnaseq,
RTCGA.RPPA, RTCGA.PANCAN12, mRNA, RTCGA.methylation可以转换成ExpressionSet对象，RTCGA.CNV可以转换成GRanges对象。
生存分析
RTCGA.clinical包中分癌症类型，有38个数据表，如BRCA.clinical, OV.clinical等等，癌症的编号详见http://gdac.broadinstitute.org/。每个数据表包括"admin.bcr"、"admin.day_of_dcc_upload"、"admin.disease_code" 等共计3703列。如下所示：
#BiocManager::install("RTCGA.clinical")
library("RTCGA.clinical")
library("tidyverse")
...
> BRCA.clinical%>%colnames()%>%head
[1] "admin.bcr" "admin.day_of_dcc_upload"
[3] "admin.disease_code" "admin.file_uuid"
[5] "admin.month_of_dcc_upload" "admin.patient_withdrawal.withdrawn"
> BRCA.clinical%>%ncol()
[1] 3703
survivalTCGA()可以从RTCGA.clinical中获取临床数据，如果不注明extract.cols参数，那么结果只有三列："times 、bcr_patient_barcode、patient.vital_status"，分别是生存时间、barcode及病人生存状态。
kmTCGA()用于绘制生存曲线，参数explanatory.names代表分组：
# 获取生存数据，不同的疾病
survivalTCGA(BRCA.clinical, OV.clinical, extract.cols = "admin.disease_code") -> BRCAOV.survInfo
# 绘图
kmTCGA(BRCAOV.survInfo, explanatory.names = "admin.disease_code", pval = TRUE)
kmTCGA(BRCAOV.survInfo, explanatory.names = "admin.disease_code", main = "",xlim = c(0,4000))
生存分析
截断x轴
# 获取生存数据，不同的治疗方式
survivalTCGA(BRCA.clinical,
extract.cols = c("patient.drugs.drug.therapy_types.therapy_type")) %>%
filter(patient.drugs.drug.therapy_types.therapy_type %in%
c("chemotherapy", "hormone therapy")) %>%
rename(therapy = patient.drugs.drug.therapy_types.therapy_type) -> BRCA.survInfo.chemo
# 绘图
kmTCGA(BRCA.survInfo.chemo, explanatory.names = "therapy", xlim = c(0, 3000), conf.int = FALSE)
image.png
表达分析
以RTCGA.rnaseq包为例，expressionsTCGA ()用于获取表达数据。
expressionsTCGA()如指定参数extract.cols，则返回特定基因在各个样本的表达量，如不指定则返回全部基因，其值的形式为“Gene symbol|Gene ID”，如 "VENTX|27287"。
PCA
library(RTCGA.rnaseq)
# 获取全部基因的表达情况
expressionsTCGA(BRCA.rnaseq, OV.rnaseq, HNSC.rnaseq) %>%
rename(cohort = dataset) %>%
filter(substr(bcr_patient_barcode, 14, 15) == "01") -> BRCA.OV.HNSC.rnaseq.cancer
BRCA.OV.HNSC.rnaseq.cancer %>%ncol()
# 20533
# 由于基因数太多了，随机选1000个基因绘图
BRCA.OV.HNSC.rnaseq.cancer%>%select(sample(colnames(.),1000),-bcr_patient_barcode,cohort) %>%
pcaTCGA("cohort")->pca.rnaseq
plot(pca.rnaseq)
分组
热图
hearmap用于绘制热图，主要有四个参数，第一个为包含所需变量的数据库，第二个和第三个分别是热图的X、Y轴分组，第四个变量为展示的变量。
# 获取MET、ZNF500、ZNF501三个基因在ACC、BLCA、BRCA、OV癌症中的表达
expressionsTCGA(ACC.rnaseq, BLCA.rnaseq, BRCA.rnaseq, OV.rnaseq,
extract.cols = c("MET|4233", "ZNF500|26048", "ZNF501|115560")) %>%
rename(cohort = dataset, MET = `MET|4233`) %>% #cancer samples
filter(substr(bcr_patient_barcode, 14, 15) == "01") %>%
mutate(MET = cut(MET,
round(quantile(MET, probs = seq(0,1,0.25)), -2),
include.lowest = TRUE,
dig.lab = 5)) -> ACC_BLCA_BRCA_OV.rnaseq
# 以癌症为第一分组、MET基因表达量为第二分组，计算各分组ZNF500、ZNF501的基因表达中位数
ACC_BLCA_BRCA_OV.rnaseq %>%
select(-bcr_patient_barcode) %>%
group_by(cohort, MET) %>%
summarise_all(median) %>%
mutate(ZNF500 = round(`ZNF500|26048`),
ZNF501 = round(`ZNF501|115560`)) -> ACC_BLCA_BRCA_OV.rnaseq.medians
# 绘制ZNF501热图
heatmapTCGA(ACC_BLCA_BRCA_OV.rnaseq.medians,
"cohort", "MET", "ZNF501", title = "Heatmap of ZNF501 expression")
热图
箱线图
boxplotTCGA用于绘制箱线图，主要的参数有data、x、y、fill，分别是数据库、x轴、y轴，染色填充，两个有意思的参数coord.flip坐标轴翻转（默认翻转），facet.names分组（facet）变量，xlab和ylab用于定义坐标轴标题，legend.title定义图例，legend定义图例的位置。
# 获取RNAseq表达量数据
expressionsTCGA(ACC.rnaseq, BLCA.rnaseq, BRCA.rnaseq, OV.rnaseq,
extract.cols = "MET|4233") %>%
rename(cohort = dataset,
MET = `MET|4233`) %>%
#cancer samples
filter(substr(bcr_patient_barcode, 14, 15) == "01") -> ACC_BLCA_BRCA_OV.rnaseq
h1 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "cohort", "MET")
h2 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "cohort", "log1p(MET)") # 数据变换，可以压缩异常值的离群趋势
h3 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq, "reorder(cohort,log1p(MET), median)", "log1p(MET)") # 调整x的顺序
h4 = boxplotTCGA(ACC_BLCA_BRCA_OV.rnaseq,"cohort", "log1p(MET)",facet.names = "cohort")
library(patchwork)
h1 + h2 + h3 + h4 +plot_layout()
箱线图
RNAseq的数据演示结束
获取miRNA数据
有些miRNA并无表达，因此选择miRNA满足以下条件：其在至少十个病人中表达量均大于1。
rm(list = ls())
options(stringsAsFactors = F)
#BiocManager::install("RTCGA.miRNASeq")
library(RTCGA.miRNASeq)
rowName <- rownames(LUAD.miRNASeq)[seq(1,nrow(LUAD.miRNASeq),by=3)]
expr<-LUAD.miRNASeq[seq(1,nrow(LUAD.miRNASeq),by=3),3:ncol(LUAD.miRNASeq)]
expr<- apply(expr,2,as.numeric)
expr<- na.omit(expr)
dim(expr)
expr <- expr[,apply(expr, 2,function(x){sum(x>1)>10})]
rownames(expr) <- rowName
dim(expr)
获取临床信息
library(RTCGA.clinical)
meta <- LUAD.clinical
meta <- meta[,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')]
数据清洗
- 调整NA值、调整数据格式（数值变量、因子变量）。
- 生存时间拆分在两列中：死亡时间和随访时间，因此需要合并到一起。
- 上面获得的miRNA数据expr共有561行，而clinical数据meta共有522行，因此meta需要剔除部分不匹配的数据（下面的match函数部分）。
group_list <- ifelse(substr(rownames(expr),14,15)=='01','tumor','normal')
#table(group_list)
expr <- expr[group_list=='tumor',]
meta[,3][is.na(meta[,3])]=0
meta[,4][is.na(meta[,4])]=0
meta$days <- as.numeric(meta[,3])+as.numeric(meta[,4])
meta <- meta[,c(1:2,5:9)]
#colnames(meta)
colnames(meta)=c('ID','event','race','age','gender','stage',"days")
meta$event=ifelse(meta$event=='alive',0,1)
meta$stage <- str_split(meta$stage,' ',simplify = T)[,2]
table(meta$stage)
meta$age <- as.numeric(meta$age)
meta$age_group <- ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')
meta$time <- meta$days/30
meta$race[is.na(meta$race)] <- "unknown"
meta$ID <- toupper(meta$ID)
meta <- meta[match(substr(rownames(expr),1,12), meta$ID),]
meta[1:4,1:4]
挑选miRNA，构建cox回归模型
- 由于生存分析绘图常用
survminer，因此使用survminer绘制cox图时有一个坑：绘图前需要手动指定好cox回归的模型。 - cox回归还可以展示森林图。具体如下，虽然总p值为0.018，但是8个miRNA均没有达到显著水平（置信区间横跨无效线）。
library(survival)
library(survminer)
e <- expr[,c('hsa-mir-31','hsa-mir-196b','hsa-mir-766','hsa-mir-519a-1',
'hsa-mir-375','hsa-mir-187','hsa-mir-331','hsa-mir-101-1')]
e <- log2(e+1)
colnames(e) <- c('miR31','miR196b','miR766','miR519a1','miR375','miR187','miR331','miR101')
dat <- cbind(meta,e)
dat$gender <- factor(dat$gender)
dat$stage <- factor(dat$stage)
colnames(dat)
s <- Surv(time, event) ~ miR31+miR196b+miR766+miR519a1+miR375+miR187+miR331+miR101
model <- coxph(s, data = dat )
summary(model,data=dat)
options(scipen=1)
fit <- survfit(model, data=dat)
fit$call$formula <- s #手动指定模型
ggsurvplot(fit, palette = "#2E9FDF", ggtheme = theme_minimal())
ggforest(model, data =dat,
main = "Hazard ratio", #标题
cpositions = c(0.10, 0.22, 0.4), #前三列距离
fontsize = 1.0, #字体大小
refLabel = "1",#相对变量的数值标签
noDigits = 4 #HR值及置信区间的有效数字位数
)
生存曲线
森林图，目前还不是太明白，需要再学习
绘制风险因子关联图
跟森林图一样都是不熟悉的东西，需要进一步学习
new_dat=dat
fp <- predict(model,new_dat,type="risk")
fp <- predict(model,new_dat,type="expected")
plot(fp,meta$days)
fp <- predict(model,new_dat)
basehaz(model)
library(Hmisc)
options(scipen=200)
with(new_dat,rcorr.cens(fp,Surv(time, event)))
library(cowplot)
library(pheatmap)
# https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
fp_dat=data.frame(s=1:length(fp),v=as.numeric(sort(fp )))
sur_dat=data.frame(s=1:length(fp),
t=meta[names(sort(fp )),'time'] ,
e=meta[names(sort(fp )),'event'] )
sur_dat$e=ifelse(sur_dat$e==0,'alive','death')
exp_dat=new_dat[names(sort(fp )),10:17]
plot.point=ggplot(fp_dat,aes(x=s,y=v))+geom_point()
plot.sur=ggplot(sur_dat,aes(x=s,y=t))+geom_point(aes(col=e))
mycolors <- colorRampPalette(c("blue", "white", "red"), bias = 1.2)(100)
tmp=t(scale(exp_dat))
tmp[tmp > 1] = 1
tmp[tmp < -1] = -1
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = F)
plot_grid(plot.point, plot.sur, plot.h$gtable,
labels = c("A", "B","C"),
align = 'v',ncol = 1)
虽然不知道怎么解释，但是觉得应该能看懂
