NBAMSeq: RNA-Seq数据的负二项式加法模型
徐仁，裴芬宽
2022-11-01
安装
安装和加载NBAMSeq
简介
高通量测序实验和差异表达分析是一种广泛应用的检测基因组生物标志物的方法。差异表达分析的一个基本步骤是建立基因计数和相关协变量之间的关联模型。NBAMSeq是一种基于广义可加性模型的灵活统计模型，允许在方差估计中跨基因共享信息。具体来说，我们将平均基因计数的对数建模为平滑函数的和，同时通过嵌套迭代估计平滑参数和系数。方差由贝叶斯收缩方法估计，以充分利用所有基因的信息。
NBAMSeq的工作流程主要包括三个步骤:
步骤1:数据输入使用
NBAMSeqDataSet；步骤2:差分表达式(DE)分析
NBAMSeq函数;步骤3:使用提取DE结果
结果函数。
下面我们分别说明这些步骤。
数据输入
用户需要提供三个部分的输入，即countData，colData,设计．
countData是RNASeq实验生成的基因计数矩阵。
countData的一个例子n =50## n代表基因数量m =20.## m代表样本量countData =矩阵（rnbinom(n＊米,μ=One hundred.，大小=1/3.)，ncol =米)+1模式(countData) =“整数”colnames(countData) =paste0（“样本”，1：米)rownames(countData) =paste0（“基因”，1：n)头(countData)sample1 sample2 sample3 sample4 sample5 sample6 sample7 sample8 sample9 gene1 3 5 2 2 61 13 212 6 107 gene2 1 106 1 22 1 335 2 36 gene3 7 98 18 1 30 64 41 35 17 gene4 41 6 1 60 143 319 11 4 gene5 10 3 1 1 312年43 3 22 gene6 1 81 449 74 362 4 19 532 sample10 sample11 sample12 sample13 sample14 sample15 sample16 sample17 gene1 43 273 71 30 14 1 1 gene2 5 1 120 6 16 2 7 805 gene3 55 120 6 59 1 11 14 572 gene4 31 344 58 12 20 330 gene5 26 109 5 267 26 76 200 166 gene6 35 69 59 243188 21 632 sample18 sample19 sample20 gene1 155 289 799 gene2 64 95 1 gene3 308 313 10 gene4 27 7 4 gene5 139 340 12 gene6 11 95 4colData是一个包含样本协变量的数据帧。样品订单colData样品订单要匹配吗countData．
colData的一个例子把=runif(m,20.，80）var1 =rnorm(m)var2 =rnorm(m)var3 =rnorm(m)var4 =as.factor（样本（c（0，1，2), m,取代=真正的）)colData =data.frame（把=把,var1 =var1,var2 =var2,var3 =var3,var4 =var4)rownames(colData) =paste0（“样本”，1：米)头(colData)Pheno var1 var2 var3 var4 sample1 36.45066 -0.02413811 1.2940361 0.57054527 0 sample2 58.61363 0.05496400 0.8194974 0.81123820 1 sample3 63.32145 -0.00449474 0.3444390 0.69527733 1 sample4 60.62626 0.26733796 0.1804601 0.04029040 2 sample5 44.44742 0.77958638 0.5657089 0.09589169 1 sample6 75.39852 0.54212324 0.8143508 0.97449953 0设计是指定如何对样本建模的公式。与其他执行DE分析的包(包括DESeq2)相比(Love, Huber, and Anders 2014),磨边机(罗宾逊，麦卡锡，和史密斯2010), NBPSeq(Di et al. 2015)和BBSeq(Zhou, Xia, and Wright 2011)， NBAMSeq通过mgcv支持协变量的非线性模型(《木材与木材》2015)．为了表示模型中的非线性协变量，用户需要使用s (variable_name)在设计公式。在我们的例子中，如果我们想建模把作为一个非线性协变量，设计公式应为:
关于。应该做几个注意事项设计公式:
支持多个非线性协变量，例如:
设计= ~ s(pheno) + s(var1) + var2 + var3 + var4；非线性协变量不能是离散变量，例如。
设计= ~ s(pheno) + var1 + var2 + var3 + s(var4)作为var4是一个因素，将一个因素建模为非线性是没有意义的;中至少应提供一个非线性协变量
设计．如果假设所有协变量对基因计数有线性影响，则使用DESeq2(Love, Huber, and Anders 2014),磨边机(罗宾逊，麦卡锡，和史密斯2010), NBPSeq(Di et al. 2015)或BBSeq(Zhou, Xia, and Wright 2011)代替。如。设计= ~ pheno + var1 + var2 + var3 + var4NBAMSeq不支持;不支持设计矩阵。
然后我们构建NBAMSeqDataSet使用countData，colData,设计：
类:NBAMSeqDataSet dim: 50 20元数据(1):拟合的assays(1):计数rownames(50): gene1 gene2…gene49 gene50 rowData names(0): colnames(20): sample1 sample2…sample19 sample20 colData names(5): pheno var1 var2 var3 var4微分表达式分析
微分表达式分析可由NBAMSeq功能:
其他几个论点NBAMSeq用户可自定义分析功能。
γ参数可用于控制非线性函数的平滑度。更高的γ意味着非线性函数将更加平滑。看到γ的观点gam在MGCV中的作用(《木材与木材》2015)获取详细信息。默认的γ是2.5;fitlin要么是真正的或假说明拟合非线性模型后是否要拟合线性模型;平行要么是真正的或假指示是否应使用parallel。如运行NBAMSeq与平行= TRUE：
提取DE结果
DE分析结果可由结果函数。对于连续协变量，的名字参数应指定，以指示感兴趣的协变量。对于非线性连续协变量，将返回基本平均值、有效自由度(edf)、检验统计量、p值和调整后的p值。
DataFrame with 6行7列baseMean edf stat pvalue padj AIC BIC <数字> <数字> <数字> <数字> <数字> <数字> <数字> gene1 90.0788 1.00143 0.000371608 0.9952494 0.995249 217.150 224.121 gene2 80.7395 1.00025 1.114631449 0.2910806 192.597 199.568 gene3 42.4583 1.00008 0.016321577 0.8985639 0.916902 206.026 gene4 81.8683 1.00031 0.201885112 0.6537538 0.806844 216.570 225.066 gene6145.3781 1.00005 0.154988158 0.6938859 0.806844 239.810 246.781对于线性连续协变量，将返回基本平均值、估计系数、标准误差、检验统计量、p值和调整后的p值。
DataFrame with 6行8列baseMean coef SE stat pvalue padj AIC < 1.0124537 0.399229 -0.311943 0.388345 -2.753602 0.00589435 0.0982391 192.597 gene3 42.4583 -0.772035 0.340774 -2.265536 0.413360 -1.064868 0.5313624 216.570 gene5 80.2695 -0.705700 0.351417 -2.008156 0.04462673 0.1716413218.095 gene6 145.3781 -0.740655 0.334364 -2.215113 0.02675232 0.1230501 239.810 BIC gene1 224.121 gene2 199.568 gene3 213.026 gene4 223.540 gene5 225.066 gene6 246.781 对于离散协变量对比参数。如。对比度= c("var4"， "2"， "0")表示比较级别2和级别0var4．
DataFrame with 6行8列baseMean coef SE stat pvalue padj AIC <数值> <数值> <数值> <数值> <数值> <数值> <数值> gene1 90.0788 1.4089402 1.058595 1.3309531 0.1832044 0.436201 217.150 gene2 80.7395 0.3354875 1.0354018 0.3240117 0.7459292 0.925890 192.597 gene3 42.4583 -0.2760074 0.900537 -0.3064920 0.925890 206.056 gene4 81.8683 0.0691645 1.093739 0.0632368 0.9495779 0.961397 216.570 gene5 80.2695 -2.0069057 0.932441 -2.1523139 0.235583 218.095gene6 145.3781 0.8913475 0.887066 1.0048265 0.3149804 0.627228 239.810 BIC gene1 224.121 gene2 199.568 gene3 213.026 gene4 223.540 gene5 225.066 gene6 246.781 可视化
我们提出了两种方法来可视化非线性关联。第一种方法是绘制拟合的负二项式加法模型的平滑分量plot.gam在MGCV中的作用(《木材与木材》2015)．这可以通过打电话来完成makeplot函数和传入NBAMSeqDataSet对象。期望用户提供感兴趣的表型phenoname论点和基因的兴趣genename论点。
此外，为了探索协变量的非线性关联，查看对数归一化计数与变量散点图也是有指导意义的。下面我们将展示如何生成这样的图。
在这里，我们探索最重要的非线性关联res1 =res1 [订单(res1＄pvalue),)topgene =rownames(res1) [1］科幻小说=getsf(德牧)获得估计的大小因子将原始计数除以大小因子以获得标准化计数countnorm =t（t(countData)/科幻小说)头(res1)DataFrame with 6行7列baseMean edf stat pvalue padj AIC BIC gene26 110.2313 1.00010 13.22154 9.13066 0.00252013 0.0630032 205.670 212.640 gene40 106.3222 1.00015 6.35095 0.01174566 0.1957609 228.787 235.757 gene5 818.1699 1.00015 5.17247 0.2102943 240.353 247.323 gene3554.7496 1.00019 5.00944 0.02523531 0.2102943 207.447 214.418 图书馆(ggplot2)setTitle =topgenedf =data.frame（把=把,logcount =log2(countnorm [topgene,]+1）)ggplot(df,aes（x =把,y =logcount))+geom_point（形状=19，大小=1）+geom_smooth（方法=“黄土”）+xlab（“把”）+ylab（"log(normcount + 1)"）+注释（“文本”，x =马克斯(df＄把)-5，y =马克斯(df＄logcount)-1，标签=paste0（法国电力公司(edf):“，signif(res1 [topgene,“法国电力公司(edf)”]，数字=4)))+ggtitle(已)+主题（文本=element_text（大小=10)，情节。title =element_text（hjust =0．5）)
会话信息
R版本4.2.1(2022-06-23)平台:x86_64-pc-linux-gnu(64位)运行环境:Ubuntu 20.04.5 LTS矩阵产品:默认BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。so LAPACK: /home/biocbuild/bbs-3.16-bio /R/lib/libRlapack。so locale: [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC=C [3] LC_TIME=en_GB 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_TELEPHONE=C [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C附加的基础包:[1]stats4 stats graphics grDevices utils datasets methods[8]基础其他附加包:[1]ggplot2_3.3.6 BiocParallel_1.32.0 [3] NBAMSeq_1.14.0 SummarizedExperiment_1.28.0 [5] Biobase_2.58.0 GenomicRanges_1.50.0 [7] GenomeInfoDb_1.34.0 IRanges_2.32.0 [9] S4Vectors_0.36.0 BiocGenerics_0.44.0 [11] MatrixGenerics_1.10.0 matrixStats_0.62.0通过命名空间加载(且未附加):[1] httr_1.4.4 sass_0.4.2 bit64_4.0.5 [4] jsonlite_1.8.3 splines_4.2.1 bslib_0.4.0 [7] assertthat_0.2.1 highr_0.9 blob_1.2.3 [10] GenomeInfoDbData_1.2.9 yaml_2.3.6 pillar_1.8.1 [13] RSQLite_2.2.18 pillar_1.8.1 [13] RSQLite_2.2.18 lattice_0.20-45 glue_1.6.2 [16] digest_0.6.30 RColorBrewer_1.1-3 XVector_0.38.0 [19] colorspace_2.0-3 htmltools_0.5.3 Matrix_1.5-1 [22] DESeq2_1.38.0 XML_3.99-0.12 pkgconfig_2.0.3 [25] genefilter_1.80.0 zlibbioc_1.44.0 xtable_1.8 [28] scales_1.2.1 tibble_3.1.8 annotate_1.76.0 [31] mgcv_1. 41KEGGREST_1.38.0 farver_2.1.1 [34] generics_0.1.3 withr_2.5.0 cachem_1.0.6 [37] cli_3.4.1 survival_3.4-0 magrittr_2.0.3 [40] crayon_1.5.2 memoise_2.0.1 evaluate_0.17 [43] fansi_1.0.3 nlme_1 -160 tools_4.2.1 [46] lifecycle_1.0.3 string_1 .4.1 locfit_1. 1.5-9.6 [49] munsell_0.5.0 DelayedArray_0.24.0 AnnotationDbi_1.60.0 [52] Biostrings_2.66.0 compiler_4.2.1 jquerylib_0.1.4 [55] rlang_1.0.6 grid_4.2.1 RCurl_1.98-1.9 [58] labeling_0.4.2 bitops_1.0-7 rmarkdown_2.17 [61] gtable_0.3.1 codetools_0.2-18DBI_1.1.3 [64] R6_2.5.1 knitr_1.40 dplyr_1.0.10 [67] fastmap_1.1.0 bit_4.0.4 utf8_1.2.2 [70] stringi_1.7.8 parallel_4.2.1 Rcpp_1.0.9 [73] vctrs_0.5.0 geneplotter_1.76.0 png_0.1-7 [76] tidyselect_1.2.0 xfun_0.34参考文献
迪，Y, DW Schafer, JS Cumbie, Chang JH . 2015。“NBPSeq: rna测序数据的负二项式模型。”R软件包版本0.3。0, URL Http://CRAN。r项目。Org/Package = NBPSeq．
《爱》，迈克尔一世，沃尔夫冈·胡贝尔，西蒙·安德斯，2014。“用Deseq2对Rna-Seq数据的折叠变化和离散度进行调节估计。”基因组生物学15(12): 550。
马克·D·罗宾逊，戴维斯·J·麦卡锡，戈登·K·史密斯，2010。EdgeR:用于数字基因表达数据差异表达分析的生物导体包。生物信息学26(1): 139-40。
西蒙·伍德，西蒙·伍德，2015。“包Mgcv’。”R包版本1: 29。
周怡辉，夏凯，弗雷德·A·赖特，2011。“一个强大而灵活的方法来分析Rna序列计数数据。”生物信息学27(19): 2672-8。