最近開發(fā)了一個R包系羞,enONE(https://github.com/thereallda/enONE),可以用于NAD-RNA-seq或者其他代謝物修飾RNA的數(shù)據(jù)分析中霸琴。相關(guān)工作也發(fā)表了椒振,歡迎大家使用并提出意見 :)
Epitranscriptome analysis of NAD-capped RNA by spike-in-based
normalization and prediction of chronological age
(https://doi.org/10.1016/j.isci.2023.108558)
前言
細胞內(nèi)的核苷類代謝物,如煙酰胺腺嘌呤二核苷酸(NAD)沈贝,能夠被添加到mRNA
5’末端杠人,形成NAD帽修飾的RNA(NAD-RNA);這類非經(jīng)典起始核苷酸(NCIN)修飾的RNA被統(tǒng)稱為NCIN帽修飾的RNA(NCIN-RNA)宋下。NCIN-RNA天然地將代謝物與基因表達緊密地聯(lián)系起來嗡善,定義了一種新穎但功能未知的表觀轉(zhuǎn)錄修飾。NCIN-RNA的表征是功能研究的前提学歧。當前NCIN-RNA的表征方法罩引,如ONE-seq和DO-seq,實現(xiàn)了代謝物帽修飾RNA在細胞和組織中的鑒定枝笨。然而袁铐,這些方法產(chǎn)出的數(shù)據(jù)會受到源于化學酶促標記及親和富集反應等實驗操作誤差的影響。因此横浑,這里提出了enONE
這種NAD-RNA表觀轉(zhuǎn)錄組的計算分析方法剔桨,實現(xiàn)NAD-RNA數(shù)據(jù)的定量和比較分析。
enONE
的R包可以在GitHub中找到(https://github.com/thereallda/enONE)
Setup
在這個教程中徙融,我們會使用人類外周血單核細胞(PBMCs)的NAD-RNA-seq數(shù)據(jù)來展示
enONE
的工作流程洒缀。
要注意的是,這里混入了三種spike-in RNAs:
- 來自模式生物黑腹果蠅的總RNA欺冀,用于
enONE
的normalization树绩; - Synthetic RNA,其中含有5%的NAD修飾形式(相對于m7G-RNA)隐轩,用于確定富集靈敏度饺饭;
- Synthetic RNA,其為100% m7G加帽的形式职车,用于確定富集特異性瘫俊。
Figure: Schematic workflow for total RNAs from PBMCs and three sets of
spike-ins.
我們首先讀取數(shù)據(jù),包括基因表達量的count matrix和樣本信息的metadata悴灵。matrix和樣本信息的军援。metadata至少包括兩列,即用于樣本生物分組的”condition”列和用于樣本富集分組的”enrich”列称勋。
library(enONE)
library(tidyverse)
library(patchwork)
# read in metadata and counts matrix
counts_mat <- read.csv("D:/Academic/data/Counts.csv", row.names = 1)
meta <- read.csv("D:/Academic/data/metadata.csv", comment.char = "#")
head(meta)
#> id sample.id condition batch enrich
#> 1 D1 Y2 Y.Input D Input
#> 2 D2 Y10 Y.Input D Input
#> 3 D3 Y14 Y.Input D Input
#> 4 D4 Y16 Y.Input D Input
#> 5 D5 M13 M.Input D Input
#> 6 D6 M14 M.Input D Input
# rownames of metadata should be consistent with the colnames of counts_mat
rownames(meta) <- meta$id
另外胸哥,metadata也可以包括其他信息,例如樣本的批次信息赡鲜。
接下來空厌,我們使用count matrix和metadata創(chuàng)建一個 Enone
對象。
Enone
對象银酬,包含了NAD-RNA-seq數(shù)據(jù)集的數(shù)據(jù)(raw and normalized counts data)和結(jié)果信息(例如嘲更,數(shù)據(jù)標準化性能的打分和NAD-RNA富集結(jié)果)。
在構(gòu)建 Enone 對象時揩瞪,我們可以提供以下參數(shù):
- spike-in 基因的前綴赋朦,例如果蠅Flybase基因id前綴為FB (
spike.in.prefix = "^FB"
); - synthetic spike-in的id (
synthetic.id = c("Syn1", "Syn2"
)),這需要與counts_mat
中的行名相同宠哄; - Input樣本分組的id (
input.id = "Input"
) 和Enrichment樣本分組的id (enrich.id = "Enrich"
)壹将,與enrich
列相同。
其中, “Syn1”代表5% NAD-RNA spike-in and “Syn2”代表100% m7G-RNA spike-in.
synthetic.id
是可選的參數(shù).
# prefix of Drosophila spike-in genes
spikeInPrefix <- "^FB"
# create Enone
Enone <- createEnone(data = counts_mat,
col.data = meta,
spike.in.prefix = spikeInPrefix,
synthetic.id = c("Syn1", "Syn2"),
input.id = "Input",
enrich.id = "Enrich"
)
Enone
#> class: Enone
#> dim: 76290 62
#> metadata(0):
#> assays(1): ''
#> rownames(76290): ENSG00000223972 ENSG00000227232 ... Syn1 Syn2
#> rowData names(3): GeneID SpikeIn Synthetic
#> colnames(62): D1 D2 ... F22 F24
#> colData names(7): id id.1 ... enrich replicate
標準流程
enONE
流程包括以下四步:
- 質(zhì)量控制(Quality control);
- 基因集選让怠(Gene set selection);
- 數(shù)據(jù)校正(Normalization procedures);
- 數(shù)據(jù)校正效果評估(Normalization performance assessment).
質(zhì)量控制(Quality control)
enONE
首先進行質(zhì)量控制步驟诽俯,通過 FilterLowExprGene
保留至少在 n
個樣本中至少具有 min.count
的基因來完成此步驟。n
由 group
中指定的最小樣本組的大小確定承粤。
Enone <- FilterLowExprGene(Enone, group = Enone$condition, min.count = 20)
Enone
#> class: Enone
#> dim: 27661 62
#> metadata(0):
#> assays(1): ''
#> rownames(27661): ENSG00000227232 ENSG00000268903 ... Syn1 Syn2
#> rowData names(3): GeneID SpikeIn Synthetic
#> colnames(62): D1 D2 ... F22 F24
#> colData names(7): id id.1 ... enrich replicate
此外暴区,可以通過OutlierTest
對異常樣本進行評估并進一步移除(return=TRUE
)。由于沒有樣本被標記為異常值辛臊,因此我們將在后續(xù)分析中使用所有樣本仙粱。
## ronser"s test for outlier assessment
OutlierTest(Enone, return=FALSE)
#> Rosner's outlier test
#> i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
#> 1 0 -1.953993e-14 51.77231 80.37523 15 1.552475 3.212165 FALSE
#> 2 1 -1.317627e+00 51.14304 78.20437 19 1.554894 3.205977 FALSE
#> 3 2 -2.642993e+00 50.50717 70.38990 41 1.445991 3.199662 FALSE
運行enONE
enONE
函數(shù)可以進行基因選擇、數(shù)據(jù)標準化和標準化效果評估彻舰。該函數(shù)在 rowData
中存儲選擇的基因集伐割,在 counts slot中返回標準化計數(shù)矩陣(當return.norm=TRUE
時),在 enone_factor
slot中返回標準化因子淹遵,在 enone_metrics
和 enone_scores
slots中返回評估指標和得分口猜。以下是這些步驟的詳細描述。
Enone <- enONE(Enone,
scaling.method = c("TC", "UQ", "TMM", "DESeq", "PossionSeq"),
ruv.norm = TRUE, ruv.k = 3,
eval.pam.k = 2:6, eval.pc.n = 3,
return.norm = TRUE
)
#> Gene set selection for normalization and assessment...
#> - The number of negative control genes for normalization: 1000
#> - Estimate dispersion & Fit GLM...
#> - Testing differential genes...
#> - The number of positive evaluation genes: 500
#> - Estimate dispersion & Fit GLM...
#> - Testing differential genes...
#> - The number of negative evaluation genes: 500
#> - Estimate dispersion & Fit GLM...
#> - Testing differential genes...
#> Apply normalization...
#> - Scaling...
#> - Regression-based normalization...
#> Perform assessment...
基因集選韧复А(Gene set selection)
enONE
定義了三組基因济炎,包括:
- 負對照(
NegControl
):默認情況下,enONE
將 果蠅 spike-ins(或其他外源生物的RNA spike-in)中 FDR 排名最低的 1,000 個基因作為負對照(不富集的基因)辐真,用于校正無關(guān)誤差须尚。 - 負評估(
NegEvaluation
):默認情況下,enONE
將樣本中 FDR 排名最低的500 個基因作為負評估基因侍咱,用于評估無關(guān)誤差耐床。 - 正評估(
PosEvaluation
):默認情況下,enONE
將樣本中 FDR 排名最高的500 個基因作為正評估基因楔脯,用于評估研究感興趣的差異撩轰。
基因集的選擇可以在 enONE
函數(shù)中使用 auto=TRUE
參數(shù)(默認值)自動定義,也可以在 neg.control
, pos.eval
, neg.eval
參數(shù)中分別提供昧廷。
可以通過getGeneSet
函數(shù)和基因集的名稱(即"NegControl"
堪嫂、"NegEvaluation"
、"PosEvaluation"
)來獲取相應的基因集木柬。
getGeneSet(Enone, name = "NegControl")[1:5]
#> [1] "FBgn0031247" "FBgn0031255" "FBgn0004583" "FBgn0031324" "FBgn0026397"
數(shù)據(jù)校正(Normalization)
enONE
整合了全局縮放(global scaling)和基于回歸的校正方法(regression-based normalization)來產(chǎn)生normalization方法皆串。
對于全局縮放方法,enONE
整合了五種方法眉枕,包括:
Total Count (TC);
Upper-Quartile (UQ);
Trimmed Mean of M Values (TMM);
DESeq;
PossionSeq.
默認情況下恶复,enONE
使用所有縮放方法怜森,但用戶可以在 scaling.method
參數(shù)中選擇要使用的縮放程序。
對于regression-based方法谤牡,enONE
利用了三種RUV的方法副硅,包括:
RUVg;
RUVs;
RUVse.
例如,您可以通過選擇 ruv.norm=TRUE, ruv.k=2
來執(zhí)行使用前兩個無關(guān)誤差因子的 RUV拓哟。
RUVse 是類似于 RUVs 的方法想许。它基于每個富集樣本或背景樣本組的重復樣本中的負對照基因來估計無關(guān)誤差因子伶授,其假設富集效應在重復之間是相對穩(wěn)定的断序。
所有使用的normalization方法都可以使用 listNormalization
來提取。
listNormalization(Enone)
#> [1] "TC" "UQ" "TMM"
#> [4] "DESeq" "PossionSeq" "Raw"
#> [7] "Raw_RUVg_k1" "Raw_RUVg_k2" "Raw_RUVg_k3"
#> [10] "Raw_RUVs_k1" "Raw_RUVs_k2" "Raw_RUVs_k3"
#> [13] "Raw_RUVse_k1" "Raw_RUVse_k2" "Raw_RUVse_k3"
#> [16] "TC_RUVg_k1" "TC_RUVg_k2" "TC_RUVg_k3"
#> [19] "TC_RUVs_k1" "TC_RUVs_k2" "TC_RUVs_k3"
#> [22] "TC_RUVse_k1" "TC_RUVse_k2" "TC_RUVse_k3"
#> [25] "UQ_RUVg_k1" "UQ_RUVg_k2" "UQ_RUVg_k3"
#> [28] "UQ_RUVs_k1" "UQ_RUVs_k2" "UQ_RUVs_k3"
#> [31] "UQ_RUVse_k1" "UQ_RUVse_k2" "UQ_RUVse_k3"
#> [34] "TMM_RUVg_k1" "TMM_RUVg_k2" "TMM_RUVg_k3"
#> [37] "TMM_RUVs_k1" "TMM_RUVs_k2" "TMM_RUVs_k3"
#> [40] "TMM_RUVse_k1" "TMM_RUVse_k2" "TMM_RUVse_k3"
#> [43] "DESeq_RUVg_k1" "DESeq_RUVg_k2" "DESeq_RUVg_k3"
#> [46] "DESeq_RUVs_k1" "DESeq_RUVs_k2" "DESeq_RUVs_k3"
#> [49] "DESeq_RUVse_k1" "DESeq_RUVse_k2" "DESeq_RUVse_k3"
#> [52] "PossionSeq_RUVg_k1" "PossionSeq_RUVg_k2" "PossionSeq_RUVg_k3"
#> [55] "PossionSeq_RUVs_k1" "PossionSeq_RUVs_k2" "PossionSeq_RUVs_k3"
#> [58] "PossionSeq_RUVse_k1" "PossionSeq_RUVse_k2" "PossionSeq_RUVse_k3"
校正后的數(shù)據(jù)可以通過 Counts
函數(shù)獲取糜烹。
head(enONE::Counts(Enone, slot="sample", method="DESeq_RUVg_k2"))[,1:5]
#> D1 D2 D3 D4 D5
#> ENSG00000227232 8.513936 11.645203 12.431441 7.453572 11.104087
#> ENSG00000268903 3.305409 8.100305 5.277670 26.325845 7.088225
#> ENSG00000269981 2.270299 4.876642 2.402698 18.118129 6.068101
#> ENSG00000279457 11.895695 9.998730 8.760526 11.204313 16.375896
#> ENSG00000225630 13.767575 16.292915 15.365490 10.140265 17.834578
#> ENSG00000237973 28.042040 37.277287 39.847772 34.078301 33.952026
數(shù)據(jù)校正效果評估
enONE
利用了八個與基因表達分布的不同方面相關(guān)的標準化性能指標來評估數(shù)據(jù)標準化的性能违诗。這八個指標分別為:
-
BIO_SIM
:生物組的相似性。通過計算前三個表達PCs(默認情況下)上的歐氏距離度量疮蹦,由condition
列定義的聚類的平均輪廓寬度诸迟。BIO_SIM
越大越好。 -
EN_SIM
:富集組的相似性愕乎。通過計算前三個表達PCs(默認情況下)上的歐氏距離度量阵苇,由enrich
列定義的聚類的平均輪廓寬度。EN_SIM
越大越好感论。 -
BAT_SIM
:批次組的相似性绅项。通過計算前三個表達PCs(默認情況下)上的歐氏距離度量,由batch
列定義的聚類的平均輪廓寬度比肄。BAT_SIM
越小越好快耿。 -
PAM_SIM
:PAM聚類組的相似性。通過計算前三個表達PCs(默認情況下)上的歐氏距離度量芳绩,由PAM 聚類(根據(jù)eval.pam.k
進行聚類)定義的聚類的最大平均輪廓寬度掀亥。PAM_SIM
越大越好。 -
WV_COR
:研究相關(guān)變量的保留程度妥色。對原始計數(shù)的正評估基因(PosEvaluation
)子矩陣的前eval.pc.n
個PCs 進行回歸的 R2 測量搪花。WV_COR
越大越好。 -
UV_COR
:無關(guān)誤差移除的程度嘹害。對原始計數(shù)的負評估基因(NegEvaluation
)子矩陣的前eval.pc.n
個PCs 進行回歸的 R2 測量撮竿。UV_COR
越小越好。 -
RLE_MED
:平均平方中位數(shù)相對對數(shù)表達(RLE)吼拥。RLE_MED
越小越好倚聚。 -
RLE_IQR
:相對于四分位距(IQR)的方差。RLE_IQR
越小越好凿可。
評估指標可以通過 getMetrics
進行獲取惑折。
getMetrics(Enone)[1:5,]
#> BIO_SIM EN_SIM BATCH_SIM PAM_SIM RLE_MED
#> DESeq_RUVg_k2 -0.02467593 0.3576650 -0.03483061 0.5713398 9.319387e-06
#> PossionSeq_RUVs_k3 0.10240479 0.3374029 -0.03965069 0.5556398 9.971803e-05
#> TMM_RUVg_k2 -0.02664678 0.3518415 -0.03524885 0.5595637 8.500478e-05
#> DESeq_RUVg_k3 -0.04651463 0.3485155 -0.04427468 0.5708689 1.370353e-05
#> DESeq_RUVs_k1 0.03352455 0.3391415 -0.04444996 0.5224180 7.861089e-06
#> RLE_IQR WV_COR UV_COR
#> DESeq_RUVg_k2 0.011996148 0.7401496 0.3258051
#> PossionSeq_RUVs_k3 0.012314397 0.7709943 0.3154878
#> TMM_RUVg_k2 0.012089861 0.7404610 0.3288269
#> DESeq_RUVg_k3 0.011133971 0.6688187 0.4166522
#> DESeq_RUVs_k1 0.009145155 0.6495402 0.3680677
評估指標通過轉(zhuǎn)換為秩值授账,并通過平均各指標的秩值來獲得特定方法的打分(score)。數(shù)據(jù)標準化效果的打分可以通過getScore
獲取惨驶。
getScore(Enone)[1:5,]
#> BIO_SIM EN_SIM BATCH_SIM PAM_SIM RLE_MED RLE_IQR WV_COR
#> DESeq_RUVg_k2 32 59 40 59 56 25 20
#> PossionSeq_RUVs_k3 60 23 47 50 33 23 23
#> TMM_RUVg_k2 30 50 42 54 39 24 21
#> DESeq_RUVg_k3 20 48 53 58 54 29 12
#> DESeq_RUVs_k1 55 32 54 25 57 32 5
#> UV_COR SCORE
#> DESeq_RUVg_k2 59 43.750
#> PossionSeq_RUVs_k3 60 39.875
#> TMM_RUVg_k2 58 39.750
#> DESeq_RUVg_k3 39 39.125
#> DESeq_RUVs_k1 52 39.000
選擇數(shù)據(jù)標準化方法
enONE
提供了biplot來探索數(shù)據(jù)校正方法的全空間白热,可以通過interactive=TRUE
來開啟交互模式。
# get performance score
enScore <- getScore(Enone)
# perform PCA based on evaluation score, excluding BAT_SIM column (3) if no batch information provided, and SCORE column (9).
# pca.eval <- prcomp(enScore[,-c(3, 9)], scale = TRUE)
pca.eval <- prcomp(enScore[,-c(9)], scale = TRUE)
# pca biplot
PCA_Biplot(pca.eval, score = enScore$SCORE, interactive = FALSE)
在這個圖中粗卜,每個點對應一個標準化程序屋确,并顏色代表性能得分(八個性能指標rank的均值)。藍色箭頭對應于性能指標的PCA載荷续扔。藍色箭頭的方向和長度可以解釋為每個指標對前兩個主成分的貢獻程度攻臀。
這里使用打分最高的標準化程序 DESeq_RUVg_k2
進行下游分析。
# select normalization
norm.method <- rownames(enScore[1,])
# get normalized counts
norm.data <- enONE::Counts(Enone, slot = "sample", method = norm.method)
# get normalization factors
norm.factors <- getFactor(Enone, slot = "sample", method = norm.method)
norm.method
#> [1] "DESeq_RUVg_k2"
需要注意的是纱昧,如果
enONE
運行的時候沒有返回校正后數(shù)據(jù)(i.e.,return.norm=FALSE
)刨啸,可以通過下面的代碼手動執(zhí)行數(shù)據(jù)標準化。# perform normalization Enone <- UseNormalization(Enone, slot = "sample", method = norm.method) # get normalized counts norm.data <- Counts(Enone, slot = "sample", method = norm.method)
數(shù)據(jù)校正的效果
這里使用PCA展示數(shù)據(jù)校正的效果识脆。
# create sample name, e.g., Y3.Input
samples_name <- paste(Enone$condition, Enone$replicate, sep=".")
# PCA for raw count
p1 <- PCAplot(enONE::Counts(Enone, slot="sample", "Raw"),
color = Enone$batch,
shape = Enone$enrich,
label = samples_name,
vst.norm = TRUE) +
ggtitle("Before normalization")
# PCA for normalized count
p2 <- PCAplot(log1p(norm.data),
color = Enone$batch,
shape = Enone$enrich,
label = samples_name,
vst.norm = FALSE) +
ggtitle("After normalization")
# combine two plots
p1 + p2
鑒定富集基因
enONE
可用于富集基因的鑒定设联。FindEnrichment
函數(shù)自動地根據(jù)condition
列中的每個生物學分組鑒定富集的基因。
默認情況下灼捂,富集基因(即NAD-RNA)為Enrichment/Input ≥ 2(logfc.cutoff = 1
)离例,F(xiàn)DR < 0.05(p.cutoff = 0.05
)的基因。
使用getEnrichment
函數(shù)可以獲取一個列表包含了各組的富集結(jié)果悉稠。
# find all enriched genes
Enone <- FindEnrichment(Enone, slot="sample", norm.method = norm.method,
logfc.cutoff = 1, p.cutoff = 0.05)
#> - Estimate dispersion & Fit GLM...
#> - Testing differential genes...
# get filtered enrichment results
res.sig.ls <- getEnrichment(Enone, slot="sample", filter=TRUE)
# count number of enrichment in each group
unlist(lapply(res.sig.ls, nrow))
#> Y.Enrich_Y.Input M.Enrich_M.Input O.Enrich_O.Input
#> 579 672 683
res.sig.ls
中每個表都是一個data.frame
宫蛆,其中行為基因,列為基因的相關(guān)信息(GeneID偎球、logFC洒扎、p值等)。表格中包含以下列:
GeneID
:基因的ID衰絮。logFC
:富集樣本和輸入樣本之間的log2倍變化袍冷。正值表示基因在富集組中更高度富集。logCPM
:所有樣本平均表達的log2 CPM(counts per million)猫牡。LR
:似然比檢驗的似然比胡诗。PValue
:來自似然比檢驗的p值。FDR
:p值的假陽性發(fā)現(xiàn)率淌友,默認使用"BH"方法煌恢。
head(res.sig.ls[[1]])
#> GeneID logFC logCPM LR PValue FDR
#> 1 ENSG00000013275 2.271643 7.518198 731.6422 3.936608e-161 5.444329e-157
#> 2 ENSG00000113387 2.643018 9.732157 631.6712 2.164465e-139 1.496728e-135
#> 3 ENSG00000104853 2.244442 7.839256 587.7484 7.738626e-130 3.567507e-126
#> 4 ENSG00000105185 2.094086 6.147852 558.6392 1.662061e-123 5.746574e-120
#> 5 ENSG00000143110 3.084209 9.115976 550.0416 1.232914e-121 3.410239e-118
#> 6 ENSG00000173915 2.324782 7.005478 483.7653 3.239173e-107 7.466295e-104
enONE
提供 reduceRes
函數(shù)將各個富集表格的結(jié)果轉(zhuǎn)換到同一個data.frame
中。接著震庭,可以使用 BetweenStatPlot
可視化各組NAD-RNA修飾水平瑰抵。
# simplify group id
names(res.sig.ls) <- c("Young", "Mid", "Old")
# logfc.col specify the name of logFC column
nad_df1 <- reduceRes(res.sig.ls, logfc.col = "logFC")
# convert the Group column as factor
nad_df1$Group <- factor(nad_df1$Group, levels = unique(nad_df1$Group))
# draw plot
bxp1 <- BetweenStatPlot(nad_df1, x="Group", y="logFC", color="Group",
step.increase = 0.6, add.p = "p",
comparisons = list(c("Young", "Mid"),
c("Young", "Old")))
bxp1
處理synthetic spike-in
由于樣本中加入了synthetic RNA茎匠,其中一個帶有5%的NAD帽局荚,另一個帶有100%的m7G帽箫津,我們可以利用這些spike-ins來確定捕獲的靈敏度和特異性俭嘁。
synEnrichment
計算具有給定標準化方法的spike-ins的富集水平。DotPlot
可用于可視化spike-in的富集水平肴颊。
# compute synthetic spike-in enrichment
syn_level <- synEnrichment(Enone, method=norm.method, log=TRUE)
# transform to long format
syn_df <- as.data.frame(syn_level) %>%
rownames_to_column("syn_id") %>%
pivot_longer(cols = -syn_id,
names_to = "id",
values_to = "logFC") %>%
left_join(meta[,c("id","condition")], by="id")
# remove suffix of condition for simplification
syn_df$condition <- gsub("\\..*", "", syn_df$condition)
# rename facet label
samples_label <- setNames(c("5% NAD-RNA", "100% m7G-RNA"),
nm=c("Syn1", "Syn2"))
# draw dotplot
DotPlot(syn_df, x="syn_id", y="logFC", fill="syn_id") +
theme(legend.position = "none") +
scale_x_discrete(labels=samples_label)
#> Bin width defaults to 1/30 of the range of the data. Pick better value with `binwidth`.
其中氓栈,具有5% NAD-RNA的Syn1有明顯富集,而100% m7G-RNA的Syn2并沒有出現(xiàn)富集婿着,表明了富集實驗的特異性授瘦。
# save Enone data
save(Enone, file="data/Enone.RData")
Session Info
sessionInfo()
#> R version 4.2.1 (2022-06-23 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22631)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=Chinese (Simplified)_China.utf8
#> [2] LC_CTYPE=Chinese (Simplified)_China.utf8
#> [3] LC_MONETARY=Chinese (Simplified)_China.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=Chinese (Simplified)_China.utf8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] patchwork_1.1.2 forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10
#> [5] purrr_0.3.5 readr_2.1.3 tidyr_1.2.1 tibble_3.1.8
#> [9] ggplot2_3.3.6 tidyverse_1.3.2 enONE_1.0.0
#>
#> loaded via a namespace (and not attached):
#> [1] googledrive_2.0.0 colorspace_2.0-3
#> [3] ggsignif_0.6.4 ellipsis_0.3.2
#> [5] class_7.3-20 modeltools_0.2-23
#> [7] EnvStats_2.7.0 mclust_5.4.10
#> [9] XVector_0.36.0 fs_1.5.2
#> [11] GenomicRanges_1.48.0 rstudioapi_0.14
#> [13] farver_2.1.1 ggpubr_0.4.0
#> [15] ggrepel_0.9.1 flexmix_2.3-18
#> [17] bit64_4.0.5 lubridate_1.8.0
#> [19] AnnotationDbi_1.58.0 fansi_1.0.3
#> [21] xml2_1.3.3 codetools_0.2-18
#> [23] splines_4.2.1 cachem_1.0.6
#> [25] robustbase_0.95-0 geneplotter_1.74.0
#> [27] knitr_1.40 jsonlite_1.8.2
#> [29] broom_1.0.1 annotate_1.74.0
#> [31] cluster_2.1.3 kernlab_0.9-31
#> [33] dbplyr_2.2.1 png_0.1-7
#> [35] compiler_4.2.1 httr_1.4.4
#> [37] backports_1.4.1 assertthat_0.2.1
#> [39] Matrix_1.5-1 fastmap_1.1.0
#> [41] lazyeval_0.2.2 gargle_1.2.1
#> [43] limma_3.52.4 cli_3.4.1
#> [45] htmltools_0.5.3 tools_4.2.1
#> [47] gtable_0.3.1 glue_1.6.2
#> [49] GenomeInfoDbData_1.2.8 Rcpp_1.0.9
#> [51] carData_3.0-5 Biobase_2.56.0
#> [53] cellranger_1.1.0 vctrs_0.5.1
#> [55] Biostrings_2.64.1 fpc_2.2-9
#> [57] xfun_0.33 rvest_1.0.3
#> [59] lifecycle_1.0.3 googlesheets4_1.0.1
#> [61] rstatix_0.7.0 XML_3.99-0.11
#> [63] edgeR_3.38.4 DEoptimR_1.0-11
#> [65] zlibbioc_1.42.0 MASS_7.3-57
#> [67] scales_1.2.1 hms_1.1.2
#> [69] MatrixGenerics_1.8.1 parallel_4.2.1
#> [71] SummarizedExperiment_1.26.1 RColorBrewer_1.1-3
#> [73] yaml_2.3.5 memoise_2.0.1
#> [75] pbapply_1.5-0 stringi_1.7.8
#> [77] RSQLite_2.2.18 highr_0.9
#> [79] genefilter_1.78.0 S4Vectors_0.34.0
#> [81] BiocGenerics_0.42.0 BiocParallel_1.30.4
#> [83] GenomeInfoDb_1.32.4 rlang_1.0.6
#> [85] pkgconfig_2.0.3 prabclus_2.3-2
#> [87] matrixStats_0.62.0 bitops_1.0-7
#> [89] evaluate_0.17 lattice_0.20-45
#> [91] labeling_0.4.2 htmlwidgets_1.5.4
#> [93] bit_4.0.4 tidyselect_1.2.0
#> [95] magrittr_2.0.3 DESeq2_1.36.0
#> [97] R6_2.5.1 IRanges_2.30.1
#> [99] generics_0.1.3 DelayedArray_0.22.0
#> [101] DBI_1.1.3 withr_2.5.0
#> [103] haven_2.5.1 pillar_1.8.1
#> [105] survival_3.3-1 KEGGREST_1.36.3
#> [107] abind_1.4-5 RCurl_1.98-1.9
#> [109] nnet_7.3-17 modelr_0.1.9
#> [111] crayon_1.5.2 car_3.1-1
#> [113] utf8_1.2.2 plotly_4.10.0
#> [115] tzdb_0.3.0 rmarkdown_2.17
#> [117] readxl_1.4.1 locfit_1.5-9.6
#> [119] grid_4.2.1 data.table_1.14.2
#> [121] blob_1.2.3 reprex_2.0.2
#> [123] digest_0.6.29 diptest_0.76-0
#> [125] xtable_1.8-4 paintingr_0.1.1
#> [127] stats4_4.2.1 munsell_0.5.0
#> [129] viridisLite_0.4.1