enONE

最近開發(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:

  1. 來自模式生物黑腹果蠅的總RNA欺冀,用于 enONE 的normalization树绩;
  2. Synthetic RNA,其中含有5%的NAD修飾形式(相對于m7G-RNA)隐轩,用于確定富集靈敏度饺饭;
  3. 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 流程包括以下四步:

  1. 質(zhì)量控制(Quality control);
  2. 基因集選让怠(Gene set selection);
  3. 數(shù)據(jù)校正(Normalization procedures);
  4. 數(shù)據(jù)校正效果評估(Normalization performance assessment).

質(zhì)量控制(Quality control)

enONE 首先進行質(zhì)量控制步驟诽俯,通過 FilterLowExprGene 保留至少在 n 個樣本中至少具有 min.count 的基因來完成此步驟。ngroup 中指定的最小樣本組的大小確定承粤。

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_metricsenone_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 整合了五種方法眉枕,包括:

  1. Total Count (TC);

  2. Upper-Quartile (UQ);

  3. Trimmed Mean of M Values (TMM);

  4. DESeq;

  5. PossionSeq.

默認情況下恶复,enONE 使用所有縮放方法怜森,但用戶可以在 scaling.method 參數(shù)中選擇要使用的縮放程序。

對于regression-based方法谤牡,enONE 利用了三種RUV的方法副硅,包括:

  1. RUVg;

  2. RUVs;

  3. 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
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市竟宋,隨后出現(xiàn)的幾起案子提完,更是在濱河造成了極大的恐慌,老刑警劉巖袜硫,帶你破解...
    沈念sama閱讀 206,013評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件氯葬,死亡現(xiàn)場離奇詭異挡篓,居然都是意外死亡婉陷,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2 382
  • 文/潘曉璐 我一進店門官研,熙熙樓的掌柜王于貴愁眉苦臉地迎上來秽澳,“玉大人,你說我怎么就攤上這事戏羽〉I瘢” “怎么了?”我有些...
    開封第一講書人閱讀 152,370評論 0 342
  • 文/不壞的土叔 我叫張陵始花,是天一觀的道長妄讯。 經(jīng)常有香客問我,道長酷宵,這世上最難降的妖魔是什么亥贸? 我笑而不...
    開封第一講書人閱讀 55,168評論 1 278
  • 正文 為了忘掉前任,我火速辦了婚禮浇垦,結(jié)果婚禮上炕置,老公的妹妹穿的比我還像新娘。我一直安慰自己男韧,他們只是感情好朴摊,可當我...
    茶點故事閱讀 64,153評論 5 371
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著此虑,像睡著了一般甚纲。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上朦前,一...
    開封第一講書人閱讀 48,954評論 1 283
  • 那天介杆,我揣著相機與錄音讹弯,去河邊找鬼。 笑死这溅,一個胖子當著我的面吹牛组民,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播悲靴,決...
    沈念sama閱讀 38,271評論 3 399
  • 文/蒼蘭香墨 我猛地睜開眼臭胜,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了癞尚?” 一聲冷哼從身側(cè)響起耸三,我...
    開封第一講書人閱讀 36,916評論 0 259
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎浇揩,沒想到半個月后仪壮,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 43,382評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡胳徽,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2 323
  • 正文 我和宋清朗相戀三年积锅,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片养盗。...
    茶點故事閱讀 37,989評論 1 333
  • 序言:一個原本活蹦亂跳的男人離奇死亡缚陷,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出往核,到底是詐尸還是另有隱情箫爷,我是刑警寧澤,帶...
    沈念sama閱讀 33,624評論 4 322
  • 正文 年R本政府宣布聂儒,位于F島的核電站虎锚,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏衩婚。R本人自食惡果不足惜窜护,卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望谅猾。 院中可真熱鬧柄慰,春花似錦、人聲如沸税娜。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽敬矩。三九已至概行,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間弧岳,已是汗流浹背凳忙。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1 260
  • 我被黑心中介騙來泰國打工业踏, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留,地道東北人涧卵。 一個月前我還...
    沈念sama閱讀 45,401評論 2 352
  • 正文 我出身青樓勤家,卻偏偏與公主長得像,于是被迫代替她去往敵國和親柳恐。 傳聞我的和親對象是個殘疾皇子伐脖,可洞房花燭夜當晚...
    茶點故事閱讀 42,700評論 2 345

推薦閱讀更多精彩內(nèi)容