简介
做生信数据分析的小伙伴一定离不开TCGA和GEO这两个数据库，但是在下载以及整理的数据的过程中就会遇到很多麻烦，比如探针id的转换，各种id的转换，数据的过滤，基本的差异分析，PCA，聚类，批量生存分析等。
真正分析的代码也就那几句，但是整理数据的过程真的是太麻烦了！！
今天介绍的这个tinyarray包就可以帮助我们解决这些问题，大幅度简化你的数据清理过程！
这个包是生信技能树团队小洁老师写的R包哟，真的是太棒啦！！
下载安装
在线安装：
install.packages("tinyarray")if(!require(devtools))install.packages("devtools")
if(!require(tinyarray))devtools::install_github(
下载zip包后本地安装：
devtools::install_local("tinyarray-master.zip",upgrade = F,dependencies = T)
简化TCGA数据分析
表达矩阵的整理
简化id转换，简化分组，一步到位，提取mRNA和lncRNA!
# 以xena下载的BRCA为例，其他也一样可以！
rm(list = ls())
library(tinyarray) # 加载包
##
##
survinfo <- data.table::fread("E:/projects/lisy/files/TCGA-BRCA.survival.tsv",data.table = F)
expr <- data.table::fread("E:/projects/lisy/files/TCGA-BRCA.htseq_fpkm.tsv.gz",data.table = F)
rownames(expr) <- expr$Ensembl_ID
expr <- expr[,-1] # 变成标准的表达矩阵格式，行是基因，列是样本
expr[1:4,1:4]
## TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A1EU-11A
## ENSG00000242268.2 0.09170787 0.000000000 0.05789928
## ENSG00000270112.3 0.01957347 0.004700884 0.01630174
## ENSG00000167578.15 2.23589760 1.863334319 1.70475318
## ENSG00000273842.1 0.00000000 0.000000000 0.00000000
## TCGA-A8-A06X-01A
## ENSG00000242268.2 0.000000
## ENSG00000270112.3 0.000000
## ENSG00000167578.15 1.947481
## ENSG00000273842.1 0.000000
分别提取mRNA和lncRNA，并转换id为gene symbol，一步到位！
expr_mrna <- trans_exp(expr,mrna_only = T) # 提取mrna，转换id为gene symbol
## 19712 of genes successfully mapping to mRNA,14805 of genes successfully mapping to lncRNA
expr_mrna[1:4,1:4]
## TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A1EU-11A TCGA-A8-A06X-01A
## RAB4B 2.2358976 1.86333432 1.70475318 1.9474808
## C12orf5 2.3219445 4.22669935 1.97575523 2.8087572
## RNF44 3.6200560 3.54611744 3.39694310 4.7232701
## DNAH3 0.3370874 0.01601615 0.04145534 0.0023613
expr_lncrna <- trans_exp(expr,lncrna_only = T) # 提取lncrna，转换id为gene symbol
## 19712 of genes successfully mapping to mRNA,14805 of genes successfully mapping to lncRNA
expr_lncrna[1:4,1:4]
## TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A1EU-11A
## RP11-368I23.2 0.09170787 0.000000000 0.05789928
## RP11-742D12.2 0.01957347 0.004700884 0.01630174
## EHD4-AS1 0.08466075 0.000000000 0.46162417
## RP11-166P13.4 0.05394113 0.000000000 0.50697074
## TCGA-A8-A06X-01A
## RP11-368I23.2 0.00000000
## RP11-742D12.2 0.00000000
## EHD4-AS1 0.08891242
## RP11-166P13.4 0.19295943
tcga样本分为癌和癌旁，也是一步到位！
group_list <- make_tcga_group(expr_mrna)
table(group_list)
## group_list
## normal tumor
## 113 1104
下面接差异分析就非常简单了！就不在演示了。
批量生存分析的简化
生存分析需要将表达矩阵和临床信息合并，或者变成顺序一样才行，直接手撸代码也是可以的，但是tinyarray包把这个过程简化了！
# 去除重复的tumor样本
expr_mrna_fi <- sam_filter(expr_mrna)
## filtered 13 samples.
# 匹配tcga的表达矩阵和临床信息
match_exp_cl(expr_mrna_fi, survinfo, "_PATIENT")
## match successfully.
names(cl_matched)[4:5] <- c("event","time") # 把生存状态和生存时间名字改一下
identical(rownames(cl_matched), colnames(exp_matched))
## [1] TRUE
可以看到现在表达矩阵和临床信息都是1181个样本，而且顺序是一致的！
下面就可以批量做生存分析了。
首先展示下最佳截点的计算方法，比较简单的就是根据基因表达量的中位数作为截断值，但是这样有时会导致p值不显著，所以可以使用其他方法。
# 挑选10个基因做，一共有19712个基因
meta <- cl_matched
exprSet <- exp_matched[1:10,]
point_cut(exprSet,meta)
## $RAB4B
## [1] 1.868292
##
## $C12orf5
## [1] 1.468261
##
## $RNF44
## [1] 3.582367
##
## $DNAH3
## [1] 0.007520641
##
## $RPL23A
## [1] 7.546291
##
## $ARL8B
## [1] 5.517558
##
## $CALB2
## [1] 3.726976
##
## $MFSD3
## [1] 2.525482
##
## $PIGV
## [1] 3.45386
##
## $ZNF708
## [1] 2.11762
可以看到每个基因的最佳截点都给你算好了！
批量KM生存分析：
surv_KM(exprSet,meta,cut.point = T) # 使用最佳截点
## MFSD3 CALB2 RAB4B RPL23A DNAH3 ZNF708
## 2.985878e-11 1.549063e-06 8.240647e-05 6.176400e-03 1.078377e-02 1.764367e-02
## C12orf5 RNF44
## 4.037767e-02 4.395910e-02
批量cox生存分析：
surv_cox(exprSet, meta, cut.point = T)
## coef se z p HR HRse
## RAB4B -0.5685845 0.1460777 -3.892343 9.928061e-05 0.5663265 0.08272767
## C12orf5 -0.4358979 0.2147148 -2.030125 4.234383e-02 0.6466837 0.13885258
## RNF44 -0.3019320 0.1504873 -2.006362 4.481766e-02 0.7393884 0.11126856
## DNAH3 0.6244150 0.2484815 2.512924 1.197353e-02 1.8671534 0.46395306
## RPL23A -0.4741161 0.1745146 -2.716770 6.592238e-03 0.6224350 0.10862400
## CALB2 0.7304669 0.1551399 4.708441 2.496186e-06 2.0760497 0.32207808
## MFSD3 -1.0713457 0.1684207 -6.361128 2.002771e-10 0.3425472 0.05769205
## ZNF708 0.4487043 0.1901424 2.359833 1.828314e-02 1.5662815 0.29781645
## HRz HRp HRCILL HRCIUL
## RAB4B -5.242182 1.586884e-07 0.4253293 0.7540644
## C12orf5 -2.544542 1.094210e-02 0.4245476 0.9850483
## RNF44 -2.342186 1.917117e-02 0.5505257 0.9930420
## DNAH3 1.869054 6.161529e-02 1.1472872 3.0387000
## RPL23A -3.475889 5.091623e-04 0.4421268 0.8762764
## CALB2 3.340959 8.348949e-04 1.5317308 2.8137988
## MFSD3 -11.395899 0.000000e+00 0.2462411 0.4765192
## ZNF708 1.901444 5.724382e-02 1.0789972 2.2736273
批量画箱线图：
boxes <- exp_boxplot(exprSet, color = c("blue","red"))
boxes[[1]]
unnamed-chunk-23-134515972拼图：
patchwork::wrap_plots(boxes, nrow = 2)
unnamed-chunk-24-134515972批量画生存分析图：
surv_plots <- exp_surv(exprSet, meta)
patchwork::wrap_plots(surv_plots, nrow = 3)
image-20220204114708988
