太好用了！批量生存分析，一步到位，还支持最佳截点~

花花写于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都可忽略