转录组数据分析笔记（7）——DESeq2差异分析

前面处理好raw count定量表达矩阵，建立样本列表矩阵后，我们就可以在rstudio里运行DESeq2包进行差异基因筛选了。代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

library("DESeq2")
mycounts <- read.csv("gene_count_matrix_0_1.csv",row.names = 1)
mycounts_1 <- mycounts[rowSums(mycounts) != 0,] # 重新定义数据集，过滤mapping数为0的基因
mymeta <- read.csv("sample_list_0_1.csv",stringsAsFactors = T)  # 载入样本分组文件，遇到字符串将其转化为因子
colnames(mycounts_1) == mymeta$id   # 检查导入的两个数据集是否匹配，返回值为F需要重新匹配
mymeta$index <- factor(mymeta$index,levels = c("0","1"))    # 把样本分组文件的分组列转换到因子，两两比对把对照组放前面！
dds <- DESeqDataSetFromMatrix(countData = mycounts_1,
                              colData = mymeta,
                              design = ~index)  #构造用于差异表达分析的数据集
dds <- DESeq(dds)
res <- results(dds)
res_1 <- data.frame(res)    # 结果res不是常规的数据，需要转化成数据框
library("dplyr")
res_1 %>%   # dplyr给数据集增加新列
  mutate(group = case_when(
    log2FoldChange >=1 & padj <=0.05 ~ "UP",
    log2FoldChange <=-1 & padj <=0.05 ~ "DOWN",
    TRUE ~ "NOT_CHANGE"
  )) -> res_2
table(res_2$group)
write.csv(res_2,file = "gene_0_1.csv",
          quote = F)    # 输出文件

详细解释一下过程：

在R里运行程序或者写代码，首先要确定好工作目录在哪里，将之前Stringtie转化的定量表达矩阵和样本列表矩阵全都放在工作目录下，这里我的表达量矩阵是transcript_count_matrix_0_1.csv，分组列表矩阵是sample_list_0_1.csv。getwd()可以查看当前工作目录，在全局设置里可以更改工作目录。

library("DESeq2") # 加载DESeq2这个R包

mycounts <- read.csv("gene_count_matrix_0_1.csv",row.names = 1) # 载入raw count矩阵，以第一列数据作为行名，读取的矩阵命名为mycounts

mycounts_1 <- mycounts[rowSums(mycounts) != 0,] # 过滤每一行mapping总数为0的基因，将数据集整理命名为mycounts_1

mymeta <- read.csv("sample_list_0_1.csv",stringsAsFactors = T) # 载入样品列表，遇到字符串将其转化为一个因子

colnames(mycounts_1) == mymeta$id # 检查raw count矩阵第一行是否与样品列表的id列是否一致（如下）。这个很重要，不一致跑DESeq2会报错。如果显示false就要调整

mymeta$index <- factor(mymeta$index,levels = c("0","1")) # 这一步同样重要，把样本分组文件的分组列转换到因子，不然会报错。我这里对照组是第0天，所以把“1”放在“0”之后，这里顺序需要特别说明！样本和定量矩阵的分组只要一一对应可以不排序，这里一定要分清楚哪个组和哪个组进行比较，否则会得出完全相反的结论！

1
2
3

dds <- DESeqDataSetFromMatrix(countData = mycounts_1,
                              colData = mymeta,
                              design = ~index) 

# 中间那一长串是DESeq2包里的函数，countData是raw count定量矩阵，colData是样品列表，design是分组信息，这步是为了构造用于差异表达分析的数据集，并将数据集命名为dds

dds <- DESeq(dds) # 分析的核心DESeq程序

res <- results(dds) # 将结果输出至res数据集

res_1 <- data.frame(res) # res不是常规的数据，我们可以用head和class命令查看一下（如下图），需要转化成常规的数据框格式才可以对其进行加减列等操作，转换格式后的数据集名字为res_1

library("dplyr") # 加载这个包是为了对数据框进行操作，我是要增加新的一列统计差异表达情况

1
2
3
4
5
6

res_1 %>%   
  mutate(group = case_when(
    log2FoldChange >=1 & padj <=0.05 ~ "UP",
    log2FoldChange <=-1 & padj <=0.05 ~ "DOWN",
    TRUE ~ "NOT_CHANGE"
  )) -> res_2

# 调用dplyr包给数据集增加新的一列group，log2FoldChange >=1，padj <=0.05，判断这个基因表达为上调，在log2FoldChange <=-1，padj <=0.05时判断这个基因表达为下调，其余情况为该基因表达情况不变。将结果输出到res_2数据集。

FoldChange表示两样品间表达量比值，是倍数变化，差异表达基因分析里，log2 fold change绝对值大于1为差异基因筛选标准。padj是调整后的p值，在p检验里，p值小于0.05是有显著差异的标志。

table(res_2$group) # 查看差异基因表达的结果，上调基因多少，下调基因有多少，不变的有多少

write.csv(res_2,file = "gene_0_1.csv", quote = F) # 输出和生成gene_0_1.csv文件，即为结果文件

我对“0”组样本（对照）和长日照1天，2天和3天这一共4组样本分别进行两两对比，做了基因表达差异分析（0_1表示0组和1组对比，1_2表示1组和2组对比，依次类推不再赘述），同时与原文献补充数据2中的差异分析结果做对比

0组与1组对比，187个基因上调，149个基因下调；原文57个基因上调，79个基因下调

0组与2组对比，51个基因上调，42个基因下调；原文21个基因上调，24个基因下调

0组与3组对比，142个基因上调，143个基因下调；原文107个基因上调，84个基因下调

1组和2组对比，26个基因上调，29个基因下调；原文3个基因上调，0个基因下调

1组和3组对比，46个基因上调，50个基因下调；原文8个基因上调，2个基因下调

2组和3组对比，11个基因上调，13个基因下调；原文2个基因上调，1个基因下调

点击这里下载原文基因表达差异总表（补充数据2）

3.1 分析

我做的差异分析总体趋势和文章类似，都是短日照组（SD组，也就是0组）与长日照1、2、3天组（LD组，分别对应123组）相比，有较多基因出现差异性表达；而长日照1、2、3天组之间相比差异表达基因较少，这样才合理，表明拟南芥茎尖分生组织在开花期长日照下，确实有不同的基因参与了光周期诱导的开花过程。如果要对这些差异表达的基因做更深入的分析（下游分析），我们就要比对GO库或者KEGG库进行代谢通路富集分析和注释。以后的笔记会说到。
至于为什么我得到的差异基因普遍比原文多，一个是差异分析之前用的分析软件不同，在过滤数据过程中我的条件比较松（只过滤了count数为0的基因），另外不同软件的组装、比对和计数的算法实现也不一样。另外，我们得到这个基因差异表达结果之后，还可以做个更直观的火山图看我们的差异基因分布是否合理，之后也会说。