转载果子学生信 https://mp.weixin.qq.com/s/Ph1O6V5RkxkyrKpVmB5ODA
前面我们从GDC下载了TCGA肿瘤数据库的数据,也能够把GDC下载的多个TCGA文件批量读入R
今天我们讲一下TCGA数据的标准化，以及差异分析，得到了标准化后的数据，我们就可以按照以前的帖子，做一系列操作
在得到了差异分析的结果后，我们可以完成热图，火山图，GO分析，KEGG分析，GSEA分析，就跟这个帖子中的一样。
来完成你的生信作业，这是最有诚意的GEO数据库教程
下面开始今天的教程：
首先加载上一次课获得的数据；
### 加载数据
load("expr_df.Rdata")
现在的数据是这个样子的
处理前
去掉ensemble ID的点号
library(tidyr)
expr_df_nopoint <- expr_df %>%
tidyr::separate(gene_id,into = c("gene_id"),sep="\.")
现在的数据是这个样子的
处理后
去掉点号，是为了用gtf文件。
gtf文件的获取和作用在这里
GTF文件有什么用啊？别的不谈，最起码能提lncRNA
加载gtf文件,这是目前我们能接触的最大文件，有260万行。
load(file = "gtf_df.Rda")
提取mRNA
mRNA_exprSet <- gtf_df %>%
dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #筛选gene,和编码指标
dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%
tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")
最终得到19668行，这是编码基因的个数，现在的数据是这个样子的
编码RNA
提取lncRNA
这里很有争议，而我的理由是，即使是编码基因，也会出现非编码转录本，而长链非编码RNA，指的是转录本，所以不能用gene的编码与否来界定
ncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA",
"processed_transcript","sense_intronic",
"bidirectional_promoter_lncRNA","non_coding",
"antisense_RNA")
LncRNA_exprSet <- gtf_df %>%
dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA) %>% #注意这里是transcript_biotype
dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>%
dplyr::distinct() %>% #删除多余行
dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%
tidyr::unite(gene_id,gene_name,gene_id,transcript_biotype,sep = " | ")
最终得到25530个非编码转录本，数据是这个样子的
非编码RNA
数据标准化
标准化和差异分析都是用Deseq2这个包来完成，首先要构建dds对象，构建这个对象需要两个文件，第一是输入数据，我们已经有了，第二个是分组文件metadata，他至少由两列构成，一列是样本名称，一列是分组信息。
首先把样本名称变成数据框格式
metadata <- data.frame(TCGA_id =colnames(expr_df)[-1])
分组信息包含在TCAG_id的第14,15字符很有用，他指示了样本是癌症还是癌旁或者是转移 病灶
官网解释如下,01-09是癌症，10-19是正常，20-29是癌旁
Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29
TCGA barcode的详细信息如下：
同时我们要注意，即使是肿瘤组织，01-09意义各不相同，比如，01代表原发灶，02代表转移灶，详细信息如下：
我们用table这个函数统计一下脑胶质瘤GBM样本的分类
table(substring(metadata$TCGA_id,14,15))
有154个是原发灶，有13个是转移灶，很奇怪是吧，没有癌旁。但是这个是能理解的，人的大脑正常组织是有用的，不同于肝脏这类奇怪多一块少一块无所谓，切取大脑正常组织是没有伦理的。实际上TCGA里面还有一部分肿瘤是没有癌旁的，比如，淋巴瘤。
这一部分没有正常对照的肿瘤如何进行差异分析呢，一种方法是，使用GTEx数据库中的正常组织，这个我们留一个坑，以后再讲。
但是，今天我们的活还是要做，我们就用复发和非复发来区分即可。
sample <- ifelse(substring(metadata$TCGA_id,14,15)=="01","cancer","recur")
## 这里的factor是为了dds的需要
metadata$sample <- as.factor(sample)
此时metadata是这个样子的
构建dds的两个文件全部准备好，我们开始下一步
mRNA标准化
这一步是为了代码复用，把counts文件统一命名
mycounts <- mRNA_exprSet
构建dds对象，如果mycounts中的TCGA_id是行名，tidy这个参数设置为FASLE
dds <-DESeqDataSetFromMatrix(countData=mycounts,
colData=metadata,
design=~sample,
tidy=TRUE)
Deseq2分析，这里面有很多步骤都自己运行了，这一步十分耗时，取决于样本数以及电脑内存大小，我的16g内存电脑运行5分钟，而我的学员们有的人要运行20个小时。甚至，如果，你分析的是乳腺癌，1000多个样本，小电脑根本过不去，此时，你可以考虑升级一下装备。
dds <- DESeq(dds)
这个数据很重要，而且有些人获得也不容易，所以，需要保存一下,方便以后使用。
save(dds,file="mRNA_exprSet_dds_sample.Rdata")
vst标准化,这一步跟上一步一样，速度取决于样本量和电脑
vsd <- vst(dds, blind = FALSE)
为什么选择vst呢？看这个
转录组的高级分析前该如何标准化数据？
Deseq2标准化的原理是什么，youtube上的StatQuest小哥视频说的特别好，可以看看这个帖子
DESeq2的标准化方法
这时候，Deseq2还内置了主成分分析来看一下样本分布
plotPCA(vsd, "sample")
从图上我们可以看出，原发灶和转移灶，并不能完美分开，生物学意义就是，转移灶不是新的类型的肿瘤，他实际上还是脑胶质瘤，后续可能发生的结果是，下游额差异分析接结果不好，可能的解决方法是，找出配对的原发灶和转移灶来分析。我们看结果来说话。
获取标准化后的数据,这一步还会自动过滤掉不符合规定的基因，这时候，数据明显被标准化了
mRNA_exprSet_vst <- as.data.frame(assay(vsd))
标准化之后
保存一下这个数据，调整一下格式，就可以用于本文开头说的那一系列操作。
save(ncRNA_exprSet_vst,file = "ncRNA_exprSet_vst.Rda")
差异分析
这里用到前面保存的dds,使用results函数提取
res <- results(dds, tidy=TRUE)
我们看到这个数据，有foldchange值，有pvlaue，那么筛选差异基因，热图，火山图，GO，KEGG分析，GSEA分析就顺理成章啦。