内容

1简介

XBSeq是一种检测RNA-seq差分表达(DE)的新算法，该算法基于观测到的信号是真实表达信号和测序噪声的卷积假设，建立了一个统计模型。在非外显子区域映射的读取被认为是序列噪声，它遵循泊松分布。给定可测量的观察信号和RNA-seq数据的背景噪声，可以描绘出真实的表达信号，假设受负二项分布支配，从而准确检测差异表达基因。XBSeq论文发表于BMC genomics[1]。我们最近还整合了用于测试差异替代多聚腺苷酸(apa)使用的roar包的功能。

2安装

XBSeq可以从Bioconductor安装

如果(!requireNamespace("BiocManager"， quiet =TRUE)) install.packages("BiocManager")::install("XBSeq")
库(“XBSeq”)

您还可以通过安装XBSeq的开发版本

库(devtools) install_github(“liuy12 / XBSeq”)

3.使用XBSeq测试差分表达式

为了使用XBSeq测试DE，序列对齐后，我们需要进行两次读数计数程序，测量映射到外显子区域(观测信号)和非外显子区域(背景噪声)的读数。有几种现有的方法可以达到这个目的。这里我们介绍了使用featu复述、HTSeq和summarizeOverlaps的读计数

3．1featureCounts

featurets是一个读取摘要程序，可用于从RNA或DNA测序技术产生的读取，它实现了高效的染色体哈希和特征阻塞技术，速度相当快，需要较少的计算机内存。featurerts可以从r中的Rsubread包中获得。基本上，你需要运行以下命令:

fc_signal <- featucounts (files = bamLists, annot.)ext = gtf_file, isGTFAnnotationFile = TRUE) fc_bg <- featurecrets (files = bamLists, annot. txt)ext = gtf_file_bg, isGTFAnnotationFile = TRUE)

用于测量观测信号和背景噪声的gtf文件可以在github的gtf文件夹中下载:https://github.com/Liuy12/XBSeq_files．我们已经为几种不同基因组的生物构建了注释。如果您想自己构造gtf文件，我们还在github中存放了用于构造gtf文件的perl脚本。关于我们用来构造背景gtf文件的过程的详细信息，可以在小插图的详细信息部分中找到。

3.2HTSeq过程

或者，你也可以使用HTSeq进行读取计数:

htseq-count [options]   > Observed_count.txt htseq-count [options]   > background_count.txt

关于HTSeq如何工作的详细信息可以在这里找到:http://www-huber.embl.de/HTSeq/doc/count.html

3．3summarizeOverlaps

你也可以使用summarizeOverlaps该功能可从GenomicRanges包获得:

features_signal <- import(gtf_file) features_signal <- split(features_signal, mcols(features_signal)$gene_id) so_signal <- summarizeOverlaps(features_signal, bamLists) ## for背景噪声features_bg <- import(gtf_file_bg) features_bg <- split(features_bg, mcols(features_bg)$gene_id) so_bg <- summarizeOverlaps(features_bg, bamLists)

3.4测试APA使用差异

选择性聚腺苷酸化(APA)是一种广泛存在的机制，其中一个基因使用替代poly(a)位点来编码多个不同3'UTR长度的mRNA转录本。用户可以直接从RNA-seq比对文件中推断不同的APA使用:

apaStats <- apaUsage(bamTreatment, bamControl, apaAnno)

在哪里bamTreatment是否有治疗条件和数据的bam对齐文件名的完整路径列表bamControl是与要考虑的控制条件的数据进行bam对齐的文件名的完整路径的列表。apaAnno吼叫包使用的apa注释的完整路径。APA注释的几种生物不同的基因组构建可以从下载在这里．有关如何构建APA注释的详细信息，请参阅详细信息部分。

3．5XBSeq测试DE

在HTSeq程序之后，我们将对每个基因进行两次测量，观测到的信号和背景噪声。在这里，我们将使用小鼠RNA-seq数据集，其中包含3个重复的野生型小鼠肝组织(WT)和3个重复的Myc转基因小鼠肝组织(Myc)。数据集来源于Gene Expression Omnibus(GSE61875)．

作为初步步骤，我们已经进行了上述HTSeq程序，为每个基因生成观测信号和背景噪声。这两个数据集可以被加载到用户的工作空间

数据(ExampleData)

我们可以先看看这两个数据集:

(观察)
样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本样本
(背景)
0610005C13Rik 512 374 466 496 ## 0610007C21Rik 50 44 40 42 ## 0610007L01Rik 33 22 35 39 ## 0610007P08Rik 21 17 22 10 ## 0610007P22Rik 16 19 26 14 # 0610007L01Rik 504 648 ## 0610007C21Rik 42 60 ## 0610007L01Rik 40 22 ## 0610007P14Rik 18 ## 0610007P14Rik 24 23 ## 0610007P22Rik 28 28

行表示映射到每个基因或相应背景区域的读值。列表示样本。微分表达式分析如下:

首先，我们需要构造一个XBSeqDataSet对象。条件是实验的设计矩阵。观察和背景是来自HTSeq的输出矩阵(记住要删除输出矩阵的汇总统计的底部几行)。

条件<- factor(c(rep(" c "， 3)， rep("T"， 3))) XB <- XBSeqDataSet(已观察到，背景，条件)

通常建议您事先检查观察到的信号和背景噪声的分布。我们提供函数' r XBplot '来实现这一点。我们建议通过将参数unit设置为“logTPM”来检查以log2 ppm (TPM)为单位的分布。基因长度信息通过“ExampleData”加载。理想情况下，还应该提供库的大小。默认情况下，使用映射到外显子区域的所有读的和。

XBplot(XB, Samplenum = 1, unit = "LogTPM"， genelth = genelth [， 2])
XBplot(XB, Samplenum = 1, unit = "LogTPM"， Genelength = Genelength [，: Libsize未提供，每个样本中映射到exonic ##区域的所有读计数的总和被用作该样本的总库大小
##警告:删除了16453行包含非有限值(stat_bin)。
##警告:删除4行包含丢失的值(geom_bar)。

然后估计初步的潜在信号，然后是归一化因子和色散估计

XB <- estimateRealCount(XB) XB <- estimateSizeFactors(XB) XB <- estimateSCV(XB, method = "pooled"， sharingMode = "maximum"， fitType = "local")

看看scv的合身信息

plotSCVEsts (XB)

使用函数XBSeqTest执行DE测试

Teststas <- XBSeqTest(XB, levels(conditions)[1L]， levels(conditions)[2L]， method ='NP')

根据测试统计数据绘制地图

mapplot (Teststas, padj = FALSE, pcuff = 0.01, lfccuff = 1)

#或者，上面所有的代码都可以用包装器函数XBSeq Teststats <- XBSeq(观察到的，背景，条件，方法= "pooled"， sharingMode = "maximum"， fitType = "local"， pvals_only = FALSE, parameterthod = "NP")

3.6将结果与DESeq进行比较

现在我们将使用DESeq对同一数据集进行DE分析，然后比较这两种方法得到的结果

如果您以前没有安装过DESeq，也可以从Bioconductor获得DESeq

BiocManager::安装(“DESeq”)

那么对DESeq进行DE分析可以通过:

库(“DESeq”)
##警告:软件包“DESeq”已弃用，将从Bioconductor ## 3.13版本中删除。请使用DESeq2
library('ggplot2') de <- newCountDataSet(Observed, conditions) de <- estimateSizeFactors(de) de <- estimateDispersions(de, method =" pooled"， fitType="local") res <- nbinomTest(de, levels(conditions)[1]， levels(conditions)[2])

然后我们可以比较XBSeq和DESeq的结果

DE_index_DESeq <——(res (pval < 0.01 & abs (log2FoldChange) > 1)) DE_index_XBSeq <——(Teststas, (pval < 0.01 & abs (log2FoldChange) > 1)) DE_index_inters < -相交(DE_index_DESeq DE_index_XBSeq) DE_index_DESeq_uniq < - setdiff (DE_index_DESeq DE_index_XBSeq) DE_plot < - MAplot (Teststas padj = FALSE, pcuff = 0.01, lfccuff = 1,形状= 16)DE_plot + geom_point (data = Teststas DE_index_inters,, aes (x = baseMean, y = log2FoldChange),颜色=“绿色”,+ geom_point(data = Teststas[DE_index_DESeq_uniq，]， aes(x = baseMean, y = log2FoldChange)， color = "blue"， shape = 16)

红点表示仅通过XBSeq鉴定的DE基因。绿点是XBSeq和DESeq的共享结果。蓝点是仅由DESeq识别的DE基因。

4细节

4．1背景区域gtf文件的构建

  • 外显子区域注释来自UCSC数据库。

  • 非外显子区域的构造遵循以下几个标准:
    1. 从UCSC数据库下载refFlat表
    2. 基因组的几个功能元件(mRNA、peudo基因等)被排除在全基因组注释之外。
    3. 构建每个基因的背景区域，使该区域与基因的外显子区域具有相同的结构或长度。

关于如何构造真实示例的背景区域注释文件的更多详细信息，可以在ExampleData的手册页和我们的XBSeq出版物中找到。

4.2构建咆哮包使用的APA注释的gtf文件

APA位点通过使用POLYAR程序[2]进行预测，该程序通过使用12种不同的先前映射的聚(A)信号(PAS)六聚体应用期望最大化(EM)方法。POLYAR预测的APA位点可分为pa -强、pa -中、pa -弱三类。只选择pa -strong类中的APA位点构建最终的APA注释。已经构建了不同版本的人类和小鼠基因组APA注释，并可从github下载:https://github.com/Liuy12/XBSeq_files

4.３关于内含子保留

我经常被问到内含子保留事件是否对XBSeq的性能有任何影响。内含子保留是控制转录组活性的一种常见剪接机制。一些文章已经证明，内含子保留的转录本会通过一种称为无意义介导的mRNA衰变(nonnonsense -mediated mRNA decay, NMD)的机制被降解[3-5]。在我看来，内含子保留不会影响XBSeq的性能，因为这种类型的转录本最终会被降解，将它们视为背景噪声是有意义的。

5XBSeq的最新更新:XBSeq2

6错误报告

将错误作为问题报告在我们的GitHub库或者你也可以直接向我的邮箱举报:liuy12@uthscsa.edu．

7会话信息

sessionInfo ()
## R版本4.0.3(2020-10-10)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 18.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas。所以## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_US。UTF-8 LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:##[1]并行stats4统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:[7] SummarizedExperiment_1.20.0 Biobase_2.50.0 ## [9] MatrixGenerics_1.2.0 matrixStats_0.57.0 ## [13] IRanges_2.24.0 S4Vectors_0.28.0 ## [15] BiocGenerics_0.36.0 BiocStyle_2.18.0 ## ##通过命名空间加载(并且没有附加):## [7] GenomeInfoDbData_1.2.4 Rsamtools_2.6.0 yaml_2.2.1 ## [10] pillar_1.4.6 RSQLite_2.2.1 glue_1.4.2 ## [13] digest_0.6.27 RColorBrewer_1.1-2 XVector_0.30.0 ## [19] XML_3.99-0.5 pkgconfig_2.0.3 magick_2.5.0 ## [22] genefilter_1.72.0 bookdown_0.21 zlibbioc_1.36.0 ## [25] purrr_0.3.4 xtable_1. 4 scales_1.1.1 ## [28] pracma_2.2.9 BiocParallel_1.24.0 ##tibble_3.0.4 ## [31] annotate_1.68.0 farver_2.0.3 generics_0.0.2 ## [34] ellipsis_0.3.1 withr_2.3.0 survival_3.2-7 ## [37] magrittr_1.5 crayon_1.3.4 memoise_1.1.0 ## [40] evaluate_0.14 tools_4.0.3 formatR_1.7 ## [43] lifecycle_0.2.0 stringr_1.4.0 munsell_0.5.0 ## [46] DelayedArray_0.16.0 AnnotationDbi_1.52.0 Biostrings_2.58.0 ## [49] compiler_4.0.3 rlang_0.4.8 grid_4.0.3 ## [52] RCurl_1.98-1.2 labeling_0.4.2 bitops_1.0-6 ## [55] rmarkdown_2.5 gtable_0.3.0 DBI_1.1.0 ## [58] R6_2.4.1 GenomicAlignments_1.26.0 knitr_1.30 ## [61] dplyr_1.0.2 rtracklayer_1.50.0 bit_4.0.4 ## [64] stringi_1.5.3 Rcpp_1.0.5 vctrs_0.3.4 ## [67] geneplotter_1.68.0 tidyselect_1.1.0 xfun_0.18

8确认

XBSeq是基于DESeq和DESeq2的源代码在R语言中实现的。

9参考文献

  1. 陈海义，刘颖，邹颖，赖志强，D. Sarkar, Huang Y.等，“基于非外显子映射reads的RNA测序数据差异表达分析”，BMC Genomics, vol. 16增刊7,p. S14, 2015年6月11日。
  2. Akhtar MN, Bukhari SA, Fazal Z, Qamar R, Shahmuradov IA: POLYAR，一种用于预测人类序列中poly(a)位点的新计算机程序。BMC基因组学2010,11:646。
  3. Jung, Hyunchul等人，内含子保留是肿瘤抑制因子失活的普遍机制。自然遗传学(2015)。
  4. brunschweig, Ulrich等，“哺乳动物中广泛的内含子保留在功能上调节转录组。”基因组研究24.11(2014):1774-1786。
  5. Lykke-Andersen, Søren和Torben Heick Jensen。“无意义介导的mRNA衰变:形成转录组的复杂机制。”自然评论分子细胞生物学(2015)。