跟着存档教程动手学RNAseq分析（一）

原始资料：https://hbctraining.github.io/DGE_workshop/

这个中文教程是原始资料的关键内容过一遍，以帮助自己的研究和学习。如果对R不熟悉，推荐学习 Introduction to R[1]。

学习目标

• 质控
• 利用DESeq2进行差异分析
• 可视化差异基因表达模式
• 功能分析

准备工作

pkgs = c("BiocManager", "RColorBrewer", "pheatmap", "ggrepel", "devtools", "tidyverse")
sapply(pkgs, function(p) install.packages(p))

# Bioconductor package
pkgs = c("DESeq2",
         "clusterProfiler",
         "DOSE",
         "org.Hs.eg.db",
         "pathview",
         "DEGreport",
         "EnsDb.Hsapiens.v86",
         "AnnotationHub",
         "ensembldb")
BiocManager::install(pkgs)

加载这些包，确定没有问题：

library(DESeq2)
library(ggplot2)
library(RColorBrewer)
library(pheatmap)
library(ggrepel)
library(clusterProfiler)
library(DEGreport)
library(org.Hs.eg.db)
library(DOSE)
library(pathview)
library(tidyverse)
library(EnsDb.Hsapiens.v86)
library(AnnotationHub)
library(ensembldb)

没有问题就Ok了，再正式学习之前感谢下教程的创作团队。

Harvard Chan Bioinformatics Core

img

基因水平差异表达分析概述

相关数据文件 📎DGE_workshop_salmon.zip[2]

概览

RNAseq差异分析的目的是评估不同实验条件下哪些基因/功能发生了改变。一个典型的RNAseq分析流程如下图所示：

img

在接下来的几节内容中，我们将带你通过使用各种R包完成端到端基因水平RNA-seq差异表达工作流程。我们将从读取Salmon获得的数据开始，将伪计数转换为计数，进行探索性数据分析以进行质量评估，并探索样本之间的关系，进行差异表达分析，并在进行下游功能分析之前可视化地研究结果。

数据介绍

数据来源于 Kenny PJ et al, Cell Rep 2014[3] 的一部分内容。

RNA-Seq在HEK293F细胞上进行，这些细胞要么转染了MOV10转基因，要么转染了siRNA以抑制MOV10的表达，要么转染了非特异性（无关的）siRNA。这导致了三个条件：Mov10 oe（过表达），Mov10 kd（敲除）和Irrelevant kd，分别。重复数量如下所示。

img

利用这些数据，我们将评估与MOV10表达干扰相关的转录模式。请注意，无关的siRNA组别将作为我们的控制条件（对照组）。

What is the purpose of these datasets? What does Mov10 do?

The authors are investigating interactions between various genes involved in Fragile X syndrome, a disease in which there is aberrant production of the FMRP protein.

FMRP is “most commonly found in the brain, is essential for normal cognitive development and female reproductive function. Mutations of this gene can lead to fragile X syndrome, mental retardation, premature ovarian failure, autism, Parkinson’s disease, developmental delays and other cognitive deficits.” - from wikipedia[4]

MOV10, is a putative RNA helicase that is also associated with FMRP in the context of the microRNA pathway.

The hypothesis the paper[5] is testing is that FMRP and MOV10 associate and regulate the translation of a subset of RNAs.

img

我们在此处探索的问题是：

1. MOV10过表达和缺失的表达模式是什么
2. 两种情况下共有的基因有哪些

设置

1. 下载Salmon处理后的结果文件。📎data.zip[6]
2. 下载注释文件，用于基因名称的映射。你可以通过 https://github.com/hbctraining/DGE_workshop_salmon/raw/master/data/tx2gene_grch38_ens94.txt 保存，也可以在仓库目录中找到。

img

将下载的Salmon文件夹内的数据文件放到data目录下：

img

我们通过下载的教程文件夹创建一个RStudio项目。

img

img

在控制台输入下面命令创建一个用于差异分析的脚本：

file.edit("de_script.R")

我们就在这个脚本中键入代码和运行查看结果。

载入包

## Gene-level differential expression analysis using DESeq2

## Setup
### Bioconductor and CRAN libraries used
library(DESeq2)
library(tidyverse)
library(RColorBrewer)
library(pheatmap)
library(DEGreport)
library(tximport)
library(ggplot2)
library(ggrepel)

载入数据

Salmon的主要输出是一个 quant.sf 文件。我们看下输出的格式：

img

NOTE: The effective gene length in a sample is the average of the transcript lengths after weighting for their relative expression. You may see effective lengths that are larger than the physical length. The interpretation would be that in this case, given the sequence composition of these transcripts (including both the sequence-specific and fragment GC biases), you’d expect a priori to sample more reads from them — thus they have a longer estimated effective length.

Salmon生成的伪计数表示为标准化TPM计数(transcripts per million)，并映射到转录本。为了执行DESeq2分析，这些需要转换为非标准化计数估计。为了使用DESeq2，我们还需要将我们的丰度估计从转录水平分解到基因水平。我们将使用R Bioconductor包tximport来完成上述所有操作，并为DESeq2进行设置。

延伸阅读：

RPKM FPKM TPM[7] | https://www.cnblogs.com/Belter/p/13205635.html

The confusion of using TPM (transcripts per million)[8]

TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-seq Data from the NCI Patient-Derived Models Repository[9]

（能用标准的normalized count尽量使用，TPM, FPKM都不适用于组间比较）

列出数据文件，解析样本名：

## List all directories containing data
samples <- list.files(path = "./data", full.names = T, pattern="salmon$")

## Obtain a vector of all filenames including the path
files <- file.path(samples, "quant.sf")

## Since all quant files have the same name it is useful to have names for each element
names(files) <- str_replace(samples, "./data/", "") %>%
  str_replace(".salmon", "")

我们的Salmon索引是由Ensembl id列出的转录序列生成的，但tximport需要知道这些转录本来自哪些基因。我们将使用下载的注释表提取转录本的基因信息。

# Load the annotation table for GrCh38
tx2gene <- read.delim("data/tx2gene_grch38_ens94.txt")

# Take a look at it
tx2gene %>% View()

tx2gene是一个三列的数据帧，将转录ID(第1列)连接到基因ID(第2列)到基因符号(第3列)。我们将把前两列作为输入到tximport。列名是不相关的，但列的顺序是(即转录ID必须是第一)。

现在，我们已经准备好运行tximport了。请注意，虽然在我们的quant.sf文件中有一列对应于每个文本的估计计数值，但这些值与有效长度相关。我们要做的是使用countsFromAbundance= "lengthScaledTPM"参数。这将使用TPM列，并计算与原始计数相同规模的数量，只是不再与跨样本的transcript长度相关。

?tximport   # let's take a look at the arguments for the tximport function

# Run tximport
txi <- tximport(files, type="salmon", tx2gene=tx2gene[,c("tx_id", "ensgene")], countsFromAbundance="lengthScaledTPM")

如果你自己分析使用的基因ID数据包含后缀，如ENSG00000265439.2，可以使用ignoreTxVersion = TRUE忽略版本号。

查看数据

txt对象是一个列表，包含以下的一些属性：

> attributes(txi)
$names
[1] "abundance"           "counts"              "length"              "countsFromAbundance"

让我们看下计数矩阵，你可能注意到有小数点，让我们将其取整到距离最近的数值，并将计数矩阵转换为数据框：

# Look at the counts
txi$counts %>% View()

# Write the counts to an object
data <- txi$counts %>% 
  round() %>% 
  data.frame()

创建元数据

跟踪我们的数据信息是非常重要的。至少，我们需要有一个文件将我们的样本映射到我们正在调查的相应样本组。我们将使用来自计数矩阵的列名作为元数据文件的行名，并使用一列来标识每个示例为“MOV10_overexpression”、“MOV10_knockdown”或“control”。

## Create a sampletable/metadata
sampletype <- factor(c(rep("control",3), rep("MOV10_knockdown", 2), rep("MOV10_overexpression", 3)))
meta <- data.frame(sampletype, row.names = colnames(txi$counts))

差异基因表达分析综述

那么这个计数数据到底代表什么呢？用于差异表达分析的计数数据表示来自特定基因的序列读数。计数越多，与该基因相关的读数就越多，这就意味着样本中该基因的表达水平较高。

img

通过差异表达分析，我们寻找两个或多个组（在元数据中定义）之间表达变化的基因。

• case 对比 control
• 表达与某些变量或临床结果的相关性

为什么不能通过根据基因在两组之间的差异程度(基于倍数变化值)来对基因进行排序来识别差异表达基因?

img

通常情况下，你的数据比你预期的要复杂得多。不同样本之间表达水平不同的基因不仅是感兴趣的实验变量的结果，也是外来因素的结果。差异表达分析的目标是确定这些效应的相对作用，并将“有趣”与“无趣”变量区分开来。

img

即使样本组之间的平均表达水平看起来有很大的不同，也有可能这种差异实际上并不显著。下图说明了“未处理”和“处理”组之间的“GeneA”表达。“治疗”组的平均基因A表达水平是“未治疗”组的两倍。但是，考虑到组内（重复）观察到的变异，组间的表达(计数)差异是否显著?

在确定基因是否有差异表达时，我们需要考虑数据中的差异（以及它可能来自哪里）。

img

RNA-seq计数分布

为了检验显著性，我们需要一个合适的统计模型，该模型可以准确地进行归一化(考虑到测序深度的差异等)和方差建模(考虑到重复数量少和动态表达范围大)。

为了确定适当的统计模型，我们需要有关计数分布的信息。为了了解RNA-seq计数是如何分布的，让我们绘制单个样本‘Mov10_oe_1’的计数图：

ggplot(data) +
  geom_histogram(aes(x = Mov10_oe_1), stat = "bin", bins = 200) +
  xlab("Raw expression counts") +
  ylab("Number of genes")

img

这个图展示了RNA-seq计数数据的通用特征：

• 低数量的计数与大比例的基因相关
• 由于没有表达上限而产生一个右长尾
• 数值变化范围很大

微阵列数据的对数强度近似于正态分布。然而，由于RNA-seq计数数据的不同性质，如整数计数而不是连续测量和非正态分布的数据，正态分布不能准确地建模RNA-seq计数（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3541212/）。

计数数据建模

一般的计数数据可以用不同的分布来建模：

1. 二项分布：给出投掷一枚硬币多次得到若干正面的概率。基于离散事件并用于特定数量的情况下。
2. 泊松分布(Poisson distribution):当案例数量非常大(即买彩票的人)，但事件发生的概率非常小(中奖的概率)时使用。泊松与二项式相似，但它是基于连续事件的。适用于均值==方差的数据。
3. 负二项（Negative binomial）分布:泊松的近似，但有一个额外的参数，调整方差独立于均值。

img

Details provided by Rafael Irizarry in the EdX class.[10]

那么我们用什么来处理RNA测序计数数据呢?

RNA-Seq数据中有非常多数目的RNA，提取到特定转录本的概率非常小。因此，使用泊松分布或负二项分布是一种合适的情况。选择一个而不是另一个将取决于我们数据中的平均值和方差之间的关系。

运行以下代码绘制' Mov10过表达'重复样本的均值与方差图：

mean_counts <- apply(data[,6:8], 1, mean)        #The second argument '1' of 'apply' function indicates the function being applied to rows. Use '2' if applied to columns 
variance_counts <- apply(data[,6:8], 1, var)
df <- data.frame(mean_counts, variance_counts)

ggplot(df) +
        geom_point(aes(x=mean_counts, y=variance_counts)) + 
        scale_y_log10(limits = c(1,1e9)) +
        scale_x_log10(limits = c(1,1e9)) +
        geom_abline(intercept = 0, slope = 1, color="red")

img

上面图中的红线代表Y=X。这里有一些重点：

1. 重复间的变异倾向于大于平均值(红线)，特别是对于平均表达水平高的基因。
2. 对于低表达的基因，我们看到相当多的分散。我们通常称之为“异方差性”。也就是说，对于一个给定的表达水平，我们在方差的数量上观察到很多变化。

这很好地说明了我们的数据不符合泊松分布。如果mRNA的比例在一个样本组的生物复制之间完全保持恒定，我们可以期望泊松分布(其中均值==方差)。或者，如果我们继续添加更多的重复样本(即> 20)，我们最终会看到散点开始减少，高表达数据点向红线靠近。所以理论上，如果我们有足够的重复样本，我们可以使用泊松分布。

然而，在实践中，大量的重复样本要么很难获得(取决于如何获得样本)，要么负担不起。更常见的情况是，数据集只有少量重复(~3-5)，并且它们之间存在适量的变异。考虑到这种重复间的变异类型，最适合的模型是负二项(NB)模型。本质上，NB模型是均值<方差的数据的一个很好的近似值，就像RNA-Seq计数数据一样。

注意：如果我们使用泊松，这将低估可变性，导致假阳性差异表达基因的增加。

利用生物重复改进平均估计

额外重复的价值在于，当你添加更多的数据时，你会得到越来越精确的群体平均估计值，并最终在区分样本类别差异的能力(即更多的DE基因)上有更大的信心。

• 生物学上的重复代表代表同一类样本的多个样本(即来自不同小鼠的RNA)。
• 技术重复代表相同的样本(即来自同一只老鼠的RNA)，但技术步骤被重复。
• 通常生物差异比技术差异大得多，所以我们不需要考虑技术差异来确定表达中的生物差异。
• 不要在技术重复上花钱——生物重复更有用。

注意:如果你正在使用细胞株，并且不确定你是否准备了生物或技术重复，请查看这个链接[11]。这是一个有用的资源，可以帮助你确定如何最好地设置你的体外实验。

下图显示了测序深度与识别出的差异表达基因数量上的重复样本数目之间的关系。

img

注意，与增加测序深度相比，重复数量的增加往往会返回更多的DE基因。因此，一般情况下，重复越多越好，但需要注意的是，低表达DE基因的检测和isoform水平的差异表达需要更高的深度。

差分表达分析工作流程

为了在进行差异表达分析时适当地建模计数，已经开发了许多用于RNA-seq数据的差异表达分析的软件包。即使新方法不断被开发，一些工具通常被推荐为最佳实践，例如**DESeq2**[12]和**EdgeR**[13]。这两种工具都使用负二项模型，使用类似的方法，并且通常产生类似的结果。它们非常严格，在敏感性和特异性之间有很好的平衡(减少假阳性和假阴性)。

Limma-Voom[14]是另一套经常用于DE分析的工具，但这种方法可能对小样品量不太敏感。当每组生物重复数增加较大(> ~ 20)时，推荐采用该方法。

许多描述这些方法之间比较的研究表明，虽然有一些一致，但不同的工具之间也有很大的差异。此外，没有一种方法可以在所有条件下都表现得最优(Soneson和Dleorenzi, 2013[15])。

我们将使用DESeq2进行DE分析，使用DESeq2的分析步骤在下面的流程图中以绿色显示。DESeq2首先将计数数据归一化，以解释样品之间文库大小和RNA组成的差异。然后，我们将使用标准化计数在基因和样本水平上为QC绘制一些图。最后一步是使用来自DESeq2包的适当函数来执行差异表达式分析。我们将在接下来的课程中深入研究这些步骤，但是关于DESeq2的更多细节和有用的建议可以在DESeq2 vignette[16]找到。

img

参考资料

[1]

Introduction to R: https://hbctraining.github.io/Intro-to-R/

[2]

📎DGE_workshop_salmon.zip: https://www.yuque.com/attachments/yuque/0/2022/zip/471931/1649675646093-6bd74a9a-b78c-4a55-af55-7aeff6dc8336.zip

[3]

Kenny PJ et al, Cell Rep 2014: http://www.ncbi.nlm.nih.gov/pubmed/25464849

[4]

wikipedia: https://en.wikipedia.org/wiki/FMR1

[5]

the paper: http://www.ncbi.nlm.nih.gov/pubmed/25464849

[6]

📎data.zip: https://www.yuque.com/attachments/yuque/0/2022/zip/471931/1649676597231-13daedf2-1273-405d-a29f-27ad7e9dd0d2.zip

[7]

RPKM FPKM TPM: https://www.jianshu.com/p/2cce4376be48

[8]

The confusion of using TPM (transcripts per million): https://bioinformatics.stackexchange.com/questions/16372/the-confusion-of-using-tpm-transcripts-per-million

[9]

TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-seq Data from the NCI Patient-Derived Models Repository: https://translational-medicine.biomedcentral.com/articles/10.1186/s12967-021-02936-w

[10]

Details provided by Rafael Irizarry in the EdX class.: https://youtu.be/fxtB8c3u6l8

[11]

链接: https://web.archive.org/web/20170807192514/http://www.labstats.net:80/articles/cell_culture_n.html

[12]

DESeq2: https://bioconductor.org/packages/release/bioc/html/DESeq2.html

[13]

EdgeR: https://bioconductor.org/packages/release/bioc/html/edgeR.html

[14]

Limma-Voom: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-2-r29

[15]

Soneson和Dleorenzi, 2013: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-14-91

[16]

DESeq2 vignette: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html