本篇学习转录组差异分析的代码，以TCGA的CHOL疾病为例。用较为基础的代码学习，以便理解。后续再学习更流畅的方法
上一篇中我们将数据下载并整理好后，现在将使用
DESeq2
、
edgeR
、
limma
这三大R包进行三次差异分析，会得到三组差异基因，最后再取交集。这样做的理由详见
转录组数据挖掘
一文
在正式介绍代码之前，先引入一个概念：函数黑匣子
函数黑匣子
如图1，想象一下这是3个盒子，里面装的就是下文要使用的
DESeq2
、
edgeR
、
limma
，只不过下文中还会添加一些细节进去
这些盒子，就是所谓的黑匣子。蓝色框是我们输入的内容exp和Group；红色框是我们得到的内容deg1、deg2、deg3；而从输入到输出，这中间发生了什么，即这些代码背后的统计学原理，就不是一眼能看明白的。（不知道背后的原理，确实不妨碍往下分析，不过若有条件还是深究一下为好）
图1
所以下面的代码几乎不用改动，我们只用输入exp和Group，然后全选运行即可。即使如此，有些细节还是需要注意，请往下看
三大R包差异分析
读取整理好的数据
load("TCGA-CHOL.Rdata")
table(Group)
1.deseq2
这部分的函数运行时间非常长，如果数据本身就很大，又要多次运行，就会很耗费时间。于是有了下面那段
if()
代码，能够将结果保存到一个RData文件里面，下次就不用运行这段函数，直接读取RData文件就行
注意事项：如果因表达矩阵有改动而需要再次运行，此时要从工作目录中删除
if()
生成的RData文件，因为
file.exists()
没法做到文件覆盖
#deseq2----
library(DESeq2)
colData condition=Group)
if(!file.exists(paste0(proj,"_dd.Rdata"))){
dds countData = exp,
colData = colData,
design = ~ condition)
dds save(dds,file = paste0(proj,"_dd.Rdata"))
}
# 用函数DESeqDataSetFromMatrix()将数据组织得符合要求，并完成了deseq2差异分析
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res "condition",rev(levels(Group))))
#constrast
c("condition",rev(levels(Group)))
class(res)
DEG1 library(dplyr)
DEG1 DEG1 = na.omit(DEG1)
head(DEG1)
#添加change列标记基因上调下调
logFC_t = 2
pvalue_t = 0.05
k1 = (DEG1$pvalue $log2FoldChange k2 = (DEG1$pvalue $log2FoldChange > logFC_t);table(k2)
DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG1$change)
head(DEG1)
2.edgeR
#edgeR----
library(edgeR)
dge dge$samples$lib.size $counts)
dge
design
dge dge dge
fit fit
DEG2=topTags(fit, n=Inf)
class(DEG2)
DEG2=as.data.frame(DEG2)
head(DEG2)
k1 = (DEG2$PValue $logFC k2 = (DEG2$PValue $logFC > logFC_t)
DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
head(DEG2)
table(DEG2$change)
3.limma
library(limma)
dge dge design v "quantile")
fit fit= eBayes(fit)
DEG3 = topTable(fit, coef=2, n=Inf)
DEG3 = na.omit(DEG3)
k1 = (DEG3$P.Value $logFC k2 = (DEG3$P.Value $logFC > logFC_t)
DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG3$change)
head(DEG3)
tj = data.frame(deseq2 = as.integer(table(DEG1$change
