m6a甲基化相关基因根据临床信息分组绘制boxplot并显示p值

作者: 生信交流平台 | 来源:发表于2021-09-23 11:15 被阅读0次

前面给大家介绍过

在上面这篇文章里面我们直接根据样本类型来分成了两类，即tumor和normal。

但更多的时候，我们其实想研究这些基因的表达跟临床特征的关系，如组织病理分期stage I, stage II, stage III, stage IV，TNM分期，T0, T1, T2, T3, T4，N0, N1, N2, N3(☞肿瘤TNM分期)等等。这个就跟前面讲的不太一样了，这些临床特征都大于两类，那么今天我们就来看看如何根据临床特征对样本进行分组绘制boxplot，并显示p值。

我们首先来捋一捋思路

1.从TCGA数据库下载想要研究肿瘤的表达谱数据和临床数据，这里以CHOL(胆管癌为例)

具体可以参考

☞如何从TCGA数据库下载RNAseq数据以及临床信息(一)

☞如何从TCGA数据库下载miRNA数据(二)

2.合并得到表达谱矩阵，并进行归一化

具体可以参考

☞如何合并TCGA表达谱数据

☞零代码合并TCGA表达谱数据

3.搜集m6A甲基化相关的基因，根据需要做基因的ID转换

这里我已经帮大家从一些相关的研究中搜集了m6A甲基化相关的16个基因，当然大家也可以根据自己的需要增，删相关基因。我们讲的主要是方法，代码都是一样的跑，多几个基因少几个基因都无所谓。甚至你还可以换成其他的任何基因列表都是可以的。

"METTL3","METTL14","METTL16","WTAP","VIRMA","RBM15","ZC3H13","YTHDF1","YTHDF2","YTHDF3","YTHDC1","YTHDC2","HNRNPA2B1","HNRNPC","FTO","ALKBH5"

这里需要注意的是TCGA中gene的ID一般是ensembl gene ID，而我们从paper中搜集的gene一般是gene名字，所以这里就需要做一个gene ID转换。这个我在前面也讲过

☞【R语言】基因ID转换

4.提取m6A基因相关的表达矩阵

这个其实就是R里面最常见的数据框取下标操作了。这里所有分析都是基于R来做的，所以了解一些R的基础知识还是很有必要的。
具体可以参考

☞R语言入门-工欲善其事必先利其器

m6A基因相关的表达矩阵如下

5.添加样本的类型信息，如肿瘤，癌旁正常对照

前面也用视频给大家详细介绍过如何从TCGA数据下载临床信息

☞如何从TCGA数据库下载RNAseq数据以及临床信息(一)

6.下载并处理临床信息

前面我们也给大家讲解过

☞玩转TCGA临床信息

也分享过如果处理多分组的临床信息

☞【R语言】临床特征分组，多分类转换成二分类

这里我们将组织病理分期分成三组之后跟m6a基因的表达值合并起来

7.对数据格式进行转换

8.根据分组绘制boxplot

9.计算p值

在m6a甲基化相关基因boxplot并显示p值一文中，我们是根据样本类型来分类的，只有两类，所以我们可以做t检验来获取p值。本文中我们将样本根据组织病理分期分成了三类（当然四类也可以），stage I，stage II和stage III/IV。所以这里就不能再使用t检验来获取p值了，而应该使用aov方差检验来获取p值。

10.将计算得到的p值转换成对应的*，添加到boxplot上

可以参考☞【R语言】P值转换成***将p值转换成相应的*