文章目录
学习目标
DESeq2差异基因表达分析流程
第一步：估计大小因子
第二步：估计基因离散（gene-wise dispersion）
第三步：拟合曲线到基因的分散估计
第四步:将基因离散估计值向曲线预测值收缩
MOV10 DE分析:探讨离散估计和评估模型拟合

学习目标

理解使用DESeq2差异表达分析过程中的不同步骤
探讨离散度在差异表达分析中的重要性，并利用离散度值的图来探讨NB模型的假设

DESeq2差异基因表达分析流程
之前，我们使用适当的设计公式创建了DESeq2对象，并使用两行代码运行DESeq2:
# DO NOT RUN

## Create DESeq2Dataset object
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)

## Run analysis
dds <- DESeq(dds)

我们用DESeq2完成了差异基因表达分析的整个工作流程。分析中的步骤输出如下:



我们将详细查看每一个步骤，以更好地理解DESeq2是如何执行统计分析的，以及我们应该检查哪些度量来探究我们的分析质量。
第一步：估计大小因子
差异表达分析的第一步是估计大小因子，这正是我们已经对原始计数进行归一化所做的。



当执行差异表达分析时，DESeq2将自动估计大小因子。但是，如果你已经使用estimateSizeFactors()生成了大小因子，就像我们前面所做的那样，那么DESeq2将使用这些值。

为了归一化计数数据，DESeq2使用先前在“计数归一化”课上讨论的比值中位数方法计算每个样本的大小因子。

MOV10 DE分析:检查大小因子

让我们快速看一下每个样本的大小因子值:
## Check the size factors
sizeFactors(dds)

Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 
 1.1150371  0.9606366  0.7493552 
Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 
 1.5634128  0.9359082  1.2257749 
Mov10_oe_2 Mov10_oe_3 
 1.1406863  0.6541689

这些数字应该与我们在运行函数estimateSizeFactors(dds)时最初生成的数字相同。看看每个样本的读段总数:
## Total number of raw counts per sample
colSums(counts(dds))

Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 
  31160785   26504972   20498243 
Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 
  45300696   26745860   32062221 
Mov10_oe_2 Mov10_oe_3 
  30025690   17285401 

这些数字与大小因子的关系是什么?
我们看到较大的大小因子对应于序列深度较高的样本，这是有意义的，因为要生成归一化的计数，我们需要将计数除以尺寸因子。这就解释了样品间测序深度的差异。
现在看一看归一化后的总深度，使用以下代码:
## Total number of normalized counts per sample
colSums(counts(dds, normalized = TRUE))

Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 
  27945962   27591050   27354509 
Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 
  28975519   28577441   26156696 
Mov10_oe_2 Mov10_oe_3 
  26322477   26423452

样本间的值如何与每个样本的总计数进行比较?

你可能期望经过归一化后，各个样本的计数是完全相同的。然而，在归一化过程中，DESeq2也解释了RNA的组成。通过使用大小因子的中位数比值，DESeq2不应该偏向于由几个DE基因吸收的大量计数；然而，这可能导致大小因子与仅仅根据测序深度预期的大小因子有很大的不同。
第二步：估计基因离散（gene-wise dispersion）
差异表达分析的下一步是对基因离散的估计。在讨论细节之前，我们应该对DESeq2中离散指的是什么有一个很好的概念。



在RNA-seq计数数据中，我们知道:

为了确定差异表达基因，考虑到组内(重复样本之间)的差异，我们需要识别组间平均表达显著不同的基因。
组内(重复样本之间)的变异需要考虑到方差随平均表达值增加的事实，如下图所示(每个黑点是一个基因)。



为了准确地识别DE基因，DESeq2需要解释方差与均值之间的关系。 我们不希望所有的DE基因都是低计数的基因，因为低表达基因的方差更低。

DESeq2没有使用方差作为数据中变异的度量衡(因为变异与基因表达水平相关)，而是使用离散度（dispersion） 作为变异度量衡，它解释了基因的变异和平均表达水平。离散度的计算方法为Var = μ + α*μ^2，其中:α = dispersion, Var = variance, 和 μ = mean，有以下关系：





对离散的影响




方差增加
离散增加


平均表达增加
离散减少


对于具有中等到高计数值的基因，离散度的平方根将等于变异系数（ \sqrtα=\frac{σ}{μ} ）。因此，0.01的离散度意味着在生物重复中围绕平均数预期的10%的变异。具有相同平均数的基因的离散度估计将仅基于其方差而有所不同。因此，离散度估计（dispersion estimates）反映了给定平均值的基因表达的方差。在下面的图中，每个黑点都是一个基因，离散度与每个基因的平均表达量（组内重复）作对比。



离散和我们的模型有什么关系?

为了精确地建立测序计数模型，我们需要对每个基因的组内变异(同一样本组重复之间的变异)做出准确的估计。由于每组只有少量(3-6)重复，对每个基因变异的估计往往是不可靠的。

为了解决这个问题，DESeq2共享基因间的信息，通过一种叫做“收缩（shrinkage）”的方法，基于基因的平均表达水平，产生更准确的变异估计。DESeq2假设具有相似表达水平的基因应该具有相似的离散度。
分别估计每个基因的离散度:

DESeq2根据基因的表达水平(组内重复的平均计数)和方差估计每个基因的离散度。
第三步：拟合曲线到基因的分散估计
工作流程的下一步是拟合曲线到基因分散估计。拟合数据曲线背后的想法是，不同的基因会有不同程度的生物变异，但是，在所有的基因中，会有一个合理的分散估计分布。



下图中这条曲线以红线表示，它绘制了给定表达强度基因的预期离散的估计值。每个黑点都是一个基因，具有相关的平均表达水平和离散度的最大似然估计(MLE)(步骤1)。


第四步:将基因离散估计值向曲线预测值收缩
工作流程的下一步是将基因离散估计值向曲线预测值收缩。



当样本量较小时，曲线可以更准确地识别差异表达的基因，每个基因的收缩强度取决于:

基因的离散度距离红色曲线的距离
样本量（样本量越大，则收缩的越少）

这种收缩方法对于减少差异表达分析中的假阳性尤其重要。 具有低离散度估计数的基因向曲线收缩，并且更精确、更大的收缩值被输出用于模型拟合和差异表达检验。这些缩小后的估计值代表了组内的变异，需要确定是否组间的基因表达有显著差异。

稍微高于曲线的离散估计也向曲线收缩，以便更好地估计离散；然而，具有极高离散值的基因则不收缩。这是由于该基因可能没有遵循建模假设，并且由于生物学或技术原因比其他基因具有更高的变异性[1]。向曲线方向收缩估计值可能会导致假阳性，所以这些值不会被收缩。这些基因被下面的蓝色圆圈所包围。



这是一个很好的绘图，可以确保你的数据很好的拟合DESeq2模型。你希望你的数据通常散布在曲线周围，随着平均表达水平的增加，离散度逐渐减小。如果你看到一片云（cloud）或不同的形状，那么你可能想要更深入探索你的数据，看看你是否有污染(线粒体等)或离群样本。请注意，对于任何低自由度的实验，在plotdispest()图中整个平均值范围的收缩程度。

令人不安的离散曲线（worrisome dispersion plots） 的例子如下:

下面的图显示了分散值的云，它们通常不遵循曲线。这将令人担忧，并表明数据与模型不符。



下一张图显示，离散度值先减小，然后随着表达值的增大而增大。根据我们的预期，较大的平均表达式值不应该有较大的离散度——我们希望离散度随着均值的增加而减少。这表明高表达基因的变异比预期的要少。这也表明在我们的分析中可能存在一个离群样本或污染。


MOV10 DE分析:探讨离散估计和评估模型拟合
让我们来看看我们的MOV10数据的离散估计:
## Plot dispersion estimates
plotDispEsts(dds)



由于我们的样本量较小，对于许多基因我们看到了相当大的收缩。你认为我们的数据能很好地拟合这个模型吗?

我们看到随着平均表达的增加离散度很好地降低了，这很好。我们还可以看到离散估计通常围绕曲线，这也是预期的。总的来说，这个图看起来不错。我们确实看到了强烈的收缩，这可能是由于我们的一个样本组只有两个重复的事实。我们拥有的重复越多，应用于离散估计的收缩就越少，可以识别的DE基因就越多。为了更好地估计变异，我们一般建议每个条件至少有4个生物重复。

Exercise
Given the dispersion plot below, would you have any concerns regarding the fit of your data to the model?

If not, what aspects of the plot makes you feel confident about your data?
If so, what are your concerns? What would you do to address them?



Answer
# Ans: Yes, there are some concerns. The data does not scatter around the fitted curve, and the distribution of normalized counts are restricted in a small range. I would double check QC of my samples to make sure that there are no contamination or outliers.

差异基因鉴定



基因表达标准化

不同样品的测序量会有差异，最简单的标准化方式是计算 

counts per million (CPM)，即原始reads count除以总reads数乘以1,000,000。

这种计算方式的缺点是容易受到极高表达且在不同样品中存在差异表达的基因的影响；这些基因的打开或关闭会影响到细胞中总的分子数目，可能导致这些基因标准化之后就不存在表达差异了，而原本没有差异的基因标准化之后却有差异了。

RPKM、FPKM和TPM是CPM按照基因或转录本长度归一化后的表达，也会受到这一影响。

calc_cpm <- function (expr_mat, spikes = NULL){
  norm_factor <- colSums(expr_mat[-spikes,])
  return(t(t(expr_mat)/norm_factor)) * 10^6
}


为了解决这一问题，研究者提出了其它的标准化方法。

量化因子 (size factor,SF)是由DESeq提出的。其方法是首先计算每个基因在所有样品中表达的几何平均值。每个细胞的量化因子(size factor)是所有基因与其在所有样品中的表达值的几何平均值的比值的中位数。由于几何平均值的使用，只有在所有样品中表达都不为0的基因才能用来计算。这一方法又被称为RLE (relative log expression)。

calc_sf <- function (expr_mat, spikes=NULL){
  geomeans <- exp(rowMeans(log(expr_mat[-spikes,])))
  SF <- function(cnts){
    median((cnts/geomeans)[(is.finite(geomeans) & geomeans >0)])
  }
  norm_factor <- apply(expr_mat[-spikes,],2,SF)
  return(t(t(expr_mat)/norm_factor))
}


上四分位数 (upperquartile,UQ)是样品中所有基因的表达除以处于上四分位数的基因的表达值。同时为了保证表达水平的相对稳定，计算得到的上四分位数值要除以所有样品中上四分位数值的中位数。

calc_uq <- function (expr_mat, spikes=NULL){
  UQ <- function(x) {
    quantile(x[x>0],0.75)
  }
  uq <- unlist(apply(expr_mat[-spikes,],2,UQ))
  norm_factor <- uq/median(uq)
  return(t(t(expr_mat)/norm_factor))
}


TMM是M-值的加权截尾均值。选定一个样品为参照，其它样品中基因的表达相对于参照样品中对应基因表达倍数的log2值定义为M-值。随后去除M-值中最高和最低的30%，剩下的M值计算加权平均值。每一个非参照样品的基因表达值都乘以计算出的TMM。这个方法假设大部分基因的表达是没有差异的。



DESeq2 差异基因鉴定一步法

为了简化差异基因的运算，易生信做了脚本封装，DESeq2.sh，只需提供原始基因表达矩阵、样品分组信息表即可进行差异基因分析和鉴定。

DESeq2.sh
OPTIONS:
    -f  Data file [A gene count matrix, NECESSARY]
        CHECK ABOVE FOR DETAILS

    -s  Sample file [A multiple columns file with header line, 
        For <timeseries>,  one <conditions> columns is needed.
        NECESSARY]
        CHECK ABOVE FOR DETAILS

    -d  The design formula for DESeqDataSetFromMatrix.
        [Default <conditions>, 
        accept <cell+time+cell:time> for example 2.]

    -D  The reduced design formula for DESeq.
        [Only applicable to <timeseries> analysis, 
        accept <cell+time> or <time> or <cell> for example 2.]

    -m  Specify the comparasion mode.
        [Default <pairwise>, accept <timeseries>,
        <pairwise> comparasion will still be done in <timeseries>
        mode.
        NECESSARY]

    -p  A file containing the pairs needed to do comparasion. 
        CHECK ABOVE FOR DETAILS
        All samples will be compared in <pairwise> mode if not specified here.

    -F  Log2 Fold change for screening DE genes.
        Default 1

    -P  FDR for screening DE genes.
        Default 0.01

    -q  FDR for screening time-series DE genes.
        Default 0.1

    -e  Execute programs
        Default TRUE

    -i  Install packages of not exist. 
        Default FALSE

Eg.
    DESeq2.sh -f matirx -s sample -d conditions
    DESeq2.sh -f matirx -s sample -p compare_pair


具体使用

cd ~/data
# ehbio_trans.Count_matrix.xls和sampleFile都是前面用到的文件
DESeq2.sh -f ehbio_trans.Count_matrix.xls -s sampleFile

[1] "Perform pairwise comparasion using <design=~conditions>"
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
[1] "Output normalized counts"
[1] "Output rlog transformed normalized counts"
[1] "Performing sample clustering"
null device 
          1 
[1] "DE genes between untrt trt"
[1] "PCA analysis"
null device 


运行结束即可获得以下结果文件

# 具体的文件内容和图的样式见后面的分步法文档
# 原始输入文件
ehbio_trans.Count_matrix.xls
sampleFile

# 所有差异基因列表
ehbio_trans.Count_matrix.xls.DESeq2.all.DE

# PCA结果
ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.pca.pdf

# 样品相关性层级聚类结果
ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.pearson.pdf

# rlog转换后的标准化后的表达结果
ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.xls

# 标准化后的表达结果
ehbio_trans.Count_matrix.xls.DESeq2.normalized.xls

# 运行脚本
ehbio_trans.Count_matrix.xls.DESeq2.r

# 差异基因结果
ehbio_trans.Count_matrix.xls.DESeq2.untrt._higherThan_.trt.id.xls
ehbio_trans.Count_matrix.xls.DESeq2.untrt._higherThan_.trt.xls
ehbio_trans.Count_matrix.xls.DESeq2.untrt._lowerThan_.trt.id.xls
ehbio_trans.Count_matrix.xls.DESeq2.untrt._lowerThan_.trt.xls

# 火山图和火山图输入数据
ehbio_trans.Count_matrix.xls.DESeq2.untrt._vs_.trt.results.xls
ehbio_trans.Count_matrix.xls.DESeq2.untrt._vs_.trt.results.xls.volcano.pdf


一步法脚本可在https://ke.qq.com/course/301387?tuin=20cd7788或点击阅读原文获得。



DESeq2 差异基因鉴定分步法



安装包 DESeq2安装方法如下

source("https://bioconductor.org/biocLite.R")
biocLite('BiocInstaller')
biocLite("DESeq2")
install.packages(c("gplots", "amap", "ggplot2"))


A distributional assumption is needed because we want to estimate the 

probability of extreme events (large fold change just appearing by 

chance) from limited replicates. The negative binomial (a.k.a. 

Gamma-Poisson) is a good choice for RNA-seq experiments because


The observed data at gene level is inherently counts or estimated 

counts of fragments for each feature and
The spread of values among biological replicates is more than given 

by a simpler, one parameter distribution, the Poisson; and it seems 

to be captured by the NB sufficiently well




加载包

library(DESeq2)

## Loading required package: S4Vectors

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, cbind, colMeans,
##     colnames, colSums, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply,
##     lengths, Map, mapply, match, mget, order, paste, pmax,
##     pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
##     rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomicRanges

## Loading required package: GenomeInfoDb

## Loading required package: SummarizedExperiment

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following object is masked from 'package:base':
## 
##     apply

library("RColorBrewer")
library("gplots")

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:IRanges':
## 
##     space

## The following object is masked from 'package:S4Vectors':
## 
##     space

## The following object is masked from 'package:stats':
## 
##     lowess

library("amap")
library("ggplot2")
library("BiocParallel")

# 多线程加速, 使用4个核
# 如果你电脑性能强大，可以把这个数变大
register(MulticoreParam(4))




读入数据

# 注意文件的位置，程序自己不知道文件在什么地方
# 如果只给文件名，不给路径，程序只会在当前目录查找文件
# 若找不到，则会报错
# 可以使用下面命令设置工作目录到文件所在目录
# 或提供文件全路径
# setwd("~/data")
data <- read.table("ehbio_trans.Count_matrix.xls", header=T, row.names=1, com='', quote='',
    check.names=F, sep="\t")

# 撇掉在多于两个样本中count<1的值，如果样本数多，这个数值可以适当增加
# 排除极低表达基因的干扰
data <- data[rowSums(data)>2,]
head(data)

##                 untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011
## ENSG00000227232           13            25            23            24
## ENSG00000278267            0             5             3             4
## ENSG00000268903            0             2             0             0
## ENSG00000269981            0             3             0             1
## ENSG00000241860            3            11             1             5
## ENSG00000279457           46            90            73            49
##                 trt_N61311 trt_N052611 trt_N080611 trt_N061011
## ENSG00000227232         12          12          22          22
## ENSG00000278267          2           4           3           1
## ENSG00000268903          0           2           0           3
## ENSG00000269981          0           2           0           0
## ENSG00000241860          3           2           0           2
## ENSG00000279457         52          46          89          31

# 读入分组信息
sample <- read.table("sampleFile", header=T, row.names=1, com='',
    quote='', check.names=F, sep="\t", colClasses="factor")

sample <- sample[match(colnames(data), rownames(sample)),, drop=F]
sample_rowname <- rownames(sample)

# 下面的可以忽略，如果没遇到错误不需要执行
# 目的是做因子转换
sample <- data.frame(lapply(sample, function(x) factor(x, levels=unique(x))))
rownames(sample) <- sample_rowname
sample

##               conditions
## untrt_N61311       untrt
## untrt_N052611      untrt
## untrt_N080611      untrt
## untrt_N061011      untrt
## trt_N61311           trt
## trt_N052611          trt
## trt_N080611          trt
## trt_N061011          trt




产生DESeq数据集

DESeq2包采用DESeqDataSet存储原始的read count和中间计算的统计量。

每个DESeqDataSet对象都要有一个实验设计formula，用于对数据进行分组，以便计算表达值的离散度和估计表达倍数差异，通常格式为~ batch + conditions 

(为了方便后续计算，最为关注的分组信息放在最后一位)。

countData: 表达矩阵

colData: 样品分组信息表

design: 实验设计信息，conditions必须是colData中的一列

ddsFullCountTable <- DESeqDataSetFromMatrix(countData = data,
        colData = sample,  design= ~ conditions)

dds <- DESeq(ddsFullCountTable)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing




获取标准化后的数据

# ?counts查看此函数功能
# normalized=T, 返回标准化的数据
normalized_counts <- counts(dds, normalized=TRUE)
head(normalized_counts)

##                 untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011
## ENSG00000227232    12.730963     21.179286    19.4979982     25.994727
## ENSG00000278267     0.000000      4.235857     2.5432172      4.332454
## ENSG00000268903     0.000000      1.694343     0.0000000      0.000000
## ENSG00000269981     0.000000      2.541514     0.0000000      1.083114
## ENSG00000241860     2.937915      9.318886     0.8477391      5.415568
## ENSG00000279457    45.048024     76.245430    61.8849507     53.072567
##                 trt_N61311 trt_N052611 trt_N080611 trt_N061011
## ENSG00000227232  13.423907   17.885811   15.750664   23.250144
## ENSG00000278267   2.237318    5.961937    2.147818    1.056825
## ENSG00000268903   0.000000    2.980969    0.000000    3.170474
## ENSG00000269981   0.000000    2.980969    0.000000    0.000000
## ENSG00000241860   3.355977    2.980969    0.000000    2.113649
## ENSG00000279457  58.170264   68.562276   63.718596   32.761567


根据基因在不同的样本中表达变化的差异程度mad值对数据排序，差异越大的基因排位越前。

normalized_counts_mad <- apply(normalized_counts, 1, mad)
normalized_counts <- normalized_counts[order(normalized_counts_mad, decreasing=T), ]

# 标准化后的数据输出
write.table(normalized_counts, file="ehbio_trans.Count_matrix.xls.DESeq2.normalized.xls",
quote=F, sep="\t", row.names=T, col.names=T)

# 只在Linux下能运行，目的是填补表格左上角内容
system(paste("sed -i '1 s/^/ID\t/'", "ehbio_trans.Count_matrix.xls.DESeq2.normalized.xls"))

# log转换后的结果
rld <- rlog(dds, blind=FALSE)
rlogMat <- assay(rld)
rlogMat <- rlogMat[order(normalized_counts_mad, decreasing=T), ]


write.table(rlogMat, file="ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.xls",
quote=F, sep="\t", row.names=T, col.names=T)

# 只在Linux下能运行，目的是填补表格左上角内容
system(paste("sed -i '1 s/^/ID\t/'", "ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.xls"))


样品层级聚类分析，判断样品的相似性和组间组内差异

# 生成颜色
hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)

# 计算相关性pearson correlation
pearson_cor <- as.matrix(cor(rlogMat, method="pearson"))

# 层级聚类
hc <- hcluster(t(rlogMat), method="pearson")

# 热图绘制
## 在命令行下运行时，需要去除下面开头的#号，把输出的图保存到文件中
## 输出到文件时，dev.off()命令是关闭输出，从而完成图片的写入。如果不做这一步，则图片则不能写入，从而不能打开
## 如果在Rstudio或其它可视化界面时，可以直接把图输出到屏幕
#pdf("ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.pearson.pdf", pointsize=10)
heatmap.2(pearson_cor, Rowv=as.dendrogram(hc), symm=T, trace="none",
col=hmcol, margins=c(11,11), main="The pearson correlation of each
sample")




#dev.off()


样品PCA分析

pca_data <- plotPCA(rld, intgroup=c("conditions"), returnData=T, ntop=5000)




差异基因分析

# 定义变量的好处是下面的代码就独立了
# 如果想比较不同的组，只需在这修改即可，后面代码都可以不动
sampleA = "untrt"
sampleB = "trt"


results函数提取差异基因分析结果，包含log2 fold changes, 

p values和adjusted p values. 

constrast用于指定比较的两组的信息，输出的log2FoldChange为log2(SampleA/SampleB)。

contrastV <- c("conditions", sampleA, sampleB)
res <- results(dds,  contrast=contrastV)
res

## log2 fold change (MLE): conditions untrt vs trt 
## Wald test p-value: conditions untrt vs trt 
## DataFrame with 26280 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE        stat    pvalue
##                  <numeric>      <numeric> <numeric>   <numeric> <numeric>
## ENSG00000227232 18.7141876     0.17695376 0.3809406  0.46451794 0.6422767
## ENSG00000278267  2.8144283     0.02802497 1.0249422  0.02734297 0.9781862
## ENSG00000268903  0.9807232    -1.71669654 2.3186888 -0.74037385 0.4590732
## ENSG00000269981  0.8256996     0.59077960 2.4538911  0.24075217 0.8097472
## ENSG00000241860  3.3713378     1.21773503 1.0789466  1.12863330 0.2590526
## ...                    ...            ...       ...         ...       ...
## 29256            0.3507226      1.8471273  3.484683  0.53007036 0.5960631
## 36308            0.9279831      0.1394489  1.598343  0.08724589 0.9304761
## 8414             0.7459738     -0.2678522  2.053434 -0.13044109 0.8962175
## 28837            1.2626343      0.2689443  1.405958  0.19128899 0.8482992
## 35320            1.0041030      0.5493562  1.717832  0.31979625 0.7491228
##                      padj
##                 <numeric>
## ENSG00000227232 0.8538515
## ENSG00000278267        NA
## ENSG00000268903        NA
## ENSG00000269981        NA
## ENSG00000241860        NA
## ...                   ...
## 29256                  NA
## 36308                  NA
## 8414                   NA
## 28837                  NA
## 35320                  NA


给DESeq2的原始输出结果增加样品平均表达信息，使得结果更容易理解和解析。

# 获得第一组数据均值
baseA <- counts(dds, normalized=TRUE)[, colData(dds)$conditions == sampleA]

if (is.vector(baseA)){
    baseMeanA <- as.data.frame(baseA)
} else {
    baseMeanA <- as.data.frame(rowMeans(baseA))
}
colnames(baseMeanA) <- sampleA
head(baseMeanA)

##                      untrt
## ENSG00000227232 19.8507436
## ENSG00000278267  2.7778822
## ENSG00000268903  0.4235857
## ENSG00000269981  0.9061570
## ENSG00000241860  4.6300269
## ENSG00000279457 59.0627429

# 获得第二组数据均值
baseB <- counts(dds, normalized=TRUE)[, colData(dds)$conditions == sampleB]
if (is.vector(baseB)){
        baseMeanB <- as.data.frame(baseB)
} else {
        baseMeanB <- as.data.frame(rowMeans(baseB))
}
colnames(baseMeanB) <- sampleB
head(baseMeanB)

##                        trt
## ENSG00000227232 17.5776316
## ENSG00000278267  2.8509744
## ENSG00000268903  1.5378607
## ENSG00000269981  0.7452421
## ENSG00000241860  2.1126487
## ENSG00000279457 55.8031756

# 结果组合
res <- cbind(baseMeanA, baseMeanB, as.data.frame(res))
head(res)

##                      untrt        trt   baseMean log2FoldChange     lfcSE
## ENSG00000227232 19.8507436 17.5776316 18.7141876     0.17695376 0.3809406
## ENSG00000278267  2.7778822  2.8509744  2.8144283     0.02802497 1.0249422
## ENSG00000268903  0.4235857  1.5378607  0.9807232    -1.71669654 2.3186888
## ENSG00000269981  0.9061570  0.7452421  0.8256996     0.59077960 2.4538911
## ENSG00000241860  4.6300269  2.1126487  3.3713378     1.21773503 1.0789466
## ENSG00000279457 59.0627429 55.8031756 57.4329592     0.08926886 0.2929300
##                        stat    pvalue      padj
## ENSG00000227232  0.46451794 0.6422767 0.8538515
## ENSG00000278267  0.02734297 0.9781862        NA
## ENSG00000268903 -0.74037385 0.4590732        NA
## ENSG00000269981  0.24075217 0.8097472        NA
## ENSG00000241860  1.12863330 0.2590526        NA
## ENSG00000279457  0.30474470 0.7605606 0.9130580

# 增加ID信息
res <- cbind(ID=rownames(res), as.data.frame(res))
res$baseMean <- rowMeans(cbind(baseA, baseB))

# 校正后p-value为NA的复制为1
res$padj[is.na(res$padj)] <- 1

# 按pvalue排序, 把差异大的基因放前面
res <- res[order(res$pvalue),]
head(res)

##                              ID     untrt       trt   baseMean
## ENSG00000152583 ENSG00000152583   77.8512  1900.412   989.1314
## ENSG00000148175 ENSG00000148175 7233.0010 19483.552 13358.2766
## ENSG00000179094 ENSG00000179094  151.7491  1380.881   766.3151
## ENSG00000134686 ENSG00000134686 1554.0390  4082.440  2818.2395
## ENSG00000125148 ENSG00000125148 1298.0392  5958.946  3628.4926
## ENSG00000120129 ENSG00000120129  773.9769  5975.522  3374.7494
##                 log2FoldChange      lfcSE      stat        pvalue
## ENSG00000152583      -4.608311 0.21090819 -21.84984 7.800059e-106
## ENSG00000148175      -1.429585 0.08639329 -16.54741  1.671380e-61
## ENSG00000179094      -3.184674 0.20065896 -15.87108  1.005030e-56
## ENSG00000134686      -1.392749 0.09190323 -15.15451  7.073605e-52
## ENSG00000125148      -2.199191 0.14796048 -14.86337  5.698869e-50
## ENSG00000120129      -2.948122 0.19931242 -14.79146  1.663007e-49
##                          padj
## ENSG00000152583 1.331002e-101
## ENSG00000148175  1.426021e-57
## ENSG00000179094  5.716612e-53
## ENSG00000134686  3.017600e-48
## ENSG00000125148  1.944910e-46
## ENSG00000120129  4.729591e-46


整体分析结果输出到文件

comp314 <- paste(sampleA, "_vs_", sampleB, sep=".")

# 生成文件名
file_base <- paste("ehbio_trans.Count_matrix.xls.DESeq2", comp314, sep=".")
file_base1 <- paste(file_base, "results.xls", sep=".")
write.table(as.data.frame(res), file=file_base1, sep="\t", quote=F, row.names=F)


提取差异表达基因

# 差异基因筛选，padj<0.1
res_de <- subset(res, res$padj<0.1, select=c('ID', sampleA,
        sampleB, 'log2FoldChange', 'padj'))
# foldchang > 1
res_de_up <- subset(res_de, res_de$log2FoldChange>=1)

file <- paste("ehbio_trans.Count_matrix.xls.DESeq2",sampleA,"_higherThan_",sampleB,
        'xls', sep=".") 
write.table(as.data.frame(res_de_up), file=file, sep="\t", quote=F, row.names=F)

res_de_dw <- subset(res_de, res_de$log2FoldChange<=(-1)*1)

file <- paste("ehbio_trans.Count_matrix.xls.DESeq2",sampleA, "_lowerThan_", sampleB, 
        'xls', sep=".") 
write.table(as.data.frame(res_de_dw), file=file, sep="\t", quote=F, row.names=F)

# 差异基因ID
res_de_up_id = data.frame(ID=res_de_up$ID, 
        type=paste(sampleA,"_higherThan_", sampleB, sep="."))
res_de_dw_id = data.frame(ID=res_de_dw$ID, 
        type=paste(sampleA,"_lowerThan_", sampleB, sep="."))
de_id = rbind(res_de_up_id, res_de_dw_id)

file <- "ehbio_trans.Count_matrix.xls.DESeq2.all.DE"
write.table(as.data.frame(de_id), file=file, sep="\t", quote=F, row.names=F, col.names=F)


绘制火山图

# 可以把前面生成的results.xls文件提交到www.ehbio.com/ImageGP绘制火山图
# 或者使用我们的s-plot https://github.com/Tong-Chen/s-plot
logCounts <- log2(res$baseMean+1)
logFC <- res$log2FoldChange
FDR <- res$padj
#png(filename=paste(file_base, "Volcano.png", sep="."))
plot(logFC, -1*log10(FDR), col=ifelse(FDR<=0.01, "red", "black"),
 xlab="logFC", ylab="-1*log1o(FDR)", main="Volcano plot", pch=".")
#dev.off()




差异基因热图

res_de_up_top20_id <- as.vector(head(res_de_up$ID,20))
res_de_dw_top20_id <- as.vector(head(res_de_dw$ID,20))

red_de_top20 <- c(res_de_up_top20_id, res_de_dw_top20_id)
red_de_top20

##  [1] "ENSG00000178695" "ENSG00000196517" "ENSG00000116584"
##  [4] "ENSG00000144369" "ENSG00000276600" "ENSG00000108821"
##  [7] "ENSG00000162692" "ENSG00000181467" "ENSG00000145777"
## [10] "ENSG00000163491" "ENSG00000183160" "ENSG00000172986"
## [13] "ENSG00000164484" "ENSG00000131389" "ENSG00000134253"
## [16] "ENSG00000124766" "ENSG00000227051" "ENSG00000146006"
## [19] "ENSG00000112837" "ENSG00000146592" "ENSG00000152583"
## [22] "ENSG00000148175" "ENSG00000179094" "ENSG00000134686"
## [25] "ENSG00000125148" "ENSG00000120129" "ENSG00000189221"
## [28] "ENSG00000109906" "ENSG00000101347" "ENSG00000162614"
## [31] "ENSG00000096060" "ENSG00000162616" "ENSG00000166741"
## [34] "ENSG00000183044" "ENSG00000164125" "ENSG00000198624"
## [37] "ENSG00000256235" "ENSG00000077684" "ENSG00000135821"
## [40] "ENSG00000164647"

# 可以把矩阵存储，提交到www.ehbio.com/ImageGP绘制火山图
# 或者使用我们的s-plot https://github.com/Tong-Chen/s-plot
red_de_top20_expr <- normalized_counts[rownames(normalized_counts) %in% red_de_top20,]
library(pheatmap)
pheatmap(red_de_top20_expr, cluster_row=T, scale="row", annotation_col=sample)






关于批次效应

见论坛讨论 http://www.ehbio.com/Esx/d/10-deseq2。

每个DESeqDataSet对象都要有一个实验设计formula，用于对数据进行分组，以便计算表达值的离散度和估计表达倍数差异，通常格式为~ batch + conditions 

(为了方便后续计算，最为关注的分组信息放在最后一位)。

如果记录了样本的批次信息，或者其它需要抹除的信息可以定义在design参数中，在下游回归的分析中会根据design formula来估计batch effect的影响，并在下游分析时减去这个影响。这是处理batch effect的推荐方式。

在模型中考虑batch effect并没有在数据矩阵中移除bacth effect，如果下游处理时，确实有需要可以使用limma包的removeBatchEffect来处理。

countData: 表达矩阵

colData: 样品分组信息表

design: 实验设计信息，batch和conditions必须是colData中的一列

ddsFullCountTable <- DESeqDataSetFromMatrix(countData = data,
        colData = sample,  design= ~ batch + conditions)

dds <- DESeq(ddsFullCountTable)

rld <- rlog(dds, blind=FALSE)
rlogMat <- assay(rld)
rlogMat <- limma::removeBatchEffect(rlogMat, c(sample$batch))


Just to be clear, there’s an important difference between removing a 

batch effect and modelling a batch effect. Including the batch in your 

design formula will model the batch effect in the regression step, which 

means that the raw data are not modified (so the batch effect is not 

removed), but instead the regression will estimate the size of the batch 

effect and subtract it out when performing all other tests. In addition, 

the model’s residual degrees of freedom will be reduced appropriately to 

reflect the fact that some degrees of freedom were “spent” modelling the 

batch effects. This is the preferred approach for any method that is 

capable of using it (this includes DESeq2). You would only remove the 

batch effect (e.g. using limma’s removeBatchEffect function) if you were 

going to do some kind of downstream analysis that can’t model the batch 

effects, such as training a classifier.

批次效应模拟

#Make some simulated data with a batch effect:
dds <- makeExampleDESeqDataSet(betaSD=1,interceptMean=10)
dds$batch <- factor(rep(c("A","B"),each=6))

#VST, remove batch effect, then plotPCA:

vsd <- vst(dds)
plotPCA(vsd, "batch")




assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$batch)
plotPCA(vsd, "batch")




SVA(批次未记录时，寻找潜在影响因子，并矫正)

dat <- counts(dds, normalized=TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx,]
mod <- model.matrix(~ dex, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
##calculating the variables
n.sv <- num.sv(dat,mod,method="leek") #gives 11
### using 4
lnj.corr <- svaBatchCor(dat,mod,mod0,n.sv=4)
co <- lnj.corr$corrected
plPCA(co)
#http://www.bioconductor.org/packages/release/bioc/vignettes/sva/inst/doc/sva.pdf
#http://bioconductor.org/help/workflows/rnaseqGene/#removing-hidden-batch-effects




参考资料


https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
http://dx.doi.org/10.1186/s13059-014-0550-8

数据准备

接上一篇《miRNA seq差异表达分析练习（一）——GEO样本数据下载》，下载的数据保存在文件夹miRNA_seq中

有60个txt.gz文件



从后缀名可以看出，这些都是压缩文件，直接打开是乱码的，解压后打开就可以看到具体内容了



事实上，用R语言是可以直接读取txt.gz文件而不需要解压的。

接下来，我将用R语言批量读取该文件夹下所有的txt.gz文件，并合并成一个包含60个样本数据的数据框。

# 读取D:\用户\桌面\\miRNA_seq文件夹下所有文件
setwd('D:\\用户\\桌面\\miRNA_seq')
fileName <- dir()  # 获取所有文本名

# 先读取第一个文件
mydata <- read.table(fileName[1], header = T, row.names = 1)
names(mydata) <- c('sample_1')

# 读取其余文件，并逐一并到第一个文件上
for(i in 2:length(fileName)){
  data <- read.table(fileName[i],header = T, row.names = 1)
  names(data) <- c(paste('sample',as.character(i),sep = '_'))
  mydata <- cbind(mydata,data)
}
names(mydata) <- c(paste('N-',1:20,sep = ''),paste('T-',1:20,sep = ''),paste('P-',1:20,sep = ''))

# 对所有数据四舍五入取整
mydata <- round(mydata,0)

# 查看mydata数据框前6行
head(mydata)

合并后的数据框mydata，行名即miRNA名称，列名为各样本名称。

因为DESeq包要求的表达矩阵数据为整数，所以在正式分析前，应将所有数据四舍五入取整数，一行代码搞定！

 

# 对所有数据四舍五入取整
mydata <- round(mydata,0)  # 0表示保留0给小数，即整数

DESeq差异表达分析

差异表达分析需要两个数据：表达矩阵和分组信息

表达矩阵简单获得，只需将前面的数据框转化为矩阵即可。

mymatrix <- as.matrix(mydata)

分组信息也很简单，60个样本，三个组织即分为3组，每组20个,其中,N组是正常组织,T和P组都属于癌变组织。

condition <- factor(c(rep('normal',20),rep('cancerous',40)))
org <- factor(c(rep('N',20),rep('T',20),rep('P',20)))
coldata <- data.frame(row.names = colnames(mydata),org)

OK，万事具备，差异分析走起！

# 差异表达分析
library(DESeq2)
dds <- DESeqDataSetFromMatrix(mymatrix,DataFrame(coldata), design = ~org)
# 标准化
dds = DESeq(dds)
# 提取结果
res <- results(dds,alpha = 0.001)
res = res[order(res$padj),]
summary(res) #查看结果
# 输出所有结果
# write.csv(res,file="D:\\用户\\桌面\\All_results.csv")

搞定！差异表达分析似乎也没有想象中那么难哈，下一篇进行结果分析与保存。

啰嗦一句，DESeq2包与其他包的安装有些不同，不过方法已经给你备好了！

source("http://bioconductor.org/biocLite.R")

biocLite("DESeq2")

安装时间还是挺久的，不过，磨刀不误砍柴工，一次安装，永久使用！

 

 


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




三者区别

 



三者区别

差异分析流程:

1 初始数据

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

为什么要标准化？

消除文库大小不同，测序深度对差异分析结果的影响

怎样标准化？

找到一个能反映文库大小的因子，利用这个因子对rawdata进行标准化

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

4 多重假设得FDR值：

检验方法：Wald、LRT

多重检验：BH

5 差异基因筛选：pvalue、padj

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


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


0. 前期工作

用到的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)