差异表达分析之Deseq2

介绍

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数据较大,为了避免重复计算浪费资源和时间,通常将结果保留至本地;
  • 设置过滤参数,主要有FDRlgFC等;
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 --------