RNA-Seq 分析流程：多時(shí)間點(diǎn)樣本分析實(shí)戰(zhàn)（一）

小魚(yú)爸

2023-07-12

簡(jiǎn)介

了解基因表達(dá)的時(shí)序動(dòng)態(tài)變化是生物學(xué)的一個(gè)基本問(wèn)題惜纸。此教程提供了多時(shí)間點(diǎn)數(shù)據(jù)的分步實(shí)戰(zhàn)流程：(1)數(shù)據(jù)集的質(zhì)量控制和標(biāo)準(zhǔn)化宪赶；(2)進(jìn)行差異表達(dá)分析子姜；(3)時(shí)序數(shù)據(jù)的聚類(lèi)巷屿；(4)用GO term和KEGG通路富集分析解釋聚類(lèi)簇漱贱。作為實(shí)戰(zhàn)流程治笨，我們應(yīng)用的數(shù)據(jù)為暴露于四種流感病毒株的小鼠的多時(shí)間點(diǎn)轉(zhuǎn)錄組數(shù)據(jù)砌们。

圖1:用于分析時(shí)間過(guò)程數(shù)據(jù)集的工作流程。

本分析流程基于先前建立的用于發(fā)育時(shí)間序列的新策略催烘。它依賴(lài)于幾個(gè)用于分析基因表達(dá)數(shù)據(jù)的標(biāo)準(zhǔn)包沥阱，其中一些特定用于多時(shí)間點(diǎn)數(shù)據(jù)。各步驟使用的主要函數(shù)在moaninR包中伊群。

安裝和設(shè)置包

moanin和timecoursedata可從bioconductor獲得考杉，并且可以使用BiocManager包中的install功能進(jìn)行安裝：

library (BiocManager)
BiocManager::install("timecoursedata")
#如果下載失敗，可以去Bioconductor官網(wǎng)（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，則可以通過(guò)以下方式安裝上面的CRAN和Bioconductor軟件包：

BiocManager::install(
     c("NMF", "ggfortify", "topGO", "biomaRt",
       "BiocWorkflowTools", "timecoursedata", ))

小鼠肺組織對(duì)流感的動(dòng)態(tài)反應(yīng)分析

該分析流程使用微陣列時(shí)序?qū)嶒?yàn)的數(shù)據(jù)進(jìn)行示范蔽午，將小鼠暴露于三種不同的流感毒株易茬，并在感染后14個(gè)時(shí)間點(diǎn)（0、3、6抽莱、9范抓、12、18食铐、24匕垫、30、36虐呻、48象泵、60 小時(shí)， 3斟叼、5 和 7 天后）收集肺組織偶惠。研究中使用的三種流感毒株是(1)低致病性季節(jié)性H1N1流感病毒(A/Kawasaki/UTK4/2009 [H1N1])；(2)2009年大流行季節(jié)的一種輕度致病性病毒(A/California/04/2009 [H1N1])朗涩；(3)高致病性H5N1禽流感病毒(A/Vietnam/1203/2004[H5N1])忽孽。給小鼠注射了每種病毒10⁵ PFU。(4)另外谢床，42只小鼠注射了較低劑量的越南禽流感病毒(10³ PFU)兄一。

通過(guò)將基因表達(dá)時(shí)序數(shù)據(jù)與病毒生長(zhǎng)數(shù)據(jù)相結(jié)合，作者發(fā)現(xiàn)识腿，肺組織的炎癥反應(yīng)是門(mén)控的（gated）出革，肺部的病毒濃度存在閾值（threshold）。一旦超過(guò)這個(gè)閾值渡讼，就會(huì)產(chǎn)生強(qiáng)烈的炎癥和細(xì)胞因子骂束。這些結(jié)果證明病理反應(yīng)受病毒濃度的非線性調(diào)節(jié)。 此外硝全，Varoquaux等人利用類(lèi)似的步驟來(lái)分析作物S. bicolor對(duì)干旱的終生轉(zhuǎn)錄組響應(yīng)的RNA-seq數(shù)據(jù)栖雾。

數(shù)據(jù)概覽

首先讓我們加載數(shù)據(jù)。moanin包包含(Shoemaker et al., 2015)的標(biāo)準(zhǔn)化數(shù)據(jù)和元數(shù)據(jù)伟众。

# Now load in the metadata
data(shoemaker2015)
meta = shoemaker2015$meta
data = shoemaker2015$data

元數(shù)據(jù)meta包含有關(guān)組、重復(fù)和每次觀察的時(shí)間點(diǎn)的信息召廷。

head(meta)

在深入探索性分析和質(zhì)量控制之前凳厢，我們先分配整個(gè)分析中使用配色方案。我們將組和時(shí)間點(diǎn)的配色方案定義為向量vectors竞慢。我們還定義了一系列標(biāo)記（或繪圖符號(hào)）來(lái)區(qū)分散點(diǎn)圖中的重復(fù)樣本先紫。我們還對(duì)描述處理的組進(jìn)行了重新排序，以便按致病性從低到高排序（對(duì)照在前）筹煮。

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)])

質(zhì)量控制和標(biāo)準(zhǔn)化

基因表達(dá)數(shù)據(jù)分析的第一步始終是對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化和質(zhì)量控制檢查遮精。通常，在歸一化之前和之后執(zhí)行兩個(gè)質(zhì)量控制和探索性分析步驟：(1)樣本的降維分析；(2)各樣本之間的相關(guān)性分析本冲。在這兩種情況下准脂，我們期望產(chǎn)生高表達(dá)的生物信號(hào)，而重復(fù)樣本應(yīng)該強(qiáng)烈聚集或彼此相關(guān)檬洞。 在進(jìn)行下游分析之前狸膏，我們只保留了高度可變的基因：在這一步，我們只保留前 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)

圖2：方差分布

讓我們首先進(jìn)行主成分分析（PCA）分析广料。

pca_data = prcomp(t(variance_filtered_data), rank=3, center=TRUE, scale=TRUE) 
percent_var = round(100 * attr(pca_data, "percentVar"))

在圖3中砾脑，我們通過(guò)繪制主成分(PC)來(lái)可視化數(shù)據(jù)，樣本按其條件或采樣時(shí)間著色艾杏，并且每個(gè)樣本復(fù)制用不同的符號(hào)區(qū)分韧衣。我們可以看到后面的時(shí)間點(diǎn)有很大的差異。

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)

圖3:主成分圖(PC)：PC2與PC1(右)和PC2與PC3(左)糜颠。

我們還在圖4中將每個(gè)樣本之間的Pearson相關(guān)性繪制為熱圖（使用NMF中的函數(shù)aheatmap）汹族。我們按組（處理）和時(shí)間點(diǎn)（采樣時(shí)間）對(duì)樣本進(jìn)行排序。

# 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")

圖 4:樣本之間相關(guān)性的熱圖

我們已經(jīng)可以從相關(guān)圖中看到有趣的結(jié)果其兴。首先顶瞒，從對(duì)照小鼠采集的樣本之間的互相關(guān)性高于其余治療組之間的互相關(guān)性。其次元旬，感染流感的小鼠直到時(shí)間點(diǎn)36才會(huì)有輕微反應(yīng)榴徐。第三，菌株的致病性越低匀归，樣品越接近對(duì)照條件坑资。第四，時(shí)間點(diǎn)120和168的越南樣本是與對(duì)照樣本差異最大的樣本穆端。

時(shí)序數(shù)據(jù)的差異表達(dá) (DE) 分析

時(shí)序數(shù)據(jù)的 DE 分析方法

基因表達(dá)分析的下一步通常是運(yùn)行差異表達(dá)分析袱贮，通常是為了找到不同條件下的差異基因。 對(duì)于時(shí)序數(shù)據(jù)体啰，有兩種不同的方法來(lái)確定差異表達(dá)基因攒巍，

1）按時(shí)間點(diǎn)分析，我們將每個(gè)時(shí)間點(diǎn)視為不同的條件荒勇，并確定哪些基因在特定時(shí)間點(diǎn)之間或單個(gè)時(shí)間點(diǎn)的特定條件之間發(fā)生變化柒莉。

2）全局分析，我們考察隨著時(shí)間的推移全局的表達(dá)模式沽翔，并推斷哪些基因在不同條件下具有不同的模式兢孝，或者隨著時(shí)間的推移具有變化的模式。該方法的第一步是為每個(gè)基因擬合一個(gè)樣條模型(Spline model, Storey et al., 2005)，然后使用該樣條模型來(lái)測(cè)試隨時(shí)間變化的不同類(lèi)型的差異表達(dá)跨蟹。

按時(shí)間點(diǎn)分析使用經(jīng)典的差異表達(dá)方法雳殊，并且通常是處理短時(shí)間點(diǎn)數(shù)據(jù)集時(shí)提倡的方法，一般只有幾個(gè)時(shí)間點(diǎn)喷市。然而相种，對(duì)于長(zhǎng)時(shí)間點(diǎn)數(shù)據(jù)集，每個(gè)時(shí)間點(diǎn)的單獨(dú)測(cè)試會(huì)導(dǎo)致多次重復(fù)比對(duì)分析品姓，例如每個(gè)時(shí)間點(diǎn)一個(gè)測(cè)試寝并，其結(jié)果很難集成。全局分析簡(jiǎn)化了對(duì)較長(zhǎng)時(shí)間點(diǎn)數(shù)據(jù)的分析腹备，每個(gè)時(shí)間點(diǎn)分析保留用于對(duì)各個(gè)時(shí)間點(diǎn)進(jìn)行特別感興趣的重點(diǎn)比較衬潦。

時(shí)序數(shù)據(jù)可以是單一條件（以識(shí)別隨時(shí)間變化的基因），也可以是多個(gè)條件（例如我們正在示范的流感數(shù)據(jù)集）植酥，這將隨著我們感興趣的問(wèn)題而改變镀岛。

moanin中的DE分析。包moanin提供了執(zhí)行這兩種類(lèi)型方法的功能友驮，我們的重點(diǎn)是全局方法漂羊，特別是通過(guò)樣條模型擬合到基因。

在這兩種情況下卸留，我們首先需要設(shè)置一個(gè)對(duì)象（moanin對(duì)象）來(lái)保存元數(shù)據(jù)走越，以及用于擬合樣條模型的信息（公式、樣條自由度數(shù)等）耻瑟。moanin對(duì)象將包含整個(gè)分析過(guò)程中使用的信息旨指，特別是每個(gè)樣本的條件和時(shí)間點(diǎn)以及用于擬合樣條模型的基礎(chǔ)矩陣。

我們首先使用create_moanin_model函數(shù)創(chuàng)建moanin對(duì)象喳整。 我們需要為函數(shù)提供兩個(gè)東西：帶有元數(shù)據(jù)的 data.frame谆构，以及函數(shù)建模中使用的樣條線的自由度數(shù)。元數(shù)據(jù)data.frame對(duì)象應(yīng)至少包含兩列：一列名為 Group框都，包含處理搬素，第二列名為Timepoint，包含時(shí)間點(diǎn)信息）魏保。

為我們的數(shù)據(jù)創(chuàng)建一個(gè)moanin對(duì)象：

moanin_model = create_moanin_model(data=data, meta=meta,
    degrees_of_freedom=6)

create_moanin_model函數(shù)的主要操作除了保存樣本的必要元數(shù)據(jù)之外蔗蹋，還為樣條擬合創(chuàng)建基礎(chǔ)矩陣。該矩陣給出了每個(gè)樣本的所有樣條基礎(chǔ)函數(shù)的評(píng)估囱淋。默認(rèn)情況下，create_moanin_model將創(chuàng)建樣條基礎(chǔ)函數(shù)餐塘，這將為每個(gè)組（由元數(shù)據(jù)data.frame中的Group變量定義）提供不同的樣條擬合妥衣。通過(guò)以下 R 公式語(yǔ)法實(shí)現(xiàn)：

formula = ~Group:ns(Timepoint, df=degrees_of_fredoom) + Group + 0

或者，可以定義自己的公式，或者簡(jiǎn)單地自己提供基礎(chǔ)矩陣税手。 當(dāng)我們打印對(duì)象時(shí)我們可以看到這些信息：

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類(lèi)擴(kuò)展了Bioconductor的SummarizedExperiment類(lèi)蜂筹，以便數(shù)據(jù)、元數(shù)據(jù)和樣條信息都保存在一個(gè)對(duì)象中芦倒。 這意味著對(duì)象可以像數(shù)據(jù)矩陣一樣被子集化艺挪，相應(yīng)的元數(shù)據(jù)和基礎(chǔ)函數(shù)評(píng)估也將被類(lèi)似地子集化：

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

每個(gè)時(shí)間點(diǎn)的差異表達(dá)分析

moanin提供了一個(gè)簡(jiǎn)單的操作來(lái)執(zhí)行逐個(gè)時(shí)間點(diǎn)的差異表達(dá)分析。 傳統(tǒng)上兵扬，組之間的比較是通過(guò)將組比較（在線性模型中稱(chēng)為對(duì)比contrasts）定義為模型系數(shù)的線性組合來(lái)完成的麻裳。 比較每個(gè)時(shí)間點(diǎn)內(nèi)的組可以創(chuàng)建許多對(duì)比，因此moanin提供了以自動(dòng)方式創(chuàng)建這些對(duì)比的功能器钟，然后在提供的對(duì)比集上調(diào)用limma(Ritchie et al., 2015)津坑。默認(rèn)情況下，moanin需要輸入RNA-Seq基因計(jì)數(shù)counts傲霸，并估計(jì)voom權(quán)重疆瑰，因此對(duì)于微陣列數(shù)據(jù)，我們將設(shè)置use_voom_weights=FALSE昙啄。

在這里穆役，我們展示了一個(gè)示例，其中我們創(chuàng)建了對(duì)照小鼠(“M”)和感染高劑量流感病毒株 A/Vietnam/1203/04 (H5N1)(“VL”)之間每個(gè)時(shí)間點(diǎn)的差異的對(duì)比梳凛。

首先耿币，創(chuàng)建兩組感興趣的所有時(shí)間點(diǎn)的對(duì)比：

# Define contrasts  
contrasts = create_timepoints_contrasts(moanin_model, "M", "VL")

這將創(chuàng)建一個(gè)要測(cè)試的對(duì)比特征向量，每個(gè)時(shí)間點(diǎn)一個(gè)伶跷，采用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將使用函數(shù)DE_timepoints對(duì)所有這些時(shí)間點(diǎn)聯(lián)合運(yùn)行差異表達(dá)分析掰读。

weekly_de_analysis = DE_timepoints( moanin_model, contrasts,
     use_voom_weights=FALSE)

輸出是一個(gè)結(jié)果表，其中每行對(duì)應(yīng)一個(gè)基因叭莫，列對(duì)應(yīng)對(duì)比組的p值 (pval)蹈集、對(duì)數(shù)倍數(shù)變化(lfc)和調(diào)整后的p值(qval)；表中基因的順序與輸入數(shù)據(jù)相同雇初。

我們將重復(fù)這一步驟拢肆，將其余三種處理中的每一種與對(duì)照(“M”)進(jìn)行比較。

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 中靖诗，我們顯示了每個(gè)時(shí)間點(diǎn)在對(duì)照和每種流感病毒株之間的差異表達(dá)基因分布郭怪。 分析結(jié)果證明了一些總體趨勢(shì)，明顯有更多的基因在稍后的時(shí)間點(diǎn)出現(xiàn)差異表達(dá)刊橘，并且越南高劑量可能比低劑量顯示出更早的發(fā)病時(shí)間鄙才。

圖5：將流感毒株與對(duì)照進(jìn)行比較時(shí)，在每個(gè)時(shí)間點(diǎn)發(fā)現(xiàn) DE 的基因數(shù)量促绵。

然而攒庵，通過(guò)每個(gè)時(shí)間點(diǎn)差異表達(dá)分析得到的基因數(shù)量的分布結(jié)果也體現(xiàn)了一些問(wèn)題嘴纺。 我們可以看到，一些時(shí)間點(diǎn)的顯著差異表達(dá)的基因很少（例如川崎菌株的時(shí)間點(diǎn)6H和18H）浓冒。 雖然某些基因在這些時(shí)間點(diǎn)之間可能存在生物學(xué)差異栽渴，但在時(shí)間點(diǎn)3H差異表達(dá)的絕大多數(shù)基因不太可能在6H停止差異表達(dá)，然后跳回到9H進(jìn)行差異表達(dá)稳懒。更可能的解釋是闲擦，6H樣本存在一些技術(shù)操作或生物實(shí)驗(yàn)誤差，這些誤差會(huì)產(chǎn)生更低的組內(nèi)一致性场梆，從而降低檢測(cè)差異基因的能力墅冷。

這種方法的另一個(gè)問(wèn)題是難以得到特定基因的一般時(shí)間序列表達(dá)模式。 例如辙谜，對(duì)于川崎菌株與對(duì)照的比較俺榆，有330個(gè)基因在時(shí)間點(diǎn)48、60装哆、120罐脊、168H是差異表達(dá)的，但在72H不是蜕琴。 其中許多基因可能沒(méi)有在72H中達(dá)到統(tǒng)計(jì)學(xué)上的顯著性萍桌，但在48H-168H之間的總體趨勢(shì)中并未顯示出真正的差異。 在圖6中凌简，我們顯示了前10個(gè)此類(lèi)基因的圖上炎。 我們可以看到，大多數(shù)基因在整個(gè)時(shí)間過(guò)程中顯示出總體變化雏搂，但在72H時(shí)總體變異性增加藕施，通常是由于單個(gè)重復(fù)表現(xiàn)為異常值：盡管表現(xiàn)出總體表達(dá)模式是不同的，這種增加的變異性導(dǎo)致72H時(shí)該時(shí)間點(diǎn)缺乏顯著性凸郑。

#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)

圖6：10個(gè)基因隨時(shí)間的表達(dá)值在時(shí)間點(diǎn)48裳食、60、120芙沥、168H發(fā)現(xiàn)DE诲祸，但在72H時(shí)間點(diǎn)不存在顯著DE。

綜上所述而昨，經(jīng)典的差異表達(dá)方法適用于常規(guī)無(wú)時(shí)序的數(shù)據(jù)救氯，但不適合處理時(shí)間序列結(jié)構(gòu)的數(shù)據(jù)。

兩組間時(shí)序差異表達(dá)分析

為了分析時(shí)間序列結(jié)構(gòu)的數(shù)據(jù)歌憨，Storey etal.(2005)提出用樣條函數(shù)對(duì)時(shí)序微陣列中的每個(gè)基因進(jìn)行建模着憨，并使用log-ratio likelihood test 來(lái)檢測(cè)差異表達(dá)的基因。

moanin進(jìn)一步擴(kuò)展了這一想法务嫡，不僅為每個(gè)基因擬合樣條函數(shù)享扔，而且還提供了比較不同條件之間的時(shí)間序列數(shù)據(jù)的功能底桂。通過(guò)函數(shù)DE_timecourse實(shí)現(xiàn)，該函數(shù)采用與DE_timepoints類(lèi)似的輸入惧眠，只是現(xiàn)在它將擬合每個(gè)基因的樣條函數(shù)并測(cè)試整個(gè)平均函數(shù)（與DE_timepoints不同，因此不需要擴(kuò)大到各個(gè)時(shí)間點(diǎn)的對(duì)比）于个。

# 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"

為了方便起見(jiàn)，我們將它們分成兩個(gè)矩陣厅篓。

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]

得到的差異表達(dá)基因數(shù)量約為12000至29000個(gè)秀存，具體取決于給予小鼠的流感病毒株和劑量(圖7)。

圖7：基于全局樣條分析的每種條件的 DE 基因分布

時(shí)序數(shù)據(jù)的對(duì)數(shù)倍變化

經(jīng)典差異表達(dá)分析的下一步通常是通過(guò)計(jì)算每個(gè)基因處理和對(duì)照之間基因表達(dá)的對(duì)數(shù)倍數(shù)變化來(lái)評(píng)估處理的效果羽氮。

人們對(duì)隨時(shí)間變化的平均對(duì)數(shù)倍數(shù)變化或累積對(duì)數(shù)倍數(shù)變化感興趣或链。有時(shí)，基因可能在時(shí)間過(guò)程數(shù)據(jù)開(kāi)始時(shí)過(guò)度表達(dá)档押，然后在實(shí)驗(yàn)結(jié)束時(shí)表達(dá)下調(diào)澳盐。因此，moanin提供了多種可能的方法來(lái)計(jì)算整個(gè)時(shí)間過(guò)程中的對(duì)數(shù)倍數(shù)變化令宿。這是通過(guò)函數(shù)estimate_log_fold_change實(shí)現(xiàn)的叼耙，該函數(shù)需要以下參數(shù)：moanin 對(duì)象、要評(píng)估的對(duì)比以及用于估計(jì)對(duì)數(shù)倍數(shù)變化的方法粒没。

單個(gè)時(shí)間點(diǎn)筛婉。第一種方法（“timely”）提供了一個(gè)簡(jiǎn)單的函數(shù)來(lái)計(jì)算每個(gè)單獨(dú)時(shí)間點(diǎn)的對(duì)數(shù)倍數(shù)變化。

log_fold_change_timepoints = estimate_log_fold_change(
    moanin_model, timecourse_contrasts,  method="timely")

然后可以使用該矩陣來(lái)可視化每個(gè)時(shí)間點(diǎn)每個(gè)對(duì)比度的對(duì)數(shù)倍變化癞松。 累積效應(yīng)爽撒。有時(shí)，每個(gè)對(duì)比的每個(gè)基因的單個(gè)值更有用响蓉，上面使用的estimate_log_fold_change也為此提供了幾個(gè)選項(xiàng)硕勿。 但請(qǐng)注意，當(dāng)基因不一致上調(diào)或下調(diào)時(shí)厕妖，方向的估計(jì)將無(wú)法準(zhǔn)確代表觀察到的變化首尼。

我們?cè)跀?shù)據(jù)上演示了其中幾種方法。

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")

返回的對(duì)象是一個(gè)矩陣言秸，其中每行對(duì)應(yīng)一個(gè)基因软能，每列對(duì)應(yīng)一個(gè)對(duì)照，每個(gè)條目對(duì)應(yīng)這對(duì)對(duì)照和基因的對(duì)數(shù)倍數(shù)變化举畸。 在圖8中查排，我們繪制了每種方法相對(duì)于彼此的對(duì)數(shù)倍數(shù)變化情況。我們看到每種方法捕獲時(shí)序數(shù)據(jù)的不同元素抄沮，例如總體變化與最大變化跋核。

圖8：計(jì)算每個(gè)基因?qū)?shù)倍變化的不同方法的比較岖瑰。

通過(guò)對(duì)數(shù)倍數(shù)變化和p值的合并，我們現(xiàn)在可以查看傳統(tǒng)的火山圖砂代。 在圖9中蹋订，我們顯示了火山圖的示例，用于使用計(jì)算對(duì)數(shù)倍數(shù)變化的“timecourse”方法來(lái)比較對(duì)照品和川崎品系刻伊。

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")

圖9：川崎菌株與對(duì)照菌株比較的火山圖露戒，其中倍數(shù)變化是用時(shí)序方法計(jì)算的。

可視化感興趣的基因

moanin包還提供了一個(gè)簡(jiǎn)單的實(shí)用函數(shù)(plot_splines_data)來(lái)可視化不同條件下的基因時(shí)間過(guò)程數(shù)據(jù)捶箱。 在圖10中智什，我們繪制了p值最小的10個(gè)基因。

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))

圖10：p 值最小的前 10 個(gè)基因

圖11丁屎，可視化具有最大絕對(duì)時(shí)序?qū)?shù)倍變化的基因荠锭。

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))

圖11：時(shí)序絕對(duì)對(duì)數(shù)倍變化最大的前10個(gè)基因

在分析這些可視化結(jié)果時(shí)，我們可以看到基因通常遵循相似的表達(dá)模式晨川，盡管每個(gè)基因的變化幅度不同证九。 我們可以利用這一觀察結(jié)果將基因聚類(lèi)成具有相似轉(zhuǎn)錄組反應(yīng)模式的clusters。

這就是我們下一篇將要展示的內(nèi)容：RNA-Seq 分析流程：多時(shí)間點(diǎn)樣本分析實(shí)戰(zhàn)（二）础爬。
敬請(qǐng)期待甫贯，反饋越熱烈，更新越快噢看蚜！