花花写于2020-7-15。今天一口气写了几个新的处理tcga数据的函数。加在tinyarray里了，1.3.1版本以上可用~
当前版本有一个小缺点，因为有一步用到了字符型和数值型的转换，会引入NA，虽然不影响最终结果的正确性，但总是会提示“NAs introduced by coercion”，有些烦人~我暂时不知道要怎么去掉它，如果知道的话可以告诉我丫。
可能我写包发推文，给一些初学者造成了误解，今天有个小孩加我好友，问处理数据是不是要自己写包~~~我说明一下，我是因为用了一年多，有些代码复制粘贴了太多次，并且刚好感兴趣才写成包的，常规使用，用别人的R包就好，不是所有人都要自己写的啊！
本文使用的数据都是tinyarray包的内置数据，最新版本的包在https://github.com/xjsun1221/tinyarray获取呀。
1.给TCGA表达矩阵分tumor-normal组
用包的内置数据做例子，写函数搞定他。
library(tinyarray)
group_list = make_tcga_group(exp_hub1);table(group_list)
## Warning in make_tcga_group(exp_hub1): NAs introduced by coercion
## group_list
## normal tumor
## 171 179
这里会提示引入NA，其实不会的，忽略就行。
2.批量计算基因表达量最佳截点
在给基因分组为高低表达量两组时，比较简单的办法就是用中位数作为截断值，后来我又知道了，可以用R包计算最佳截点，我把他写成了函数，直接提供一个表达矩阵和meta信息，返回最佳截点值组成的向量。
仍然是内置数据集。
dim(exprSet_hub1)
## [1] 8 177
point_cut(exprSet_hub1,meta1)
## CXCL8 FN1 COL3A1 ISG15 COL1A2 CXCL10 ICAM1 MMP9
## 310 31190 44889 744 53889 275 2021 326
3.KM批量生存分析
gdcrnatools虽然也可以，但他不支持最佳截点，我把使用最佳截点作为一个参数加上去了。
返回的结果是log rank test计算的p值
surv_KM(exprSet_hub1,meta1)
## CXCL10
## 0.00423786
surv_KM(exprSet_hub1,meta1,cut.point = T)
## CXCL8 COL3A1 CXCL10 COL1A2 ISG15 FN1
## 0.0005655718 0.0011140684 0.0011210275 0.0018228110 0.0031000411 0.0042871122
## ICAM1 MMP9
## 0.0142872301 0.0390594720
4.cox批量生存分析
也是支持最佳截点
surv_cox(exprSet_hub1,meta1,cut.point = T)
## coef se z p HR HRse HRz HRp HRCILL HRCIUL
## CXCL8 1.601 0.514 3.115 0.002 4.957 2.548 1.553 0.120 1.810 13.575
## FN1 1.249 0.464 2.689 0.007 3.486 1.619 1.535 0.125 1.403 8.664
## COL3A1 1.544 0.518 2.981 0.003 4.682 2.425 1.519 0.129 1.697 12.919
## ISG15 1.405 0.514 2.733 0.006 4.074 2.094 1.468 0.142 1.488 11.154
## COL1A2 1.356 0.466 2.909 0.004 3.879 1.807 1.593 0.111 1.556 9.668
## CXCL10 0.680 0.213 3.198 0.001 1.974 0.420 2.320 0.020 1.301 2.995
## ICAM1 0.589 0.244 2.419 0.016 1.803 0.439 1.828 0.068 1.118 2.906
## MMP9 0.795 0.395 2.012 0.044 2.214 0.875 1.388 0.165 1.021 4.804
5.批量画箱线图
提供感兴趣基因的TCGA表达矩阵即可。我一般是用log后的count值来画
k = exp_boxplot(log2(exp_hub1+1));k[[1]];length(k)
## Warning in make_tcga_group(exp_hub): NAs introduced by coercion
image.png
## [1] 8
patchwork::wrap_plots(k,nrow = 2)
image.png
6.批量画生存分析图
提供感兴趣基因的表达矩阵和对应的临床信息即可出图
tmp = exp_surv(exprSet_hub1,meta1);length(tmp )
## [1] 8
patchwork::wrap_plots(tmp,nrow = 2)
image.png
7.箱线图和生存分析图拼图
生存分析的R包是基于ggplot2扩展来的，需要将ggplot对象提取出来，才能用patchwork拼图，我已经在函数里把这一步搞定了。
还调整了一下比例~
k = box_surv(log2(exp_hub1+1),exprSet_hub1,meta1);patchwork::wrap_plots(k,nrow = 2)
## Warning in make_tcga_group(exp_hub): NAs introduced by coercion
所有的warning都可忽略
网友评论