edgeR/limma/DESeq2差异基因分析→ggplot2作火山图→biomaRt转换ID并注释

 

请一定看这里：写下来只是为了记录一些自己的实践，当然如果能对你有所帮助那就更好了，欢迎大家和我交流

三者区别

 

三者区别

差异分析流程:

1 初始数据

2 标准化(normalization)：DESeq、TMM等

为什么要标准化？
消除文库大小不同，测序深度对差异分析结果的影响
怎样标准化？
找到一个能反映文库大小的因子，利用这个因子对rawdata进行标准化

3 根据模型检验求p value：泊松分布(poisson distribution)、负二项式分布(NB)等

4 多重假设得FDR值：

检验方法：Wald、LRT
多重检验：BH

5 差异基因筛选：pvalue、padj

---

💥💥💥一、 edgeR的使用💥💥💥

因为目前没有合适的数据，所以数据来源于这里 参考这篇：刘尧科学网博客

0. 前期工作

1. 用到的gene.txt文件，内容如下

   

   gene.txt文件内容

   
   c1表示为一组，c2表示为另一组，.后为第几个样本
   读取数据并设置分组，要保证样本名称和分组名称的顺序是一一对应的。

#加载edgeR包
library(edgeR)
#读进来文件
targets <- as.matrix(read.delim("你的路径/gene.txt", sep = '\t', row.names = 1))
#分组，这里是两组，每组5个样本
group <- rep(c('c1','c2'),each = 5)

1. 构建DGEList对象

根据基因表达量矩阵以及样本分组信息构建DGEList对象

dgelist <- DGEList(counts = targets, group = group)

2. 过滤低表达的基因

DESeq2能够自动识别这些低表达量的基因的，所以使用DESeq2时无需手动过滤。
edgeR推荐根据CPM（count-per-million）值进行过滤，即原始reads count除以总reads数乘以1,000,000,使用此类计算方式时，如果不同样品之间存在某些基因的表达值极高或者极低，由于它们对细胞中分子总数的影响较大（也就是公式中的分母较大), 有可能导致标准化之后这些基因不存在表达差异，而原本没有差异的基因在标准化之后却显示出差异

这里参考这篇：当我们在说RNA-seq reads count标准化时，其实在说什么？
为了解决上述问题，BSM(between-sample normalization)类分出control set去评估测序深度而不是用所有数据，主要分三种：
(1) TMM(trimmed mean of M-values): TMM是M-值的加权截尾均值，即选定一个样品为参照，其它样品中基因的表达相对于参照样品中对应基因表达倍数的log2值定义为M-值。随后去除M-值中最高和最低的30%，剩下的M值计算加权平均值，权重来自Binomial data的delta方法 (Robinson and Oshlack, 2010)。
(2) RLE (relative log expression):首先计算每个基因在所有样品中表达的几何平均值。然后再计算该值与每个样品的比值的中位数，也叫被称为量化因子scale factor (Anders and Huber 2010)。
(3) UQ (upper quartile): 上四分位数 (upper quartile, UQ)是样品中所有基因的表达除以处于上四分位数的基因的表达值。同时为了保证表达水平的相对稳定，计算得到的上四分位数值要除以所有样品中上四分位数值的中位数。
以上三种方法效果大同小异，通常比较流行的是TMM和DESeq normalization

CPM 按照基因或转录本长度归一化后的表达即 RPKM (Reads Counts Per Million)、FPKM (Fragments Per Kilobase Million)和 TPM (Trans Per Million)，推荐使用TPM######

1)直接选某个值过滤

keep <- rowSums(dgelist$counts) >= 50
dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]

2)利用cpm过滤

keep <- rowSums(cpm(dgelist) > 1 ) >= 2
dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]

实际的数据分析中，还需多加尝试，选择一个合适的过滤条件

3. 标准化

calcNormFactors()函数对数据标准化，以消除由于样品制备或建库测序过程中带来的影响。
这里选的是TMM标准化方法，还有其他的可以?calcNormFactors进行查看

TMM法

 

dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')

dgelist_norm

plotMDS()是limma包中的方法，绘制MDS图，使用无监督聚类方法展示出了样品间的相似性（或差异）。可据此查看各样本是否能够很好地按照分组聚类，评估试验效果，判别离群点，追踪误差的来源等。

plotMDS(dgelist_norm, col = rep(c('red', 'blue'), each = 3))

可以考虑一下出现较大偏差的原因

4. 估算离散值

负二项分布（negative binomial，NB）模型需要均值和离散值两个参数。
edgeR中，分组矩阵使用model.matrix()获得，并可以通过estimateDisp()估算离散值。

design <- model.matrix(~group)    #构建分组矩阵
dge <- estimateDisp(dgelist_norm, design, robust = TRUE) #估算离散值
 plotBCV(dge) #作图查看

design中用0和1表示是哪一组，比如第二列1表示属于c2组

❗需要注意，标识好0、1类型后，后面的差异分析将以分组1的基因表达量相较于分组0是上调还是下调为准进行统计。因此在本示例中，后续分析获得的基因表达量上调/下调均为分组c2相较于分组c1而言的。实际的分析中，切记不要搞反了。(有时会出现两组顺序相反的情况，还没找到方法怎么实现)

estimateDisp()实际上是个组合函数，可以一步得到多个计算结果，例如在上文我们使用分组矩阵design通过estimateDisp()估算了3个值，其实就是estimateGLMTagwiseDisp()、estimateGLMCommonDisp()和estimateGLMTrendedDisp()这3个结果的组合。如果不指定分组矩阵，则将得到estimateCommonDisp()和estimateTagwiseDisp()的结果组合。

一定要记住是谁较谁

5. 差异基因分析

前面都是准备工作，现在可以开始正式分析了！

1) 负二项式广义对数线性模型（edgeR）

首先拟合负二项式广义对数线性模型（negative binomial generalized log-linear model），获取差异基因。这种方法大致可以这样理解，如果某个基因的表达值偏离这个分布模型，那么该基因即为差异表达基因。

使用edgeR包中的函数glmFit()和glmLRT()实现，其中glmFit()用于将每个基因的read count值拟合到模型中，glmLRT()用于对给定系数进行统计检验。

fit <- glmFit(dge, design, robust = TRUE)     #拟合模型
lrt <- glmLRT(fit)   #统计检验
topTags(lrt) #查看前10行 -n可修改查看前几行

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmLRT.csv', quote = FALSE) #输出主要结果
dge_de <- decideTestsDGE(lrt, adjust.method = 'fdr', p.value = 0.05)  #查看默认方法获得的差异基因
summary(dge_de)
plotMD(lrt, status = dge_de, values = c(1, -1), col = c('blue', 'red'))  #作图观测
abline(h = c(-1, 1), col = 'gray', lty = 2) #在图后加辅助线

decideTestsDGE() 的结果

decideTestsDGE()可用于统计差异基因数量。屏幕输出了其默认值（供参考，大多数情况下我们还是优先根据Fold Change值以及p值等手动去筛选，而不会在意这个程序自己判断的数值）
-1表示下调基因数量，1表示上调基因数量，0表示无差异基因数量。注意，对于这里的示例数据，基因表达量上调/下调均为“分组c2”相较于“分组c1”而言的。

输出的 glmLRT.csv

logFC即log2转化后的 Fold Change值，但是要注意，这里不是简单的将基因的read count值直接对比，而是分别计算了基因在两组中的CPM值，然后据此计算的logFC
logCPM是log2转化后的CPM值
LR，似然比统计
PValue，差异表达的p值
FDR，FDR校正后的p值

最后结果，也可以画火山图

 

2) 类似然负二项式广义对数线性模型（edgeR）

对于类似然负二项式广义对数线性模型（quasi-likelihood negative binomial generalized log-linear model），使用edgeR包中的函数glmQLFit()和glmQLFTest()实现，同样地，glmQLFit()用于将每个基因的read count值拟合到模型中，glmQLFTest()用于对给定系数进行统计检验，如果某个基因的表达值偏离这个分布模型，那么该基因即为差异表达基因。
相较于上一个模型，作者提到这个方法更严格一些。当然实际分析中还得视情况考虑了。

fit <- glmQLFit(dge, design, robust = TRUE)        #拟合模型
lrt <- glmQLFTest(fit)    #统计检验
topTags(lrt) #查看默认前10行

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmQLFTest.csv', quote = FALSE)  #输出主要结果
dge_de <- decideTestsDGE(lrt, adjust.method = 'fdr', p.value = 0.05)  #查看默认方法获得的差异基因
summary(dge_de)

summary(dge_de)

plotMD(lrt, status = dge_de, values = c(1, -1), col = c('blue', 'red'))     #作图观测
abline(h = c(-1, 1), col = 'gray', lty = 2)

跟第一种方法只有细微差别，大部分都是一样的

3） 配对检验（edgeR）

除了拟合模型的方法外，在edgeR中还可使用exactTest()直接执行两组负二项分布count之间基因均值差异的精确检验。

dge_et <- exactTest(dge) #检验
topTags(dge_et)
write.csv(topTags(dge_et, n = nrow(dgelist$counts)), 'exactTest.csv', quote = FALSE) #输出主要结果

topTags(dge_et)

dge_de <- decideTestsDGE(dge_et, adjust.method = 'fdr', p.value = 0.05)   #查看默认方法获得的差异基因
summary(dge_de)

summary(dge_de)

detags <- rownames(dge)[as.logical(dge_de)]
plotSmear(dge_et, de.tags = detags, cex = 0.5)      #作图观测
abline(h = c(-1, 1), col = 'gray', lty = 2)

因limma包的plotMD()函数无法在此处适用，这里使用的作图函数plotSmear()是edgeR包中的方法
图中纵轴为log2 Fold Change值；横轴为log2 CPM值，反映了基因表达量信息；红色的点表示差异基因（未使用颜色进一步区分上调/下调），黑色的点为无差异基因。

结果是这样

 

4) voom线性建模（limma）

limma包可以说是处理RNA-seq数据上的“老大”了，功能强大自然无需多说。因此也很容易得知，limma包中同样提供了多种差异基因分析的方法，其中最常用的就是voom方法（请允许我这么称呼它）
我们仍可以基于上文前几步获得的预处理结果（DGEList对象、标准化数据、估算的离散值等），继续使用limma包voom方法来完成后续的差异基因分析

将read count数据转换为log2-counts per million（logCPM），通过估计均值-方差（mean-variance）关系并使用它来计算合适的observation-level weights，然后，数据就可以进行线性建模。好吧具体它怎么工作的咱也看不懂（voom参考文献来源）……不过搞懂它的分析流程，以及结果怎么解读，还是可以的

limma_voom <- voom(dgelist_norm, design, plot = TRUE)

limma_voom

fit <- lmFit(limma_voom, design)  #拟合
fit <- eBayes(fit)
topTable(fit, coef = 2)
write.csv(topTable(fit, coef = 2, number = nrow(dgelist$counts)), 'limma_voom.csv', quote = FALSE) #输出主要结果

topTable(fit, coef = 2)

---

💥💥💥二、DESeq2的简洁使用💥💥💥

参考这篇

很慢，可以下这个

devtools::install_github('mikelove/DESeq2@ae7c6bd')

如果跟已安装的包冲突的话，

remova.packages('xxx')
BiocManager::install('xxx')

开始：

library(DESeq)
x <- as.matrix(read.delim("你的路径/gene.txt", sep = '\t', header = T, row.names = 1))

分组，这里是两组，每组5个样本

group <- rep(c('c1','c2'),each = 5)

由于DESeq包要求接下来的count data必须要整数型，因此我们需要对数据进行取整,然后将数据x和分组信息group读入到cds对象中

database <- round(as.matrix(x))
cds <- newCountDataSet(database,group)

有生物学重复

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,"c1","c2")

部分有生物学重复，其实同上

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,"c1","c2")

没有生物学重复

cds <- estimateDispersions(cds, method="blind", sharingMode="fit-only" )
res <- nbinomTest(cds,"c1","c2")

查看符合阈值的基因

table(res$padj <0.05)
res <- res[order(res$padj),]
sum(res$padj<=0.01,na.rm = T)
write.csv(res,"路径")

res.csv结果是这样的

---

💥💥💥三、DESeq2的详细使用💥💥💥

参考这篇: DESeq2做差异分析

0. 一些前期准备：

“gene.txt”，是一个基因表达量数据矩阵，包含10列样本，10个样本中前5个样本属于control组（c），后5个样本属于treat组（t）

0.1 构建基因表达矩阵countdata，即读入数据 read.delim()

❗colData的行名要与countData的列名一致！！

gene <- read.delim('C:/Users/wang/Desktop/gene.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

剔除全为0值的行

all <- apply(gene, 1, function(x) all(x==0) )
newdata <- gene[!all,]

指定分组因子顺序：差异基因分析需要指定比较分组的先后顺序，以确定谁相对于谁的表达量上调/或下调。
···第一种方式是在读取分组文件后，将分组列转换为因子类型（factor），并指定因子（分组）顺序，因子顺序指定对照在前处理在后；
···第二种方式是在后续使用results()获取差异结果时，指定比较的分组（推荐这种）
注意要保证表达矩阵中的样本顺序和这里的分组顺序是一一对应的，前5列为一组，后5列为一组

0.2 构建colData,

❗colData的行名要与countData的列名一致！！

colData <- data.frame(group = factor(rep(c('control', 'treat'), each = 5)))
colData <- data.frame(row.names=colnames(gene), colData)

两者的内容，参考这篇(https://www.jianshu.com/p/3a0e1e3e41d0)

1. 构建 DESeqDataSet对象，标准化reads count值，并用于存储输入值、中间计算和差异分析的结果

1.1 构建 DESeqDataSet 对象 dds = DESeqDataSet Object

①预处理，将所有样本基因表达量之和小于1的基因过滤掉（这步？）

dds <- dds[ rowSums(counts(dds))>1, ]

②差异分析

dds <- DESeqDataSetFromMatrix(countData = gene, colData = colData, design = ~group)

1.2 查看归一化后的 count 值分布

boxplot(log10(assays(dds)[['cooks']]), range = 0, las = 2)
plotDispEsts(dds)

 

cooks距离，详见http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

 

但是报错了

 

我看了一下这个显示NULL

第二天先运行了1.3的内容后，再运行这里就可以了，不明原因 :-(

boxplot()结果

 

plotDispEsts(dds)

 

1.3 vst标准化，获取归一化的基因表达矩阵norm_matrix, blind = FALSE指定实验设计不直接用于转换

vsd <- assay(vst(dds, blind = FALSE))
head(vsd, 10)
write.table(vsd, 'norm_matrix.txt', sep = '\t', col.names = NA, quote = FALSE)

vsd

2. 差异基因分析

之后直接运行默认的DESeq2差异分析流程就可以了
函数DESeq()是一个包含因子大小估计、离散度估计、负二项模型拟合、Wald统计等多步在内的过程，结果将返回至DESeqDataSet对象。这步比较耗时，特别是数据量较大时，新旧版DESeq2的运算效率差距极为明显
通过result()可获得最终计算的log2倍数变化和校正后p值等信息
①contrast参数用于指定比较的分组顺序，即谁相对于谁的表达量上调/或下调
②pAdjustMethod设定p值校正方法
③alpha为显著性水平，这里0.05为校正后p值小于0.05即为显著

2.1 标准方法

dds <- DESeq(dds, parallel = FALSE)   #若 parallel = TRUE 将启用多线程模式
suppressMessages(dds)
res <- results(dds, contrast = c('group', 'treat', 'control'), pAdjustMethod = 'fdr', alpha = 0.05)
write.csv(res, "你的路径/res.csv")

summary(res)

plotMA(res, alpha = 0.05, ylim = c(-3, 3))

dds过程

 

suppressMessages(dds)

通过summary()，可以根据预先设定的校正后p值<0.05水平（alpha=0.05，由results()指定），输出所比较两组间的上调/下调基因数量。这个结果可供参考，在后续也可以自己根据log2FC和校正后p值自定义作筛选

summary()

 

plotMA()

 

到这儿我发现我和例子中的结果有些差别了，但是还没找到原因，先过完流程吧 :-(

2.2 an alternate analysis: likelihood ratio test 似然比检验

ddsLRT <- DESeq(dds, test = 'LRT', reduced = ~ 1)
suppressMessages(ddsLRT)
resLRT <- results(ddsLRT, contrast = c('group', 'treat', 'control'), pAdjustMethod = 'fdr', alpha = 0.05)
write.csv(resLRT, "你的路径/ .csv")

summary(resLRT)

plotMA(resLRT, alpha = 0.05, ylim = c(-3, 3))

差异分析结果保存在res中，可通过as.data.frame()直接转化为数据框类型
包含了基因id、标准化后的基因表达值的平均值baseMean、log2FoldChange值、显著性p值pvalue以及校正后p值padj等主要信息

res结果

 

可以先大概看一些差异基因的数目

table(res$padj<0.05)

table()

2.3 可以先按校正和 p 值由小到大排个序，方便查看

deseq_res <- as.data.frame(res[order(res$padj), ])

将行名写在gene_id列中，这个时候它是最后一列

deseq_res$gene_id <- rownames(deseq_res)

先输出第7列，再输出前6列

write.table(deseq_res[c(7, 1:6)], '你的路径/DESeq2-test.txt', row.names = FALSE, sep = '\t', quote = FALSE)

最后的结果

3 ggplot2对差异基因作图

3.1 读进来最后的差异基因结果并进行分类

library(ggplot2)
deseq_res <- read.delim('你的路径 / DESeq2-test.txt', sep = '\t')

|log2FC| >= 1 & FDR p-value < 0.05 定为差异

deseq_res[which(deseq_res$padj %in% NA),'sig'] <- 'no diff'
deseq_res[which(deseq_res$log2FoldChange >= 1 & deseq_res$padj < 0.05),'sig'] <- 'up (p.adj < 0.05, log2FC >= 1)'
deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj < 0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)'
deseq_res[which(abs(deseq_res$log2FoldChange) < 1 | deseq_res$padj >= 0.05),'sig'] <- 'no diff'

也可以获取上调up /下调down 的差异表达基因（padjust < 0.05，并且|log2(foldchange)|>1)

diff_up = subset(deseq_res,padj < 0.05 & (log2FoldChange > 1))
write.csv(diff_up,file="diff_up.csv",row.names = F)

diff_down = subset(deseq_res,padj < 0.05 & (log2FoldChange > 1))
write.csv(diff_down,file="diff_down.csv",row.names = F)

3.2 画火山图

①纵轴为-log10(pvalue)，横坐标为log2FoldChange，差异基因展示为不同颜色

volcano_p <- ggplot(deseq_res, aes(log2FoldChange, -log(padj, 10))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1) +
  scale_color_manual(values = c('blue2', 'gray30', 'red2')) +
  theme(panel.grid = element_blank(), 
      panel.background = element_rect(color = 'black', fill = 'transparent'), 
      legend.position = c(0.26, 0.92)) +
  theme(legend.title = element_blank(), 
        legend.key = element_rect(fill = 'transparent'), 
        legend.background = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) +
  geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.25) +
  labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) +
  xlim(-5, 5)

volcano_p
ggsave('你的路径/volcano_p.png', volcano_p, width = 5, height = 6)

sig映射到color中
背景中fill = 'transparent'，使背景变为透明色
geom_vline和geom_hline为在x轴和y轴添加辅助线
labs在x轴和y轴添加横纵坐标名称
xlim限定x轴的显示范围
ggsave保存图片，或者可以直接Export

火山图

 

②纵轴为log2FoldChange，横坐标展示为标准化后的基因表达量的平均值 Average log10 baseMean，差异基因用不同颜色显示

volcano_count <- ggplot(deseq_res, aes(y = log2FoldChange, x = log10(baseMean))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1) +
  scale_color_manual(values = c('blue2', 'gray30', 'red2')) +
  theme(panel.grid = element_blank(), 
                   panel.background = element_rect(color = 'black', fill = 'transparent'), 
                   legend.position = c(0.2, 0.9)) +
  theme(legend.title = element_blank(), 
                   legend.key = element_rect(fill = 'transparent'), 
                   legend.background = element_rect(fill = 'transparent')) +
  geom_hline(yintercept = c(-1, 1), color = 'gray', size = 0.25) +
  labs(y = 'log2 Fold Change', x = 'Average log10 baseMean') +
  ylim(-5, 5)

volcano_count

ggsave('你的路径/volcano_count.png', volcano_p, width = 5, height = 6)

火山图

4 用biomaRt注释基因

参考这篇

4.1 我们利用useMart()函数选择“ENSEMBL_MART_ENSEMBL”，并将其赋值给my_mart对象

library('biomaRt')
library("curl")

my_mart <-useMart("ensembl")

在ensembl数据库中包含了77个数据集，可用下面这样的方式查看

datasets <- listDatasets(my_mart)
View(datasets)

datasets

4.2 选择一个数据集datasset，这里选人类的

my_dataset <- useDataset("hsapiens_gene_ensembl",
                         mart = my_mart)

4.3 💥根据ensembl ID获取基因名、描述或染色体信息

💥💥💥这里前半部分有误！请一定往下看解决办法

my_newid <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description","chromosome_name"),
                  filters = "ensembl_gene_id",
                  values = newinput,
                  mart = my_dataset)

 

image.png

这里一直报错，并且输出的为内容为0行
找到原因是：EBI数据库没有小数点，所以需要进一步替换为整数的形式。需要把小数点去掉！！这个很重要，所以需要加一个步骤

 

①还是将差异文件的行名提取出来

inputdata <- as.data.frame(row.names(deseq_res))

②这里将匹配到的.以及后面的数字连续匹配并替换为空，并重新赋值，一定要是data.frame格式

newinput <- as.data.frame(gsub("\\.\\d*", "", inputdata[,1]))

③getBM()转换ID

1）attributes参数：用来指定输出的数据类型，就是你要什么，比如entrezgene，hgnc_id。忘记的话可以用listAttributes(你自定义的dataset)
2）filters参数：用来指定数据的输入类型，比如你的原始信息是基因的ensembl ID，并且有这些基因的染色体位置信息，那么此处的filter就是ensembl_gene_id和chromosome_name等。
3）values参数：就是你待转换ID的数据
4）mart参数：此前定义的数据库，此处就是my_dataset

那么在我这里：
attributes ：我想要输出"ensembl_gene_id",转换后的"external_gene_name",转换后的"description",还可以有"chromosome_name"
filters：我的原始数据"ensembl_gene_id"
mart：之前建立的数据库

listAttributes(你的dataset) 可以查看可供选择的attributes

listAttributes(my_dataset)

my_result <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description"),
                  filters = "ensembl_gene_id",
                  values = newinput,
                  mart = my_dataset)

 

ID转换成功后

这样就完成了对ensembl_id的转化和注释

 

4.4 最后需要把结果文件deseq_res和注释文件my_result两者merge起来

merge需要有相同的gene_id
💥但是一定要看看自己文件里的gene_id是不是一致，如果有一个为小数，就要再添加一列取整后的gene_id

① deseq_res中gene_id有小数点 所以再加一列变成new_deseq_res，新增加的列名为gene_new_id

new_deseq_res <- as.data.frame(deseq_res)
new_deseq_res$gene_new_id <- gsub("\\.\\d*", "", deseq_res$gene_id)

② 修改一下列名,把含有小数点的列命名为gene_all_id,取整后的为gene_id，这一步是为了方便merge

colnames(new_deseq_res) <- c('baseMean', 'log2FoldChange','lfcSE','stat','pvalue','padj','gene_all_id','gene_id')

new_deseq_res

③ merge两个文件，即new_deseq_res和my_resullt，生成final_res文件

by = intersect(names(x), names(y)) 为取两个文件所有列名中列名相同的那列！

final_res <- merge(my_result, new_deseq_res, by = intersect(names(my_result), names(new_deseq_res)))
write.table(final_res, 'C:/Users/wang/Desktop/final_result.txt',row.names = FALSE, sep = '\t', quote = FALSE)

结果文件

4.5 还可以找到某个基因所在的通路GO号

参考这篇

① 选出要查找的基因

#举个例子
entrez = c("673", "837")

② 利用ensembl构建my_mart和my_dataset

my_mart <-useMart("ensembl") 

#`listDatasets()`可以查看可用的`datasets`
datasets <- listDatasets(my_mart)
View(datasets)

#构建`my_dataset`
my_dataset <- useDataset("hsapiens_gene_ensembl",
                         mart = my_mart)

③ 查看可输出的attributes

listAttributes(my_dataset)

④ 查找GOid

GOid <- getBM(attributes = c('entrezgene_id', 'go_id'),
              filters = 'entrezgene_id',
              values = entrez,
              mart = my_dataset)

结果

4.6 与4.5相反，可以通过所在的通路GO号找到某个基因

getBM(attributes = c('entrezgene_id', 'ensembl_gene_id'),
      filters = 'go',
      values = 'GO:0005524',
      mart = my_dataset)