前面我们从GDC下载了TCGA肿瘤数据库的数据,也能够把GDC下载的多个TCGA文件批量读入R 今天我们讲一下TCGA数据的标准化,以及差异分析,得到了标准化后的数据,我们就可以按照以前的帖子,做一系列操作 在得到了差异分析的结果后,我们可以完成热图,火山图,GO分析,KEGG分析,GSEA分析,就跟这个帖子中的一样。 下面开始今天的教程: ### 加载数据 现在的数据是这个样子的 去掉ensemble ID的点号
现在的数据是这个样子的 去掉点号,是为了用gtf文件。 加载gtf文件,这是目前我们能接触的最大文件,有260万行。 load(file = 'gtf_df.Rda') 提取mRNA
最终得到19668行,这是编码基因的个数,现在的数据是这个样子的 提取lncRNA这里很有争议,而我的理由是,即使是编码基因,也会出现非编码转录本,而长链非编码RNA,指的是转录本,所以不能用gene的编码与否来界定 ncRNA <>'sense_overlapping','lincRNA','3prime_overlapping_ncRNA', 最终得到25530个非编码转录本,数据是这个样子的 数据标准化标准化和差异分析都是用Deseq2这个包来完成,首先要构建dds对象,构建这个对象需要两个文件,第一是输入数据,我们已经有了,第二个是分组文件metadata,他至少由两列构成,一列是样本名称,一列是分组信息。 首先把样本名称变成数据框格式
分组信息包含在TCAG_id的第14,15字符很有用,他指示了样本是癌症还是癌旁或者是转移 病灶 官网解释如下,01-09是癌症,10-19是正常,20-29是癌旁
TCGA barcode的详细信息如下: 同时我们要注意,即使是肿瘤组织,01-09意义各不相同,比如,01代表原发灶,02代表转移灶,详细信息如下: 我们用table这个函数统计一下脑胶质瘤GBM样本的分类 table(substring(metadata$TCGA_id,14,15)) 有154个是原发灶,有13个是转移灶,很奇怪是吧,没有癌旁。但是这个是能理解的,人的大脑正常组织是有用的,不同于肝脏这类奇怪多一块少一块无所谓,切取大脑正常组织是没有伦理的。实际上TCGA里面还有一部分肿瘤是没有癌旁的,比如,淋巴瘤。 这一部分没有正常对照的肿瘤如何进行差异分析呢,一种方法是,使用GTEx数据库中的正常组织,这个我们留一个坑,以后再讲。 但是,今天我们的活还是要做,我们就用复发和非复发来区分即可。
此时metadata是这个样子的 构建dds的两个文件全部准备好,我们开始下一步 mRNA标准化这一步是为了代码复用,把counts文件统一命名 mycounts <> 构建dds对象,如果mycounts中的TCGA_id是行名,tidy这个参数设置为FASLE
Deseq2分析,这里面有很多步骤都自己运行了,这一步十分耗时,取决于样本数以及电脑内存大小,我的16g内存电脑运行5分钟,而我的学员们有的人要运行20个小时。甚至,如果,你分析的是乳腺癌,1000多个样本,小电脑根本过不去,此时,你可以考虑升级一下装备。 dds <>q(dds) 这个数据很重要,而且有些人获得也不容易,所以,需要保存一下,方便以后使用。
vst标准化,这一步跟上一步一样,速度取决于样本量和电脑 vsd - vst(dds, blind = FALSE) 为什么选择vst呢?看这个 这时候,Deseq2还内置了主成分分析来看一下样本分布
从图上我们可以看出,原发灶和转移灶,并不能完美分开,生物学意义就是,转移灶不是新的类型的肿瘤,他实际上还是脑胶质瘤,后续可能发生的结果是,下游额差异分析接结果不好,可能的解决方法是,找出配对的原发灶和转移灶来分析。我们看结果来说话。 获取标准化后的数据,这一步还会自动过滤掉不符合规定的基因,这时候,数据明显被标准化了 mRNA_exprSet_vst - as.data.frame(assay(vsd)) 保存一下这个数据,调整一下格式,就可以用于本文开头说的那一系列操作。
差异分析这里用到前面保存的dds,使用results函数提取 res - results(dds, tidy=TRUE) 我们看到这个数据,有foldchange值,有pvlaue,那么筛选差异基因,热图,火山图,GO,KEGG分析,GSEA分析就顺理成章啦。 |
|