RNA-Seq 分析流程：多时间点样本分析实战（一）
小鱼爸
2023-07-12
简介
了解基因表达的时序动态变化是生物学的一个基本问题。此教程提供了多时间点数据的分步实战流程：(1)数据集的质量控制和标准化；(2)进行差异表达分析；(3)时序数据的聚类；(4)用GO term和KEGG通路富集分析解释聚类簇。作为实战流程，我们应用的数据为暴露于四种流感病毒株的小鼠的多时间点转录组数据。
本分析流程基于先前建立的用于发育时间序列的新策略。它依赖于几个用于分析基因表达数据的标准包，其中一些特定用于多时间点数据。各步骤使用的主要函数在moaninR包中。
安装和设置包
moanin和timecoursedata可从bioconductor获得，并且可以使用BiocManager包中的install功能进行安装：
library (BiocManager)
BiocManager::install("timecoursedata")
#如果下载失败，可以去Bioconductor官网（http://www.bioconductor.org/packages/release/data/experiment/html/timecoursedata.html）下载到本地，然后安装。
BiocManager::install("moanin")
此工作流程需要以下附加包：
# From Github
library(moanin)
library(timecoursedata)
# From CRAN
library(NMF)
library(ggfortify)
# From Bioconductor
library(topGO)
library(biomaRt)
#library(KEGGprofile)
library(BiocWorkflowTools) #Needed for compiling the .Rmd script
library(pander) #Needed for compiling the .Rmd script
如果安装了Bioconductor，则可以通过以下方式安装上面的CRAN和Bioconductor软件包：
BiocManager::install(
c("NMF", "ggfortify", "topGO", "biomaRt",
"BiocWorkflowTools", "timecoursedata", ))
小鼠肺组织对流感的动态反应分析
该分析流程使用微阵列时序实验的数据进行示范，将小鼠暴露于三种不同的流感毒株，并在感染后14个时间点（0、3、6、9、12、18、24、30、36、48、60 小时， 3、5 和 7 天后）收集肺组织。研究中使用的三种流感毒株是(1)低致病性季节性H1N1流感病毒(A/Kawasaki/UTK4/2009 [H1N1])；(2)2009年大流行季节的一种轻度致病性病毒(A/California/04/2009 [H1N1])；(3)高致病性H5N1禽流感病毒(A/Vietnam/1203/2004[H5N1])。给小鼠注射了每种病毒105 PFU。(4)另外，42只小鼠注射了较低剂量的越南禽流感病毒(103 PFU)。
通过将基因表达时序数据与病毒生长数据相结合，作者发现，肺组织的炎症反应是门控的（gated），肺部的病毒浓度存在阈值（threshold）。一旦超过这个阈值，就会产生强烈的炎症和细胞因子。这些结果证明病理反应受病毒浓度的非线性调节。 此外，Varoquaux等人利用类似的步骤来分析作物S. bicolor对干旱的终生转录组响应的RNA-seq数据。
数据概览
首先让我们加载数据。moanin包包含(Shoemaker et al., 2015)的标准化数据和元数据。
# Now load in the metadata
data(shoemaker2015)
meta = shoemaker2015$meta
data = shoemaker2015$data
元数据meta包含有关组、重复和每次观察的时间点的信息。
head(meta)
| Group | Replicate | Timepoint | ||
|---|---|---|---|---|
| GSM1557140 | 0 | K | 1 | 0 |
| GSM1557141 | 1 | K | 2 | 0 |
| GSM1557142 | 2 | K | 3 | 0 |
| GSM1557143 | 3 | K | 1 | 12 |
| GSM1557144 | 4 | K | 2 | 12 |
| GSM1557145 | 5 | K | 3 | 12 |
在深入探索性分析和质量控制之前，我们先分配整个分析中使用配色方案。我们将组和时间点的配色方案定义为向量vectors。我们还定义了一系列标记（或绘图符号）来区分散点图中的重复样本。我们还对描述处理的组进行了重新排序，以便按致病性从低到高排序（对照在前）。
group_colors = c(
"M"="dodgerblue4",
"K"="gold",
"C"="orange",
"VH"="red4",
"VL"="red2")
time_colors = grDevices::rainbow(15) [1:14]
names(time_colors) = c(0, 3, 6, 9,12, 18, 24, 30, 36, 48, 60, 72, 120, 168)
# Combine all color schemes into one named lists.
ann_colors = list(
Timepoint=time_colors,
Group=group_colors
)
replicate_markers = c(15, 17, 19)
names(replicate_markers) = c(1, 2, 3)
ann_markers = list(
Replicate=replicate_markers)
# Reorder the conditions such that:
# - Control is before any influenza treatment
# - Each treatment is ordered from low to high pathogeny
meta$Group = factor(meta$Group, levels(meta$Group)[c(3, 2, 1, 5, 4)])
质量控制和标准化
基因表达数据分析的第一步始终是对数据进行标准化和质量控制检查。通常，在归一化之前和之后执行两个质量控制和探索性分析步骤：(1)样本的降维分析；(2)各样本之间的相关性分析。在这两种情况下，我们期望产生高表达的生物信号，而重复样本应该强烈聚集或彼此相关。 在进行下游分析之前，我们只保留了高度可变的基因：在这一步，我们只保留前 50% 最可变的基因。在图2中，我们绘制了所有基因的方差分布。
variance_cutoff = 0.5
# Filter genes by median absolute deviation (mad)
variance_per_genes = apply(data, 1, mad)
min_variance = quantile(variance_per_genes, c(variance_cutoff))
variance_filtered_data = data[variance_per_genes > min_variance,]
hist(variance_per_genes, breaks=100, col="black", xlab="Variance",
ylab="# probes", main="")
abline(v=min_variance, col="#AB0000", lw=2)
让我们首先进行主成分分析（PCA）分析。
pca_data = prcomp(t(variance_filtered_data), rank=3, center=TRUE, scale=TRUE)
percent_var = round(100 * attr(pca_data, "percentVar"))
在图3中，我们通过绘制主成分(PC)来可视化数据，样本按其条件或采样时间着色，并且每个样本复制用不同的符号区分。我们可以看到后面的时间点有很大的差异。
par(mfrow=c(2, 2))
par(mar=c(4.5, 4.5, 1.5, 1.5))
plot(
pca_data$x[, c("PC1","PC2")],
col=ann_colors$Group[meta$Group],
pch=ann_markers$Replicate[as.factor(meta$Replicate)], bty="l")
mtext(at=275, text="Colored by Group", line=0.5, font=2, adj=0.5)
legend(x=250, y=25, legend=names(ann_colors$Group),
fill=ann_colors$Group, bty="n", xpd=NA, yjust=0.5)
par(mar=c(4.5, 5.5, 1.5, .1))
plot(
pca_data$x[, c("PC3","PC2")],
col=ann_colors$Group[meta$Group],
pch=ann_markers$Replicate[as.factor(meta$Replicate)],
ylab="", bty="l", yaxt="n")
axis(2, labels=FALSE)
par(mar=c(4.5, 4.5, 1.5, 1.5))
plot(
pca_data$x[, c("PC1","PC2")],
col=ann_colors$Timepoint[as.factor(meta$Timepoint)],
pch=ann_markers$Replicate[as.factor(meta$Replicate)], bty="l")
mtext(at=275, text="Colored by Timepoint", line=0.5, font=2, adj=0.5)
legend(x=250, y=-15, legend=names(ann_colors$Timepoint),
fill=ann_colors$Timepoint, bty="n", xpd=NA, yjust=0.5)
par(mar=c(4.5, 5.5, 1.5, .1))
plot(
pca_data$x[, c("PC3","PC2")],
col=ann_colors$Timepoint[as.factor(meta$Timepoint)],
pch=ann_markers$Replicate[as.factor(meta$Replicate)],
ylab="", bty="l", yaxt="n")
axis(2, labels=FALSE)
我们还在图4中将每个样本之间的Pearson相关性绘制为热图（使用NMF中的函数aheatmap）。我们按组（处理）和时间点（采样时间）对样本进行排序。
# Reorder genes on condition, time, and replicate
ord = order(
meta$Group,
meta$Timepoint,
meta$Replicate)
variance_filtered_data = variance_filtered_data[, ord]
data_corr = cor(variance_filtered_data, method="pearson")
data_corr_meta = meta[ord, ]
data_corr_meta$Timepoint = as.factor(data_corr_meta$Timepoint)
# use aheatmap from the NMF package for a heatmap
aheatmap(
data_corr,
Colv=NA, Rowv=NA,
annCol=data_corr_meta[, c("Group", "Timepoint")],
annRow=data_corr_meta[, c("Group", "Timepoint")],
annLegend=TRUE,
annColors=ann_colors,
main="Correlation plot")
我们已经可以从相关图中看到有趣的结果。首先，从对照小鼠采集的样本之间的互相关性高于其余治疗组之间的互相关性。其次，感染流感的小鼠直到时间点36才会有轻微反应。第三，菌株的致病性越低，样品越接近对照条件。第四，时间点120和168的越南样本是与对照样本差异最大的样本。
时序数据的差异表达 (DE) 分析
时序数据的 DE 分析方法
基因表达分析的下一步通常是运行差异表达分析，通常是为了找到不同条件下的差异基因。 对于时序数据，有两种不同的方法来确定差异表达基因，
1）按时间点分析，我们将每个时间点视为不同的条件，并确定哪些基因在特定时间点之间或单个时间点的特定条件之间发生变化。
2）全局分析，我们考察随着时间的推移全局的表达模式，并推断哪些基因在不同条件下具有不同的模式，或者随着时间的推移具有变化的模式。该方法的第一步是为每个基因拟合一个样条模型(Spline model, Storey et al., 2005)，然后使用该样条模型来测试随时间变化的不同类型的差异表达。
按时间点分析使用经典的差异表达方法，并且通常是处理短时间点数据集时提倡的方法，一般只有几个时间点。然而，对于长时间点数据集，每个时间点的单独测试会导致多次重复比对分析，例如每个时间点一个测试，其结果很难集成。全局分析简化了对较长时间点数据的分析，每个时间点分析保留用于对各个时间点进行特别感兴趣的重点比较。
时序数据可以是单一条件（以识别随时间变化的基因），也可以是多个条件（例如我们正在示范的流感数据集），这将随着我们感兴趣的问题而改变。
moanin中的DE分析。包moanin提供了执行这两种类型方法的功能，我们的重点是全局方法，特别是通过样条模型拟合到基因。
在这两种情况下，我们首先需要设置一个对象（moanin对象）来保存元数据，以及用于拟合样条模型的信息（公式、样条自由度数等）。moanin对象将包含整个分析过程中使用的信息，特别是每个样本的条件和时间点以及用于拟合样条模型的基础矩阵。
我们首先使用create_moanin_model函数创建moanin对象。 我们需要为函数提供两个东西：带有元数据的 data.frame，以及函数建模中使用的样条线的自由度数。元数据data.frame对象应至少包含两列：一列名为 Group，包含处理，第二列名为Timepoint，包含时间点信息）。
为我们的数据创建一个moanin对象：
moanin_model = create_moanin_model(data=data, meta=meta,
degrees_of_freedom=6)
create_moanin_model函数的主要操作除了保存样本的必要元数据之外，还为样条拟合创建基础矩阵。该矩阵给出了每个样本的所有样条基础函数的评估。默认情况下，create_moanin_model将创建样条基础函数，这将为每个组（由元数据data.frame中的Group变量定义）提供不同的样条拟合。通过以下 R 公式语法实现：
formula = ~Group:ns(Timepoint, df=degrees_of_fredoom) + Group + 0
或者，可以定义自己的公式，或者简单地自己提供基础矩阵。 当我们打印对象时我们可以看到这些信息：
show(moanin_model)
## Moanin object on 209 samples containing the following information:
## Group variable given by 'Group' with the following levels:
## M K C VL VH
## 42 42 42 42 41
## Time variable given by 'Timepoint'
## Basis matrix with 35 basis_matrix functions
## Basis matrix was constructed with the following spline_formula
## ~Group + Group:splines::ns(Timepoint, df = 6) + 0
##
## Information about the data (a SummarizedExperiment object):
## class: SummarizedExperiment
## dim: 39544 209
## metadata(0):
## assays(1): ''
## rownames(39544): NM_009912 NM_008725 ... NM_010201.1 AK078781
## rowData names(0):
## colnames(209): GSM1557140 GSM1557141 ... GSM1557347 GSM1557348
## colData names(5): '' Group Replicate Timepoint WeeklyGroup
moanin类扩展了Bioconductor的SummarizedExperiment类，以便数据、元数据和样条信息都保存在一个对象中。 这意味着对象可以像数据矩阵一样被子集化，相应的元数据和基础函数评估也将被类似地子集化：
dim(moanin_model)
## [1] 39544 209
moanin_model[1:10,1:3]
## Moanin object on 3 samples containing the following information:
## Group variable given by 'Group' with the following levels:
## M K C VL VH
## 0 3 0 0 0
## Time variable given by 'Timepoint'
## Basis matrix with 35 basis_matrix functions
## Basis matrix was constructed with the following spline_formula
## ~Group + Group:splines::ns(Timepoint, df = 6) + 0
##
## Information about the data (a SummarizedExperiment object):
## class: SummarizedExperiment
## dim: 10 3
## metadata(0):
## assays(1): ''
## rownames(10): NM_009912 NM_008725 ... NM_013782 NM_028622
## rowData names(0):
## colnames(3): GSM1557140 GSM1557141 GSM1557142
## colData names(5): '' Group Replicate Timepoint WeeklyGroup
每个时间点的差异表达分析
moanin提供了一个简单的操作来执行逐个时间点的差异表达分析。 传统上，组之间的比较是通过将组比较（在线性模型中称为对比contrasts）定义为模型系数的线性组合来完成的。 比较每个时间点内的组可以创建许多对比，因此moanin提供了以自动方式创建这些对比的功能，然后在提供的对比集上调用limma(Ritchie et al., 2015)。默认情况下，moanin需要输入RNA-Seq基因计数counts，并估计voom权重，因此对于微阵列数据，我们将设置use_voom_weights=FALSE。
在这里，我们展示了一个示例，其中我们创建了对照小鼠(“M”)和感染高剂量流感病毒株 A/Vietnam/1203/04 (H5N1)(“VL”)之间每个时间点的差异的对比。
首先，创建两组感兴趣的所有时间点的对比：
# Define contrasts
contrasts = create_timepoints_contrasts(moanin_model, "M", "VL")
这将创建一个要测试的对比特征向量，每个时间点一个，采用limma所需的格式：
contrasts[1:5]
## [1] "M.0-VL.0" "M.3-VL.3" "M.6-VL.6" "M.9-VL.9" "M.12-VL.12"
然后moanin将使用函数DE_timepoints对所有这些时间点联合运行差异表达分析。
weekly_de_analysis = DE_timepoints( moanin_model, contrasts,
use_voom_weights=FALSE)
输出是一个结果表，其中每行对应一个基因，列对应对比组的p值 (pval)、对数倍数变化(lfc)和调整后的p值(qval)；表中基因的顺序与输入数据相同。
我们将重复这一步骤，将其余三种处理中的每一种与对照(“M”)进行比较。
contrasts = create_timepoints_contrasts(moanin_model,"M", "K")
K_weekly_de_analysis = DE_timepoints(
moanin_model, contrasts,
use_voom_weights=FALSE)
contrasts = create_timepoints_contrasts(moanin_model,"M", "C")
C_weekly_de_analysis = DE_timepoints(
moanin_model, contrasts,
use_voom_weights=FALSE)
contrasts = create_timepoints_contrasts( moanin_model,"M", "VH")
VH_weekly_de_analysis = DE_timepoints(
moanin_model, contrasts,
use_voom_weights=FALSE)
在图 5 中，我们显示了每个时间点在对照和每种流感病毒株之间的差异表达基因分布。 分析结果证明了一些总体趋势，明显有更多的基因在稍后的时间点出现差异表达，并且越南高剂量可能比低剂量显示出更早的发病时间。
然而，通过每个时间点差异表达分析得到的基因数量的分布结果也体现了一些问题。 我们可以看到，一些时间点的显着差异表达的基因很少（例如川崎菌株的时间点6H和18H）。 虽然某些基因在这些时间点之间可能存在生物学差异，但在时间点3H差异表达的绝大多数基因不太可能在6H停止差异表达，然后跳回到9H进行差异表达。更可能的解释是，6H样本存在一些技术操作或生物实验误差，这些误差会产生更低的组内一致性，从而降低检测差异基因的能力。
这种方法的另一个问题是难以得到特定基因的一般时间序列表达模式。 例如，对于川崎菌株与对照的比较，有330个基因在时间点48、60、120、168H是差异表达的，但在72H不是。 其中许多基因可能没有在72H中达到统计学上的显著性，但在48H-168H之间的总体趋势中并未显示出真正的差异。 在图6中，我们显示了前10个此类基因的图。 我们可以看到，大多数基因在整个时间过程中显示出总体变化，但在72H时总体变异性增加，通常是由于单个重复表现为异常值：尽管表现出总体表达模式是不同的，这种增加的变异性导致72H时该时间点缺乏显着性。
#Code for counting the number of unique combinations of time points that are DE
getTimepoint<-function(x){sapply(strsplit(gsub("_qval","",x),"\\."),.subset2,3)}
qval_colnames = colnames(K_weekly_de_analysis)[
grepl("qval", colnames(K_weekly_de_analysis))]
signifCombos<-apply(K_weekly_de_analysis[, qval_colnames], 1,
function(x){paste(getTimepoint(qval_colnames[which(x<0.05)]),collapse=",")})
signifCombos<-signifCombos[signifCombos!=""]
tabCombos<-table(signifCombos)
exampleGenes<-names(signifCombos[signifCombos=="48,60,120,168"][1:10])
plot_splines_data(moanin_model, subset_data=exampleGenes,
colors=ann_colors$Group, smooth=TRUE)
综上所述，经典的差异表达方法适用于常规无时序的数据，但不适合处理时间序列结构的数据。
两组间时序差异表达分析
为了分析时间序列结构的数据，Storey etal.(2005)提出用样条函数对时序微阵列中的每个基因进行建模，并使用log-ratio likelihood test 来检测差异表达的基因。
moanin进一步扩展了这一想法，不仅为每个基因拟合样条函数，而且还提供了比较不同条件之间的时间序列数据的功能。通过函数DE_timecourse实现，该函数采用与DE_timepoints类似的输入，只是现在它将拟合每个基因的样条函数并测试整个平均函数（与DE_timepoints不同，因此不需要扩大到各个时间点的对比）。
# Differential expression analysis
timecourse_contrasts = c("M-K", "M-C", "M-VL", "M-VH")
# The function takes the data (data.frame or named matrix), the meta data
# (data.frame containing a timepoint and group column, the first corresponding
# to the time-course information, the latter corresponding to the
# treatment).
DE_results = DE_timecourse( moanin_model, timecourse_contrasts,
use_voom_weights=FALSE)
DE_timecourse的输出是每次比较的p值和Benjamini-Hochberg校正q值的矩阵。
names(DE_results)
## [1] "M-K_stat" "M-K_pval" "M-K_qval" "M-C_stat" "M-C_pval" "M-C_qval"
## [7] "M-VL_stat" "M-VL_pval" "M-VL_qval" "M-VH_stat" "M-VH_pval" "M-VH_qval"
为了方便起见，我们将它们分成两个矩阵。
pval_columns = colnames(DE_results)[
grepl("pval", colnames(DE_results))]
qval_columns = colnames(DE_results)[
grepl("qval", colnames(DE_results))]
pvalues = DE_results[, pval_columns]
qvalues = DE_results[, qval_columns]
得到的差异表达基因数量约为12000至29000个，具体取决于给予小鼠的流感病毒株和剂量(图7)。
时序数据的对数倍变化
经典差异表达分析的下一步通常是通过计算每个基因处理和对照之间基因表达的对数倍数变化来评估处理的效果。
人们对随时间变化的平均对数倍数变化或累积对数倍数变化感兴趣。有时，基因可能在时间过程数据开始时过度表达，然后在实验结束时表达下调。因此，moanin提供了多种可能的方法来计算整个时间过程中的对数倍数变化。这是通过函数estimate_log_fold_change实现的，该函数需要以下参数：moanin 对象、要评估的对比以及用于估计对数倍数变化的方法。
单个时间点。第一种方法（“timely”）提供了一个简单的函数来计算每个单独时间点的对数倍数变化。
log_fold_change_timepoints = estimate_log_fold_change(
moanin_model, timecourse_contrasts, method="timely")
然后可以使用该矩阵来可视化每个时间点每个对比度的对数倍变化。 累积效应。有时，每个对比的每个基因的单个值更有用，上面使用的estimate_log_fold_change也为此提供了几个选项。 但请注意，当基因不一致上调或下调时，方向的估计将无法准确代表观察到的变化。
我们在数据上演示了其中几种方法。
log_fold_change_timecourse = estimate_log_fold_change(
moanin_model, timecourse_contrasts, method="timecourse")
log_fold_change_sum = estimate_log_fold_change(
moanin_model, timecourse_contrasts, method="sum")
log_fold_change_max = estimate_log_fold_change(
moanin_model, timecourse_contrasts, method="max")
log_fold_change_min = estimate_log_fold_change(
moanin_model, timecourse_contrasts, method="min")
返回的对象是一个矩阵，其中每行对应一个基因，每列对应一个对照，每个条目对应这对对照和基因的对数倍数变化。 在图8中，我们绘制了每种方法相对于彼此的对数倍数变化情况。我们看到每种方法捕获时序数据的不同元素，例如总体变化与最大变化。
通过对数倍数变化和p值的合并，我们现在可以查看传统的火山图。 在图9中，我们显示了火山图的示例，用于使用计算对数倍数变化的“timecourse”方法来比较对照品和川崎品系。
pvalue = DE_results[, "M-K_pval"]
names(pvalue) = row.names(DE_results)
lfc_timecourse = log_fold_change_timecourse[, "M-K"]
names(lfc_timecourse) = row.names(log_fold_change_timecourse)
plot(lfc_timecourse, -log10(pvalue), pch=20, main="Volcano plot",
xlim=c(-2.5, 2), xlab="Timecourse lfc")
可视化感兴趣的基因
moanin包还提供了一个简单的实用函数(plot_splines_data)来可视化不同条件下的基因时间过程数据。 在图10中，我们绘制了p值最小的10个基因。
top_DE_genes_pval = names(sort(pvalue)[1:10])
plot_splines_data(moanin_model,subset_data=top_DE_genes_pval,
colors=ann_colors$Group, smooth=TRUE,
mar=c(1.5,2.5,2,0.1))
图11，可视化具有最大绝对时序对数倍变化的基因。
top_DE_genes_lfc = names(
sort(abs(lfc_timecourse),
decreasing=TRUE)[1:10])
plot_splines_data(moanin_model, subset_data=top_DE_genes_lfc,
colors=ann_colors$Group, smooth=TRUE,
mar=c(1.5,2.5,2,0.1))
在分析这些可视化结果时，我们可以看到基因通常遵循相似的表达模式，尽管每个基因的变化幅度不同。 我们可以利用这一观察结果将基因聚类成具有相似转录组反应模式的clusters。
这就是我们下一篇将要展示的内容：RNA-Seq 分析流程：多时间点样本分析实战（二）。
敬请期待，反馈越热烈，更新越快噢！
