内容

1简介

的mpra软件包提供了用于分析大规模并行报告分析(MPRA)数据的工具。具体来说，它包含中描述的功能(Myint et al. 2017)。主要的分析目的是使活动度量的差异分析成为可能，但是该包也可以用于生成在序列特征的活动分数的回归分析中有用的精确权重。的主要主力mpra包装是mpralm ()该功能借鉴了之前提出的用于RNA-seq分析的voom框架(Law et al. 2014)。在本文档中，我们将从一项比较episomal和慢病毒版本MPRA的研究中查看MPRA数据(Inoue et al. 2017)。

1．1如何引用

如果您正在使用此套餐，请注明(Myint et al. 2017)。如果您正在使用mpralm ()函数，也可以引用boom框架(Law et al. 2014)。

2依赖关系

该文档具有以下依赖项

库(mpra)

3.创建一个MPRASet对象

该包中包含MPRA数据MPRASet对象。因为MPRA数据没有通用的规定格式，所以必须手动创建这些对象。在本节中，我们将演示如何做到这一点。

MPRASet物体必须包含DNA和RNA计数信息，因为这是用于量化所测元素的活性水平的信息。DNA和RNA计数信息应指定为\(K \乘以S\)整数计数矩阵，其中\ (K \)如果提供条形码级别的信息，则为所有元素上的条形码总数;如果提供元素级别的信息，则为假定的监管元素(PREs)的总数。\ (\)是样本的数量(通常是独立转染的数量)。

MPRASet对象还必须包含元素标识信息。这应该作为长度的字符向量提供\ (K \)，即DNA和RNA计数矩阵中的行数。这些是用于描述/识别正在检测的唯一PREs的任何字符串。

可选地，可以将条形码序列和PRE序列指定为长度\ (K \)特征向量。

在下一节中，我们将提供如何为两种常见的差异分析设置(组织和等位基因比较)指定此信息的具体示例。虽然我们显示的是模拟数据，但这些信息通常是从文本文件中读取的。

3．1组织的比较

在组织比较研究中，同一组PREs在两种或两种以上的细胞类型中进行分析。在下面的例子中，实验观察四个PREs，每个PREs有三个条形码。研究了两个组织(肝脏和肾脏)，每个组织有4个重复(每个重复4个独立转染)。

RNA和DNA计数矩阵如下所示:

E <- 4 #元素数量B <- 3 #条形码数量s <- 4 #每个组织样本nt <- 2 #组织数量set.seed(434) rna <- matrix(rpois(E*B*s*nt, lambda = 30)， nrow = E*B, ncol = s*nt) dna <- matrix(rpois(E*B*s*nt, lambda = 30)， nrow = E*B, ncol = s*nt) rn <- as。character(outer(paste0("barcode_"， seq_len(B)， "_")， paste0("elem_"， seq_len(E))， FUN = "paste0") cn <- c(paste0("liver_"， seq_len(s))， paste0("肾_"，seq_len(s))) rownames(rna) <- rn rownames(dna) <- rn colnames(rna) <- cn colnames(dna) <- cn rna
# # liver_1 liver_2 liver_3 liver_4 kidney_1 kidney_2 kidney_3 # # barcode_1_elem_1 31日26日36 18 27 33 28 # # barcode_2_elem_1 33 28 29 32 28 33 34 # # barcode_3_elem_1 22 43 25 41 28 34 30 # # barcode_1_elem_2 28日27日31日20 37 33 24 # # 31 barcode_2_elem_2 32 32 21 38 22 27 # # barcode_3_elem_2 35 36 28 38 26日23日30 # # barcode_1_elem_3 30 34 26 34 40 26 26 # # barcode_2_elem_3 37 27 37 29 25 29 38 # # barcode_3_elem_3 34 30 20 26 30 35 33 # # barcode_1_elem_4 28 34 26日27日28 29 30 25 # # barcode_2_elem_4 31 323.6 31 24 29 ## barcode_3_elem_4 24 28 34 31 37 29 34 ## kidney_4 ## barcode_1_elem_1 30 ## barcode_2_elem_1 30 ## barcode_3_elem_1 37 ## barcode_1_elem_2 29 ## barcode_2_elem_2 29 ## barcode_3_elem_2 29 ## barcode_1_elem_3 22 ## barcode_2_elem_3 28 ## barcode_3_elem_3 38 ## barcode_1_elem_4 30 ## barcode_2_elem_4 32 ## barcode_3_elem_4 36
dna
# # liver_1 liver_2 liver_3 liver_4 kidney_1 kidney_2 kidney_3 # # barcode_1_elem_1 34 28 29日27日27日25 30 # # 43 barcode_2_elem_1 37 34 33 43 26 26 # # barcode_3_elem_1 35 23 32 27日19日29日29 # # barcode_1_elem_2 23 30 34 27 32 18 24 # # barcode_2_elem_2 32 29 42 21 37 34 32 # # barcode_3_elem_2 26日22日23日24日32 41 24 # # barcode_1_elem_3 29 38 34 27日31日27日31 # # barcode_2_elem_3 30 36 32 21 34 29 26 # # barcode_3_elem_3 34 28 42 35 26 32 32 # # barcode_1_elem_4 34 25 29 30 44 26日23 # # barcode_2_elem_4 29 26 373.7 28 37 34 ## barcode_3_elem_4 24 28 33 40 31 30 19 ## kidney_4 ## barcode_1_elem_1 31 ## barcode_2_elem_1 33 ## barcode_3_elem_1 25 ## barcode_1_elem_2 25 ## barcode_2_elem_2 31 ## barcode_3_elem_2 23 ## barcode_1_elem_3 33 ## barcode_2_elem_3 39 ## barcode_3_elem_3 26 ## barcode_1_elem_4 33 ## barcode_2_elem_4 31 ## barcode_3_elem_4 26

PRE标识字符串如下所示。当在条形码级别提供计数时，开斋节字符向量将有重复的元素。

eid <- rep(paste0("elem_"， seq_len(E))， each = B) eid
# #[1]“elem_1”“elem_1”“elem_1”“elem_2”“elem_2”“elem_2”“elem_3”“elem_3”# #[9]“elem_3”“elem_4”“elem_4”“elem_4”

我们也可以有PRE序列如下所示。这些序列必须在长度相同的字符向量中指定开斋节和相同的行数核糖核酸和dna。

eseq <- replication (E, paste(sample(c("A"， "T"， " c "， "G")， 10, replace = TRUE)， collapse = "")) eseq <- rep(eseq, each = B) eseq . copy (E, paste(sample(c("A"， "T"， " c "， "G")
##[1]“tgaccatacc”“tgaccatacc”“tgaccatacc”“acgcccaatt”“acgcccaatt”##[6]“acgcccaatt”“cggcgagggg”“cggcgagggg”“cggcgagggg”“cactgtacta”##[11]“cactgtacta”“cactgtacta”

以上作品(核糖核酸，dna，开斋节,eseq类的参数MPRASet构造函数的功能如下所示。如果条形码(条形码序列)或eseq(PRE序列)不提供，必须指定为零。

mpraset_example <- mprasetet (DNA = DNA, RNA = RNA, eid = eid, eseq = eseq, barcode = NULL
##类:MPRASet ## dim: 12 8 ##元数据(0):## assays(2): DNA RNA ## rownames(12): barcode_1_elem_1 barcode_2_elem_1…barcode_2_elem_4 ## barcode_3_elem_4 ## rowData names(2): eid eseq ## colnames(8): liver_1 liver_2…kidney_3 kidney_4 ## colData names(0): ##不存在条形码

3.2等位基因比较

在等位基因比较研究中，存在两个或两个以上等位基因的PREs被分析。所有PREs的等位基因版本都在同一样本中进行了检测。因为需要等位基因之间的活性比较，所以这些计数必须分开到不同的列中。在下面的例子中，实验观察四个PREs，每个PREs有三个条形码。每个PRE有两个等位基因，有四个重复(总共四个独立转染)。

RNA和DNA计数矩阵如下所示:

E <- 4 #元素数目B <- 3 #条形码数目s <- 4 #样本总数目nalleles <- 2 #等位基因数目set.seed(434) rna <- matrix(rpois(E*B*s*nalleles, lambda = 30)， nrow = E*B*s*nalleles) dna <- matrix(rpois(E*B*s*nalleles, lambda = 30)， nrow = E*B, ncol = s*nalleles) rn <- as。character(outer(paste0("barcode_"， seq_len(B)， "_")， paste0("elem_"， seq_len(E))， FUN = "paste0") cn <- c(paste0("allele1_sample"， seq_len(s))， paste0("allele2_sample"， seq_len(s))) rownames(rna) <- rn rownames(dna) <- rn colnames(rna) <- cn colnames(dna) <- cn colnames(dna) <- cn rna
# # # # allele1_sample1 allele1_sample2 allele1_sample3 barcode_1_elem_1 31日26日36 # # barcode_2_elem_1 33 28 29 # # barcode_3_elem_1 22 43 25 # # barcode_1_elem_2 28 27 31 # # barcode_2_elem_2 32 31 32 # # barcode_3_elem_2 35 36 28 # # barcode_1_elem_3 30 34 26 # # barcode_2_elem_3 37 27 37 # # barcode_3_elem_3 34 30 20 # # barcode_1_elem_4 28 27 34 # # barcode_2_elem_4 31 32 28 # # barcode_3_elem_4 24 28 34 # # allele1_sample4 allele2_sample1 allele2_sample2 # # barcode_1_elem_1 18 27 33 # # barcode_2_elem_132 28 33 # # barcode_3_elem_1 41 28 34 # # barcode_1_elem_2 20 37 33 # # barcode_2_elem_2 21 38 22 # # barcode_3_elem_2 38 26日23 # # barcode_1_elem_3 34 40 26 # # barcode_2_elem_3 29 25 29 # # barcode_3_elem_3 26 30 35 # # barcode_1_elem_4 26 29 30 # # barcode_2_elem_4 36 31 24 # # barcode_3_elem_4 31 37 29 # # allele2_sample3 allele2_sample4 30 # # # # barcode_1_elem_1 28 barcode_2_elem_1 34 30 # # barcode_3_elem_1 30 37 # # barcode_1_elem_2 24 29 # # barcode_2_elem_2 27 29 # # barcode_3_elem_2 30 29 # #Barcode_1_elem_3 26 22 ## barcode_2_elem_3 38 28 ## barcode_3_elem_3 33 38 ## barcode_1_elem_4 25 30 ## barcode_2_elem_4 29 32 ## barcode_3_elem_4 34 36
dna
# # allele1_sample1 allele1_sample2 allele1_sample3 # # barcode_1_elem_1 29 34 28 # # barcode_2_elem_1 43 37 34 # # barcode_3_elem_1 35 23 32 # # barcode_1_elem_2 23 30 34 # # barcode_2_elem_2 32 29 42 # # barcode_3_elem_2 26日23日22 # # barcode_1_elem_3 29 38 34 # # barcode_2_elem_3 30 36 28 42 # # 32 # # barcode_3_elem_3 34 barcode_1_elem_4 34 25 29 # # barcode_2_elem_4 29 26 37 # # barcode_3_elem_4 24 28 33 # # allele1_sample4 allele2_sample1 allele2_sample2 # # 25 # # barcode_2_elem_1 barcode_1_elem_1 27 2733 43 26 # # barcode_3_elem_1 27日19日29 # # barcode_1_elem_2 27 32 18 # # barcode_2_elem_2 21 37 34 # # barcode_3_elem_2 24 32 41 # # barcode_1_elem_3 27日31日27 # # barcode_2_elem_3 21 34 29 # # barcode_3_elem_3 35 26 32 # # barcode_1_elem_4 30 44 26 # # barcode_2_elem_4 37 28 37 # # barcode_3_elem_4 40 31 30 # # allele2_sample3 allele2_sample4 30 31 # # # # barcode_1_elem_1 barcode_2_elem_1 26 33 # # barcode_3_elem_1 29 25 # # barcode_1_elem_2 24 25 # # barcode_2_elem_2 32 31 # # barcode_3_elem_2 24 23 # #Barcode_1_elem_3 31 33 ## barcode_2_elem_3 26 39 ## barcode_3_elem_3 32 26 ## barcode_1_elem_4 23 33 ## barcode_2_elem_4 34 31 ## barcode_3_elem_4 19 26

我们可以像上面那样定义PRE标识字符串和序列，并在MPRASet构造函数类似:

mpraset_example2 <- mpraset_example2 (DNA = DNA, RNA = RNA, eid = eid, eseq = eseq, barcode = NULL
##类:MPRASet ## dim: 12 8 ##元数据(0):## assays(2): DNA RNA ## rownames(12): barcode_1_elem_1 barcode_2_elem_1…barcode_2_elem_4 ## barcode_3_elem_4 ## rowData names(2): eid eseq ## colnames(8): allele1_sample1 allele1_sample2…allele2_sample3 ## allele2_sample4 ## colData names(0): ##不存在条形码

4分析

而上面的部分演示了如何创建MPRASet对象，我们将使用预先构建的对象，其中包含来自MPRA的episomal和慢病毒版本比较的数据(Inoue et al. 2017)。

数据(mpraSetExample)

我们创建了带有episomal(突变整合酶)样本指标的设计矩阵，并将精确加权线性模型与mpralm。在MPRA实验中，活性测量被量化为RNA计数与DNA计数的log2比。当存在条形码级别的信息时(如本实验中所示)，有多种方法汇总条形码上的信息，以计算用于后续统计建模的最终元素和特定样本的日志比率。

我们已经明确了聚合= "均值"以表明元素和样本特定的对数比将通过首先计算每个条形码的RNA计数与DNA计数的对数比，然后对特定元素和样本中的条形码取平均值来计算。

我们已经明确了normalize = TRUE对RNA和DNA文库进行总计数归一化。它将所有库扩展为具有1000万次读取的共同大小。

因为这个实验是在两种不同的细胞条件下观察一组PREs，不同的样本(列MPRASet对象)是独立的。因此我们指定Model_type = "indep_groups"执行未配对分析。相反，如果我们正在进行等位基因比较，我们将指定Model_type = "corr_groups"进行了配对分析(表明的不同列MPRASet对象是链接的)。

最后我们指定情节=真实绘制对数比变异性与元素拷贝数之间的关系。

design <- data.frame(intcpt = 1, episomal = grepl("MT"， colnames(mpraSetExample))) mpralm_fit <- mpralm(object = mpraSetExample, design = design, aggregate = "mean"， normalize = TRUE, model_type = "indep_groups"， plot = TRUE)

生成的适合对象可用于topTable从limma包中。

toptab <- topTable(mpralm_fit, coef = 2, number = Inf)

因为这个数据集的元素代码相当长，我们使用了一些技巧来打印顶部的差分元素:

rownames (toptab6)
## [3] "C:SLEA_hg18:chr2:210861483-210861650|5:V_AHRARNT_02:GGGGATCGCGTGCCAGCCC; 45:V_HNF1_C:AGTTAATGATTAACCAA;64:V_AHRARNT_02: ggggatcgtgccagccc;85:V_HNF1_C:AGTTAATGATTAACCAA;104:V_AHRARNT_02: ggggatcgtgccagccc;125:V_HNF1_C:AGTTAATGATTAACCAA;144:V_AHRARNT_02:GGGGATCGCGTGCCAGCCC" ## [2] "R: ep300 - nomod_chr9:12814543-12814714]" ## [3]“C:SLEA_hg18:chr2:210861483-210861650|4:V_AHRARNT_02:GGGGATCGCGTGCCAGCCC; 46:V_AHRARNT_02: ggggatcgtgccagccc;67:V_HNF1_C:AGTTAATGATTAACCAA;86:V_AHRARNT_02: ggggatcgtgccagccc;107:V_HNF1_C:AGTTAATGATTAACCAA;126:V_AHRARNT_02: ggggatcgtgccagccc;147:V_HNF1_C:AGTTAATGATTAACCAA”## [4]"C:SLEA_hg18:chr2:210861483-210861650|6:V_AHRARNT_02: ggggatcgtgccagccc;27:V_HNF1_C:AGTTAATGATTAACCAA; 65:V_AHRARNT_02: ggggatcgtgccagccc;86:V_HNF1_C:AGTTAATGATTAACCAA;105:V_HNF1_C:AGTTAATGATTAACCAA;124:V_HNF1_C:AGTTAATGATTAACCAA;143:V_AHRARNT_02: ggggatcgtgccagccc " ## [5] "R:FOXA1_FOXA2-ChMod_chr1:38000211-38000353_[chr1:38000196-38000367]" ## [6] "R: ep300 - nomod_chr16:3070958-3071129]"
rownames(toptab6) <- NULL toptab6
## logFC AveExpr t P.Value adj.P.Val B ## 1 -1.2076627 1.284363 -26.99214 1.450176e-36 3.538429e-33 73.05825 ## 2 -1.7316415 1.668562 -23.50738 4.357926e -30 64.87553 ## 3 -1.1603766 1.497371 -23.11093 1.150733e-32 9.359292e-30 64.12459 ## 4 -0.9560848 1.016586 -18.57491 - 18.136470e -27 1.303247e-24 52.09394 ## 5 -1.1713718 1.730330 -18.25135 5.494116e-27 2.681128e-24 51.09648 ## 6 -0.7639104 0.733524 -17.71599 2.689472e-26 1.093719e-23 49.54888

5Sessioninfo

R版本3.5.0(2018-04-23)平台:x86_64-pc-linux-gnu(64位)运行在Ubuntu 16.04.4 LTS下

矩阵产品:默认BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas。so LAPACK: /home/biocbuild/bbs-3.7-bioc/R/lib/libRlapack.so

locale: [1] LC_CTYPE=en_US。utf - 8 LC_NUMERIC = C
[3]而= en_US。utf - 8 LC_COLLATE = C
[5] LC_MONETARY = en_US。utf - 8LC_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]方法基础

其他附件:[1]mpra_1.2.0 limma_3.36.0
[3] SummarizedExperiment_1.10.0 DelayedArray_0.6.0
[5] BiocParallel_1.14.0 matrixStats_0.53.1
[7] Biobase_2.40.0 genome icranges_1 .32.0
[9] GenomeInfoDb_1.16.0 IRanges_2.14.0
[11] S4Vectors_0.18.0 BiocGenerics_0.26.0
[13] BiocStyle_2.8.0

通过命名空间加载(且未附加):[1]Rcpp_0.12.16 compiler_3.5.0 plyr_1.8.4
[4] XVector_0.20.0 bitops_1.0-6 tools_3.5.0
[7] zlibbioc_1.26.0 statmod_1.4.30 digest_0.6.15
[10] evaluate_0.10.1 lattice_0.20-35 Matrix_1.2-14
[13] yaml_2.1.18 xfun_0.1 GenomeInfoDbData_1.1.0 [16] string_1 .3.0 knitr_1.20 rprojroot_1.3-2 .2 .
[19] grid_3.5.0 rmarkdown_1.9 bookdown_0.7
[22] magrittr_1.5 backports_1.1.2 scales_0.5.0
[25] htmltools_0.3.6 colorspace_1.3-2 stringi_1.1.7
[28] RCurl_1.95-4.10 munsell_0.4.3 .

参考文献

Inoue, Fumitaka, Martin Kircher, Beth Martin, Gregory M Cooper, Daniela M Witten, Michael T McManus, Nadav Ahituv和Jay Shendure。2017。“系统比较揭示了增强子活性的染色体编码与插曲编码的实质性差异。”基因组研究27:38-52。https://doi.org/10.1101/gr.212092.116。

陈云顺，史伟，高登·K·史密斯。2014。“Voom:精确权重解锁RNA-seq读取计数的线性模型分析工具。”基因组生物学15: R29。https://doi.org/10.1186/gb-2014-15-2-r29。

Myint, Leslie, Dimitrios G Avramopoulos, Loyal A. Goff, Kasper D Hansen, 2017。“线性模型在大规模并行报告分析中实现强大的微分活性分析。”bioRxiv, 196394年。https://doi.org/10.1101/196394。