Hua's Cabin

做好当下，不犹未来

差异表达分析之Deseq2

Posted on In , Valine:
Symbols count in article: 6.6k Reading time ≈ 6 mins.

介绍

Deseq2包的差异分析方法是基于负二项广义线性模型(negative binomial generalized linear models)；数据离散度和对数倍数变化。

A basic task in the analysis of count data from RNA-seq is the detection of differentially expressed genes. The count data are presented as a table which reports, for each sample, the number of sequence fragments that have been assigned to each gene. Analogous data also arise for other assay types, including comparative ChIP-Seq, HiC, shRNA screening, mass spectrometry. An important analysis question is the quantification and statistical inference of systematic changes between conditions, as compared to within-condition variability. The package DESeq2 provides methods to test for differential expression by use of negative binomial generalized linear models; the estimates of dispersion and logarithmic fold changes incorporate data-driven prior distributions

原理

文库均一化问题

数据处理时候遇到的样本表达值的均一化问题，如测序深度不同导致表达值不同。通常使用RPKM、FPKM和TPM等方法，但是DESeq2则是使用TMM方法均一化不同样本reads count数目。前三类方法均无法处理以下两种情况。

  • 测序深度不同造成不同样本的表达值不一样(reads depth)；
  • 不同组织存在特异的基因表达(library compostion)；

均一化步骤

DESeq2标准化的过程会解决上述提到两个问题。具体步骤如下：

  1. 对reads count去自然对数；取对数有助于降低离群值对数据整体影响，当然也存在0值的结果是Inf值；

  2. 求所有样本中，相同基因对数的均值：对同一基因在所有样本的数值取均值（先取对数，后用对数求均值的方法叫几何平均数）；

  3. 去掉Infinity: 第二步计算几何平均数时候，含有-inf的基因其最后的均值仍然是-inf，该步骤最后的目的是剔除样本中reads数据为0的基因。

    如果我们比较的是肝脏和脾脏的转录组，按这种剔除方法，我们会把那些所有只在肝脏（或脾脏）中转录的基因都给剔除掉，从理论上讲，最终剩下的基因基本上就是管家基因（house keeping genes）了，也就是说在不同组织类型中都表达相似的基因

  4. 矩阵减均值：上一步去除了Inf值，这一步把每个样本中的每个值减去该基因在所有样本中的均值，这里把差值可以理解成偏离值；

  5. 计算每个样本的中位数：计算每个样本中所有基因的对数的中位数，中位数均有较好的平衡离群值的作用，上一步的均值会受到离群值较大的影响即使我们做了对数转换。

    这里使用中位数主要是为了排除一些极端表达基因的影响，极端表达指的是表达量极高或极低，它们能够对均值造成很大的影响，而那些表达量差异极大的基因对于中位数的影响，并不比那些表达量差异较小的基因对中位数的影响大。再加上那些表达量差异极大的基因数量通常情况下很少，因此，更多情况下，中等程度表达差异的基因和管家基因对中位数的影响更大一些.

  6. 将中位数转换为真数，计算每个样本最终的标准化因子: 将中位数转换成真数，即$e^{median}$值;

  7. 原始reads count数除以标准化因子： 把原始reads数目初一上一步得到的标准化因子数目，后续得到的商是整型，采取四舍五入的方法取整。

总结

  • log转换会剔除那些只在某个样本类型中的转录的所有基因（例如肝脏与脾脏），这种处理也会消除原始reads数的异常值（通过几何均数）;

  • 中位数的处理会进一步降低那些高数值reads数基因的影响，从而关注那些中等表达程度的基因;

实例

该部分是从安装到分析的实操解读，使用Rstudio和R的组合。

安装Deseq2包

安装该包最好使用如下命令，Deseq2包因为功能强大，因此它的依赖包也巨多，安装起来很费时间。

1
2
3
4
5
6
7
8
9
10
11
12
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")

# 查看源文件及帮助文档命令
browseVignettes("DESeq2")

# 加载本次分析使用到的R包
library(data.table)
library(dplyr)
library(tibble)
library(DESeq2)

数据准备

分析输入数据有两类

  • 一类数值是count data的样本表达谱数据，一般情况下colnames是样本ID，rownames是基因ID等，我们用countData代替；
  • 另一类是样本ID及其对应的标签分组信息等，我们用colData代替；

所需数据已经上传至百度网盘: https://pan.baidu.com/s/1Xy-i8ApdScuDiXqki0jLpg 提取码：q5ls

  • 这次使用的数据是miRNA表达谱数据，也即是行名从gene换成miRNA。
1
2
3
4
5
6
7
8
library(data.table)
library(dplyr)
library(tibble)
library(DESeq2)

# miRNA 
mirna.phe <- read.csv("miRNA_Phenotype_20200409.csv")
mirna.prf <- fread("miRNA_counts_20200409.csv")

实操运行

在处理数据的过程中，我们一般会先对这两类数据进行预处理，主要有如下操作：

  • 根据样本ID对两类数据取交集；
  • 去除低出现率和低表达的miRNA；

Deseq2分析主要有如下步骤：

  • 构建dds对象：
    • countData：非负整型数值的矩阵；
    • colData：带有样本分组信息的矩阵或者数据框；
    • design参数是样本ID的分组信息；
    • tidy参数判断矩阵第一列是否是列名；
  • DESeq函数的差异表达分析；
  • results函数获取结果信息；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
Deseq2fun <- function(phen, prof, occurence=0.5, ncount=10){
  # phen <- mirna.phe
  # prof <- mirna.prf
  # occurence <- 0.5
  # ncount <- 10
  
  sid <- intersect(phen$SampleID, colnames(prof))
  phe <- phen %>% filter(SampleID%in%sid) %>%
    arrange(SampleID)
  
  # quality control
  prf <- prof %>% column_to_rownames("gene_id") %>% 
    dplyr::select(as.character(phe$SampleID)) %>% 
    rownames_to_column("Type") %>% 
    filter(apply(dplyr::select(., -one_of("Type")), 1, function(x){sum(x != 0)/length(x)}) > occurence) %>%
            data.frame(.) %>% 
    column_to_rownames("Type")
  prf.cln <- prf[rowSums(prf) > ncount, ]
  
  # judge no row of profile filter
	if (nrow(prf.cln) == 0) {
		stop("No row of profile to be choosed\n")
	}
  
  # determine the right order between profile and phenotype 
  for(i in 1:ncol(prf.cln)){ 
    if (!(colnames(prf.cln) == phe$SampleID)[i]) {
      stop(paste0(i, " Wrong"))
    }
  }
  
  metadata <- phe %>% dplyr::select(SampleID, group1)
  countData <- prf.cln %>% rownames_to_column("gene_id")
  dds <- DESeqDataSetFromMatrix(countData=countData, 
                              colData=metadata, 
                              design=~group1, tidy = TRUE)
  
  dds <- DESeq(dds)
  res <- results(dds)
  res <- res[order(res$padj), ]
  return(res)
}

保存结果文件以及设置过滤参数

  • 一般countData数据较大，为了避免重复计算浪费资源和时间，通常将结果保留至本地；
  • 设置过滤参数，主要有FDR和lgFC等；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# cutoff or threshold 
FDR <- 0.05
lgFC <- 2

# run DEseq2
if(!file.exists("../../Result/Differential_expression/RNA_DEseq2.RData")){
  # miRNA
  miRNA.des2 <- Deseq2fun(mirna.phe, mirna.prf)
  # RNA 
  RNA.des2 <- Deseq2fun(rna.phe, rna.prf)
  save(miRNA.des2, RNA.des2, 
       file = "../../Result/Differential_expression/RNA_DEseq2.RData")
  
}else{
  load("../../Result/Differential_expression/RNA_DEseq2.RData")
}

# filter by FDR & lgFC
diff_miRNA_deseq2 <- subset(miRNA.des2, padj < FDR & abs(log2FoldChange) > lgFC)

结果解读

Deseq2的结果矩阵没有包含富集方向列信息，这是我认为最不好的地方，它仅仅将该信息以对象形式保存在结果文件中，但又因为后续以write.csv保存数据时会丢失该信息。head(miRNA.des2)的表头有分组比较以及检验方法信息，下方则是所有的miRNA比较结果:

1
2
3
4
5
6
7
head(miRNA.des2, 1)
#log2 fold change (MLE): group1 RB vs PB 
#Wald test p-value: group1 RB vs PB 
#DataFrame with 6 rows and 6 columns
#                        baseMean   log2FoldChange             lfcSE             stat
#                       <numeric>        <numeric>         <numeric>        <numeric>
#hsa-mir-106b    4948.33354145209 1.76406198826099 0.119754045761411 14.7307089046125

除此之外，还可以使用summary查看综合结果：此处未进行任何质控。

1
2
3
4
5
6
7
8
summary(miRNA.des2)
#out of 548 with nonzero total read count
#adjusted p-value < 0.1
#LFC > 0 (up)       : 230, 42%
#LFC < 0 (down)     : 144, 26%
#outliers [1]       : 0, 0%
#low counts [2]     : 0, 0%
#(mean count < 1)

质控并添加富集方向，最后保存结果成csv文件；

1
2
3
4
5
miRNA.res.origin <- miRNA.des2 %>% 
  as.data.frame() %>%
  rownames_to_column("miRNA_id") %>%
  mutate(Enrichment=ifelse(log2FoldChange > 0, "RB", "PB"))
write.csv(miRNA.res.origin, "../../Result/Differential_expression/miRNA_DEseq2_origin.csv", row.names = F)

另外如果是能注释到基因的mRNA等数据，还可以继续对gene_id进行注释，这又涉及到繁多的gene_id命名方式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
library(biomaRt)
library(curl)
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
# origin 
ensembl_gene_id.origin <- rownames(RNA.des2)
gene_symbols.origin <- getBM(attributes=c('ensembl_gene_id','external_gene_name',
                                   "transcript_biotype","description"),
                                    filters='ensembl_gene_id',
                                    values=ensembl_gene_id.origin, mart = mart)
# group1 RB vs PB
diff_RNA_symbol_deseq2.origin <- inner_join(
  data.frame(RNA.des2) %>% rownames_to_column("gene_id"),
  gene_symbols.origin, by=c("gene_id"="ensembl_gene_id")) %>%
  mutate(Enrichment=ifelse(log2FoldChange > 0, "RB", "PB")) 

write.csv(diff_RNA_symbol_deseq2.origin, "../../Result/Differential_expression/RNA_symbol_DEseq2_origin.csv", row.names = F)

后续分析

基于差异表达数据，主要的分析就是可视化这些结果：

  • 火山图；
  • 热图；

引用

  1. DESeq2 tutorial
  2. DESeq2包的基本原理

参考文章如引起任何侵权问题，可以与我联系，谢谢。


------------- The End Thanks for reading --------