1简介

在这个工作流程中，我们检查来自RNA-seq和ATAC-seq数据的子集Alasoo等(2018）该研究涉及用干扰素γ (IFNg)治疗来自许多人类供体的巨噬细胞系，沙门氏菌感染，或者两种治疗结合。Alasoo等(2018）检测了86个成功分化诱导多能干细胞(iPSC)系的基因表达和染色质可及性，并比较了特定数量性状位点(QTL)染色质可及性和基因表达的基线和响应。作者发现许多刺激特异性表达QTL已经在未培养细胞中作为染色质QTL检测到，并进一步假设了与刺激反应相关的转录因子的性质和作用。

我们将进行一个比在Alasoo等(2018）他们使用公开的RNA-seq和ATAC-seq数据(忽略基因型)。我们将研究IFNg刺激对基因表达和染色质可达性的影响，并观察差异表达(DE)基因附近是否有差异可达(DA) ATAC-seq峰的富集。这是合理的，因为IFNg刺激的转录组反应可能是通过调控蛋白与可达区域的结合来介导的，这种结合可能会增加这些区域的可达性，从而可以被ATAC-seq检测到。

流畅的基因组学工作流程概述。首先，我们将数据作为summarizeexperiment对象导入，它支持与下游分析包的互操作性。然后，我们使用现有的Bioconductor包DESeq2和limma对我们的分析数据进行建模。我们采取我们的模型的结果，每个化验相对于他们的基因组坐标，并整合他们。首先，我们计算每个检测结果之间的重叠，然后在组合的基因组区域上进行聚合，最后总结比较差异表达基因与非差异表达基因的富集程度。最终输出可用于下游可视化或进一步转换。

图1.1:基因组学流畅工作流程概述首先,我们进口数据作为SummarizedExperiment对象，它支持与下游分析包的互操作性。然后,我们模型我们的检测数据，使用现有的Bioconductor包DESeq2而且limma．我们采取我们的模型的结果，每个化验相对于他们的基因组坐标，和集成他们。首先，我们计算每个检测结果之间的重叠，然后在组合的基因组区域上进行聚合，最后总结比较差异表达基因与非差异表达基因的富集程度。最终输出可用于下游可视化或进一步转换。

在整个工作流中(图1．1)，我们将使用现有的Bioconductor基础设施来理解这些数据集。我们将特别强调Bioconductor包的使用plyranges而且tximeta．的plyranges使用移动、窗口构建、重叠检测等操作，将绑定到基因组范围的数据进行流畅地转换。它被描述为李、库克和劳伦斯(2019）并利用底层核心Bioconductor基础设施劳伦斯等。2013；Huber等人。2015）和tidyverse设计原则维克汉姆等人(2019）．

的tximeta包由爱等人。2019）用于将RNA-seq定量数据读取到R/Bioconductor中，从而将转录本范围及其来源自动附加到包含表达值和差异表达结果的对象上。

1．1实验数据

此工作流中使用的数据可从两个包获得巨噬细胞Bioconductor实验数据包和工作流包fluentGenomics．

的巨噬细胞包包含来自24个RNA-seq样本的RNA-seq定量，RNA-seq样本的一个子集生成和分析Alasoo等(2018）．对末端读数进行量化大马哈鱼(Patro et al。2017），使用Gencode 29人类参考转录本(Frankish, gencode - consortium，和Flicek2018）．有关量化的更多细节，以及所使用的确切代码，请参阅巨噬细胞包中。该包还包含Snakemake文件，用于分发大马哈鱼集群的量化任务(Köster和拉赫曼2012）．

的fluentGenomics包包含下载和生成缓存的功能SummarizedExperiment对象提供的规范化ATAC-seq数据阿拉索及加夫尼(2017）．该对象包含所有实验条件下的所有145个ATAC-seq样本Alasoo等(2018）．数据也可以直接从Zenodo沉积。

下面的代码加载缓存数据文件的路径，如果不存在，则创建缓存并生成SummarizedExperiment使用theBiocFileCache包(谢博德和摩根2019）．

然后，我们可以读取缓存的文件，并将其分配给一个名为atac．

图书馆(fluentGenomics)atac < -readRDS（cache_atac_se（))

对我们如何得到这个的精确描述SummarizedExperiment对象可以在section中找到2．2．

2导入数据SummarizedExperiment

2.1使用tximeta导入RNA-seq量化数据

首先，我们指定一个目录dir，其中存放量化文件。你可以简单地指定这个目录:

dir < -“/路径/ /定量/文件”

其中路径相对于当前R会话。但是，在这种情况下，我们将文件分发到巨噬细胞包中。相关目录和相关文件可以使用执行．

dir < -执行（“extdata”，包=“巨噬细胞”）

关于实验的信息包含在coldata.csv文件。我们利用dplyr而且readr包(作为tidyverse)将该文件读入R维克汉姆等。2019）．我们稍后会看到plyranges扩展这些包以适应基因组范围。

图书馆(dplyr)

## ##附加包:'dplyr'

以下对象从'package:stats'中屏蔽:## ## filter, lag

以下对象从'package:base'中屏蔽:## ## intersect, setdiff, setequal, union

图书馆(readr)colfile < -file.path(dir,“coldata.csv”）coldata < -read_csv(colfile)% > %dplyr：：选择（的名字,id =sample_id,行=line_id,条件=condition_name）% > %dplyr：：变异（文件=file.path(dir,“宽客”、名称、“quant.sf.gz”)，行=因素(线),条件=relevel（因素(条件),“幼稚”））

##行:24列:13

列规范# #──────────────────────────────────────────────────────────# #分隔符:”、“# #杆(10):姓名、sample_id, line_id, condition_name macrophage_harvest, sal…## dbl (3): replication, ng_ul_mean, rna_auto ## ##ℹ使用' spec() '检索该数据的完整列规范。##ℹ指定列类型或设置' show_col_types = FALSE '来关闭此消息。

coldata

## #小猫咪:24×5 # # # #名称标识线条件文件<空空的> <空空的> < fct > < fct > < >从而向# # 1 SAMEA103885102 diku_A diku_1天真/home/biocbuild/bbs - 3.16 - bioc / R / lib…# # 2 SAMEA103885347 diku_B diku_1 IFNg /home/biocbuild/bbs - 3.16 - bioc / R / lib…# # 3 SAMEA103885043 diku_C diku_1 SL1344 /home/biocbuild/bbs - 3.16 - bioc / R / lib…# # 4 SAMEA103885392 diku_D diku_1 IFNg_SL1344 /home/biocbuild/bbs - 3.16 - bioc / R / lib…# # 5 SAMEA103885182 eiwy_A eiwy_1天真/home/biocbuild/bbs - 3.16 - bioc / R / lib…# # 6 SAMEA103885136 eiwy_B eiwy_1## 7 SAMEA103885413 eiwy_C eiwy_1 SL1344 /home/biocbuild/bbs-3.16-bioc/R/lib…## 8 SAMEA103884967 eiwy_D eiwy_1 IFNg_SL1344 /home/biocbuild/bbs-3.16-bioc/R/lib…## 9 SAMEA103885368 fikt_A fikt_3 naive /home/biocbuild/bbs-3.16-bioc/R/lib…## 10 SAMEA103885218 fikt_B fikt_3 IFNg /home/biocbuild/bbs-3.16-bioc/R/lib

在我们读过之后coldata.csv文件中，我们从这个表中选择相关的列，创建一个名为文件，并改造现有的行而且条件列成因子。在这种情况下条件，我们指定“幼稚”细胞系作为参考级别。的文件列指向每个观察结果的量化-这些文件已被gzip，但如果使用from，通常不会有' gz '结尾大马哈鱼直接。另一件需要注意的事情是管道操作符的使用，% > %，可读为“then”，即先读取数据，然后选择列,然后变异。

现在我们有一个总结实验设计和量化的位置的表格。以下几行代码为分析人员做了大量工作:导入RNA-seq量化(删除推论复制在本例中)，定位相关的参考转录组，将转录范围附加到数据上，并获取基因组信息。推论复制对于进行转录水平分析特别有用，但在这里我们将对每个基因计数使用点估计并进行基因水平分析。

结果是SummarizedExperiment对象。

suppressPackageStartupMessages（图书馆(SummarizedExperiment))图书馆(tximeta)se < -tximeta(coldatadropInfReps =真正的）

导入量化

##读取read_tsv文件

# # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21日22日23日24 # #发现匹配与转录组:# # (GENCODE——智人——发布29)# # useHub = TRUE:检查TxDb通过“AnnotationHub”# # snapshotDate(): 2022-10-26 # #发现匹配TxDb通过“AnnotationHub”从缓存加载# # # #加载所需的包:GenomicFeatures # #加载所需的包:AnnotationDbi # # # #附加包:“AnnotationDbi”掩盖了# # # #以下对象从“包:dplyr”:## ##选择## ##生成记录范围

se

##类:rangedsummarizeexperiment ## dim: 205870 24 ## metadata(6): tximetaInfo quantInfo…txomeInfo txdbInfo ## assays(3):计数丰度长度## rownames(205870): ENST00000456328.2 ENST00000450305.2…## ENST00000387460.2 ENST00000387461.2 ## rowData names(3): tx_id gene_id tx_name ## colnames(24): SAMEA103885102 SAMEA103885347…SAMEA103885308 ## SAMEA103884949 ## colData names(4): names id行条件

在有internet连接的机器上，上面的命令不需要任何额外的步骤就可以工作tximeta函数通过FTP获取任何必要的元数据，除非它已经在本地缓存。的tximeta包也可以在没有互联网连接的情况下使用，在这种情况下，链接的转录组可以直接从一个大马哈鱼索引和gtf。

makeLinkedTxome（indexDir =file.path(dir,“gencode.v29_salmon_0.12.0”)，源=“Gencode”，生物=“智人”，释放=“29”，基因组=“GRCh38”，fasta =“ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.transcripts.fa.gz”，gtf =file.path(dir,“gencode.v29.annotation.gtf.gz”)，#本地版本写=假）

因为tximeta知道正确的参考转录组，我们可以问tximeta用的方法将转录水平的数据总结到基因水平Soneson, Love and Robinson (2015）．

gse < -summarizeToGene(se)

##加载现有TxDb创建:2022-11-01 15:03:29

从数据库中获取转录到基因的图谱

##生成基因范围

##总结丰盛

##汇总计数

##汇总长度

最后要注意的是开始正链基因和结束负链基因的遗传程度现在是由基因的同工型的基因组程度决定的开始而且结束还原性的农庄)．另一种选择是对转录丰度进行操作，并对转录进行差异分析(因此避免定义一组异构体的TSS)，或者使用基因水平的总结表达，但根据异构体表达选择最具代表性的TSS。

2．2导入ATAC-seq数据SummarizedExperiment对象

的SummarizedExperiment包含ATAC-seq峰值的对象可以从以下以制表符分隔的文件中创建阿拉索及加夫尼(2017）：

示例元数据:ATAC_sample_metadata.txt.gz(< 1米)
标准化读计数矩阵:ATAC_cqn_matrix.txt.gz(109米)
注释的峰值:ATAC_peak_metadata.txt.gz(5.6米)

首先，我们读入样例元数据，遵循与生成coldataRNA-seq实验表:

atac_coldata < -read_tsv（“ATAC_sample_metadata.txt.gz”）% > %选择（sample_id,捐赠,条件=condition_name）% > %变异（条件=relevel（因素(条件),“幼稚”）)

ATAC-seq计数已经标准化cqn(Hansen, Irizarry, and Wu2012）log2的变换。加载cqn-normalized matrix of log2转换后的读计数需要大约30秒，并加载一个大约370 Mb的对象。我们设置列名，使第一列包含矩阵的行名，其余列是来自atac_coldata对象。

atac_mat < -read_tsv（“ATAC_cqn_matrix.txt.gz”，跳过=1，col_names =c（“rownames”, atac_coldata [[“sample_id”]]))rownames < -atac_mat [[“rownames”]]atac_mat < -as.matrix(atac_mat [,-1]）rownames(atac_mat) < -rownames

我们读取峰值元数据(基因组中的位置)，并将其转换为a农庄对象。的as_granges ()函数自动转换data.frame成一个农庄对象。从这个结果中，我们提取peak_id列，并将基因组信息设置为构建“GRCh38”。我们从Zenodo条目．

图书馆(plyranges)peaks_df < -read_tsv（“ATAC_peak_metadata.txt.gz”，col_types =c（“cidciicdc”））peaks_gr < -peaks_df% > %as_granges（seqnames =科)% > %选择（peak_id =gene_id)% > %set_genome_info（基因组=“GRCh38”）

最后，我们构造了一个SummarizedExperiment对象。我们将矩阵作为命名列表放入assays槽中，注释的峰值放入按行排列的范围槽中，样本元数据放入按列排列的数据槽中:

atac < -SummarizedExperiment（化验=列表（cqndata =atac_mat),rowRanges =peaks_gr,colData =atac_coldata)

3.模型分析

3．1RNA-seq差异基因表达分析

我们可以很容易地运行微分表达式分析DESeq2使用以下代码块(Love, Huber, and Anders2014）．设计公式表明我们希望控制供体基线(行)，并测试基因表达的差异。有关Bioconductor中DE工作流的更全面的讨论，请参阅爱等人。2016）而且法律等(2018）．

图书馆(DESeq2)dds < -DESeqDataSet(gse,～行+条件)

##使用tximeta的计数和平均转录长度

#过滤低表达基因#在至少6个样本中至少有10个计数保持< -rowSums（计数(dds)> =10）> =6dds < -dds [,]

该模型适用于以下代码行:

dds < -DESeq(dds)

估算尺寸因素

##使用'avgTxLength'从分析(dds)，纠正库大小

估算离散度

##基因分散估计

均值-色散关系

最终的离散度估计

装配模型和测试

下面我们在条件变量上设置对比，表明我们正在估计\ \ log_2 \ ()IFNg刺激细胞系对未培养细胞系的折叠变化(LFC)。我们感兴趣的是LFC大于1，名义错误发现率(FDR)为1%。

res < -结果(dds,对比=c（“条件”，“IFNg”，“幼稚”)，lfcThreshold =1，α=0.01）

为了查看表达式分析的结果，我们可以生成一个汇总表和一个MA图:

总结(res)

## ##出17806非零总读计数##调整p值< 0.01 ## LFC > 1.00(上升):502,2.8% ## LFC < -1.00(下降):247,1.4% ##异常值[1]:0,0% ##低计数[2]:0,0% ##(平均计数< 3)##[1]参见'cooksCutoff'参数?results ##[2]参见'independentFiltering'参数?results

DESeq2：：plotMA(res,ylim =c（-10，10）)

图3.1:可视化DESeq2结果为“MA图”。有调整的基因假定值低于0.01则为红色。

现在将结果输出为a农庄对象，由于约定plyranges，我们构造一个名为gene_id从结果的行名。现在每一行都包含基因组区域(seqnames，开始，结束，链)以及相应的元数据列gene_id以及测试结果)。请注意,tximeta已正确识别出参考基因组为“hg38”，这也已添加到农庄沿着结果列。一旦执行重叠操作，这种簿记是至关重要的，以确保plyranges不是在不兼容的基因组之间进行比较。

suppressPackageStartupMessages（图书馆(plyranges))de_genes < -结果(dds,对比=c（“条件”，“IFNg”，“幼稚”)，lfcThreshold =1，格式=“农庄”）% > %names_to_column（“gene_id”）de_genes

与17806范围和7 # #农庄对象元数据列:# # seqnames范围链| gene_id baseMean # # < Rle > < IRanges > < Rle > | <人物> <数字> # # [1]chrX 100627109 - 100639991 | ENSG00000000003.14 171.571 # # [2] chr20 50934867 - 50958555 | ENSG00000000419.12 967.751 # # [3] chr1 169849631 - 169894267 | ENSG00000000457.13 682.433 # # [4] chr1 169662007 - 169662007 + | ENSG00000000460.16 262.963 # # [5] chr1 27612064 - 27635277 | ENSG00000000938.12 2660.102  ## ... ... ... ... . ... ...# # [17802] chr10 84167228 - 84172093 | ENSG00000285972.1 10.04746 # # [17803] chr6 63572012 - 63572012 + | ENSG00000285976.1 4586.34617 # # [17804] chr16 57177349 - 57177349 + | ENSG00000285979.1 14.29653 # # [17805] chr8 103398658 - 103501895 | ENSG00000285982.1 27.76296 # # [17806] chr10 12563151 - 12563151 + | ENSG00000285994.1 6.60409 # # log2FoldChange lfcSE stat pvalue padj # # <数字> <数字> <数字> <数字> <数字> # # [1]-0.2822450 0.3005710 0.00000 1.0000000 1.000000 0.0391223 # # [2][4] -1.4718762 0.2186916 -2.15772 0.0309493 0.409728 ## [5] 0.6754781 0.2360530 0.00000 1.0000000 1.000000 ## ... ... ... ... ... ...## [17802] 0.5484518 0.444319 0 1 1 ## [17803] -0.0339296 0.188005 0 1 1 ## [17804] 0.3123477 0.522700 0 1 1 ## [17805] 0.9945187 1.582373 0 1 1 ## [17806] 0.2539975 0.595751 0 1 1 ## ------- seqinfo:来自hg38基因组的25个序列(1个循环)

由此，我们可以将结果限制为满足FDR阈值的结果，并选择(并重命名)我们感兴趣的元数据列:

de_genes < -de_genes% > %过滤器(padj<0.01）% > %选择(gene_idde_log2FC =log2FoldChange,de_padj =padj)

我们现在希望提取出有证据表明LFC存在的基因不大。我们通过指定LFC阈值和替代假设(altHypothesis)表示LFC的绝对值小于阈值。要可视化此测试的结果，可以运行结果没有格式= "农庄组织"，并将此对象传递给plotMA像以前一样。我们把这些基因标记为other_genes后来被称为“非de基因”，以便与我们的基因进行比较de_genes集。

other_genes < -结果(dds,对比=c（“条件”，“IFNg”，“幼稚”)，lfcThreshold =1，altHypothesis =“lessAbs”，格式=“农庄”）% > %过滤器(padj<0.01）% > %names_to_column（“gene_id”）% > %dplyr：：选择(gene_idde_log2FC =log2FoldChange,de_padj =padj)

3.2ATAC-seq峰差丰度分析

下面的部分描述了我们用来生成农庄对象的差异峰从ATAC-seq数据Alasoo等(2018）．

本节其余部分的代码块将不运行。

为了评估不同的可访问性，我们运行limma(史密斯2004），并生成lfc和调整后的峰值p值的汇总:

图书馆(limma)设计< -model.matrix（～捐赠+条件下,colData(atac))符合< -lmFit（分析(atac),设计)符合< -易趣(适合)idx < -哪一个（colnames(适合＄系数)= =“conditionIFNg”）tt < -topTable(健康,系数=idx,sort.by =“没有”，n =nrow(atac))

现在我们用rowRanges的SummarizedExperiment并附上lfc和调整后的p值limma，使我们可以考虑与微分表达式的重叠。请注意，我们将基因组构建设置为“hg38”，并将染色体信息重新设置为“UCSC”样式(例如“chr1”，“chr2”等)。再一次，我们从Zenodo条目的ATAC-seq数据中知道了基因组的构建。

atac_peaks < -rowRanges(atac)% > %remove_names（）% > %变异（da_log2FC =tt＄logFC,da_padj =tt＄adj.P.Val）% > %set_genome_info（基因组=“hg38”）seqlevelsStyle(atac_peaks) < -“加州”

最后一个农庄包含DA峰值的对象包含在工作流包中，可以按如下方式加载:

图书馆(fluentGenomics)山峰

## GRanges对象，296220范围和3个元数据列:## seqnames范围strand | peak_id da_log2FC ##    | <字符> > ## [1]chr1 9979-10668 * [1] ATAC_peak_1 0.266185 ## [2] chr1 10939-11473 * | ATAC_peak_2 0.322177 ## [3] chr1 15505-15729 * | ATAC_peak_3 -0.574160 ## [4] chr1 21148-21481 * | ATAC_peak_4 -1.147066 ## [5] chr1 21864-22067 * | ATAC_peak_5 -0.896143 ## ... ... ... ... . ... ...# # [296216] chrX 155896572 - 155896572 * | ATAC_peak_296216 -0.834629 # # [296217] chrX 155958507 - 155958507 * | ATAC_peak_296217 -0.147537 # # [296218] chrX 156016760 - 156016760 * | ATAC_peak_296218 -0.609732 # # [296219] chrX 156028551 - 156028551 * | ATAC_peak_296219 -0.347678 # # [296220] chrX 156030135 - 156030135 * | ATAC_peak_296220 0.492442 # # da_padj # # <数字> 9.10673 e-05 # # # # [1] [2] 2.03435 e-05 e-08 # # 3.41708 # # [3] [4] 8.22299 e-26 # # 4.79453 e-11 [5]  ## ... ...## [296216] 1.33546e-11 ## [296217] 3.13015e-01 ## [296218] 3.62339e-09 ## [296219] 6.94823e-06 ## [296220] 7.07664e-13 ## ------- ## seqinfo:来自hg38基因组的23个序列;没有seqlengths

4集成的范围

4．1找到与plyranges

我们已经用过plyranges多次以上，才能过滤器，变异,选择在农庄对象，以及确保正确的基因组注释和风格已被使用。的plyranges包提供了执行基因组数据转换的语法(李，库克，劳伦斯2019）．的组合所产生的计算plyranges中的底层高度优化的范围操作执行“动词”GenomicRanges包劳伦斯等。2013）．

对于重叠分析，我们过滤注释的峰值，使其名义FDR界为1%。

da_peaks < -山峰% > %过滤器(da_padj<0.01）

现在我们有农庄含有DE基因的对象，DE信号不强的基因，DA峰值。我们准备回答这个问题:与没有足够DE信号的基因相比，DE基因附近是否存在DA ATAC-seq峰的富集?

4.2向下采样非差异表达基因

作为plyranges是建在上面的dplyr，它实现了for的许多动词的方法农庄对象。这里我们可以用片的行进行随机抽样other_genes．的sample.int函数将生成大小等于中行数中的de -基因数的随机样本other_genes：

大小< -长度(de_genes)片(other_genessample.int（n()、大小))

与749范围和3 # #农庄对象元数据列:# # seqnames范围链| gene_id de_log2FC # # < Rle > < IRanges > < Rle > | <人物> <数字> # # [1]chr3 129430944 - 129440179 | ENSG00000129071.9 -0.278242 # # [2] chr20 63272785 - 63272785 + | ENSG00000101199.12 0.375979 # # [3] chr19 12224686 - 12294887 | ENSG00000197857.13 0.193328 # # [4] chr6 28224886 - 28224886 + | ENSG00000137185.12 0.116432 # # [5] chr7 157335381 - 157335381 + | ENSG00000105993.14 -0.152200  ## ... ... ... ... . ... ...# # [745] chr1 103525691 - 103525691 + | ENSG00000185946.15 0.161025 # # [746] chr3 11272309 - 11272309 + | ENSG00000197548.12 -0.169212 # # [747] chr1 235327350 - 235327350 + | ENSG00000152904.11 0.315958 # # [748] chr17 76712832 - 76726799 | ENSG00000070495.14 0.253207 # # [749] chr5 178204507 - 178204507 + | ENSG00000197451.11 -0.559576 # # de_padj # # <数字> 1.00756 e 03 # # # # [1] [2] 8.39617 e-06 e-09 # # 5.83820 # #[3][4] 1.19755平台以及# # 3.90094 e-08 [5]  ## ... ...## [745] 2.48758e-13 ## [746] 1.02457e-11 ## [747] 3.22481e-07 ## [748] 2.46375e-10 ## [749] 3.25575e-03 ## ------- seqinfo:来自hg38基因组的25个序列(1个循环)

我们可以重复多次，通过创建多个样本复制．通过多次重复子采样，我们最小化了由采样过程引起的富集统计量的方差。

#为结果设置种子set.seed（2019-08-02）boot_genes < -复制（10，片(other_genessample.int（n(),大小)),简化=假）

的列表农庄对象作为列表，我们可以使用bind_ranges函数。这个函数在结果上创建一个名为“resample”的新列，用于标识每个输入农庄对象:

boot_genes < -bind_ranges(boot_genes.id =“重新取样”）

类似地，我们可以将boot_genes农庄，与DE农庄对象。由于DE上没有重新采样列农庄对象，这是给定的一个缺失的值，我们重新编码为0使用变异()

all_genes < -bind_ranges（德=de_genes,not_de =boot_genes,.id =“起源”）% > %变异（起源=因素(来源、c（“not_de”，“德”)),重新取样=ifelse（is.na(重新取样),0 l,as.integer(重新取样))）all_genes

与8239年# #农庄对象范围和5元数据列:# # seqnames范围链| gene_id de_log2FC # # < Rle > < IRanges > < Rle > | <人物> <数字> # # [1]chr1 196651878 - 196651878 + | ENSG00000000971.15 4.98711 # # [2] chr6 46129993 - 46129993 + | ENSG00000001561.6 1.92722 # # [3] chr4 17577192 - 17577192 + | ENSG00000002549.12 2.93373 # # [4] chr7 150800403 - 150800403 + | ENSG00000002933.8 3.16722 # # [5] chr4 15778275 - 15778275 + | ENSG00000004468.12 5.40894  ## ... ... ... ... . ... ...# # [8235] chr17 43527844 - 43579620 | ENSG00000175832.12 -0.240918 # # [8236] chr17 18260534 - 18260534 + | ENSG00000177427.12 -0.166059 # # [8237] chr20 63895182 - 63895182 + | ENSG00000101152.10 0.250539 # # [8238] chr1 39081316 - 39081316 + | ENSG00000127603.25 -0.385054 # # [8239] chr8 41577187 - 41577187 + | ENSG00000158669.11 0.155922 # # de_padj重新取样起源# # <数字> <整数> <因素> # # 1.37057 [1]e-13 de # # 3.17478 [2] e-05 0 de # # 2.01310 [3] e-11 de # # 1.07360 [4] e-08 0 de # # [5]4.82905e-18 0 de ## ... ... ... ...## [8235] 9.91611e-03 10 not_de ## [8236] 9.12051e-05 10 not_de ## [8237] 1.74085e-09 10 not_de ## [8238] 2.65539e-03 10 not_de ## [8239] 2.96375e-17 10 not_de ## -------来自hg38基因组的25个序列(1个循环)

4.３扩展转录起始位点周围的基因组坐标

现在我们想修改我们的基因范围，使它们在转录起始位点(TSS)的两侧包含10个碱基。有很多方法可以做到这一点，但我们更喜欢通过中的锚定方法plyranges．因为a的起始，结束，宽度和链之间是相互依赖的农庄对象时，我们定义锚来修复其中的一个开始而且结束，同时修改宽度．举个例子，为了只提取TSS，我们可以锚定在范围的5 '端，并将范围的宽度修改为1。

all_genes < -all_genes% > %anchor_5p（）% > %变异（宽度=1）

锚定5 '结束的范围将固定结束负链范围，并固定开始正链范围。

然后我们可以重复同样的模式，但这次使用anchor_center ()告诉plyranges我们将TSS设为总宽度为20kb的范围的中点，或将TSS的上游和下游均设为10kb。

all_genes < -all_genes% > %anchor_center（）% > %变异（宽度=2＊1 e4）

4.4使用重叠连接来查找相对富集

我们现在准备计算RNA-seq基因(我们的DE集和bootstrap集)和ATAC-seq峰值之间的重叠。在plyranges，重叠被定义为两个之间的连接农庄对象:左和一个正确的农庄对象。在重叠连接中，匹配是对象上的任何范围左农庄它与正确的农庄．重叠连接的一个强大的方面是，结果维护每个节点的所有(元数据)列左而且正确的范围，这使得下游摘要很容易计算。

为了将DE基因与DA峰值结合起来，我们执行左重叠连接。这就回到了all_genes范围(可能存在重复)，但是使用来自重叠DA峰值的元数据列。对于任何没有重叠的基因，DA峰列将有重叠NA的年代。

genes_olap_peaks < -all_genes% > %join_overlap_left(da_peaks)genes_olap_peaks

与27766范围和8 # #农庄对象元数据列:# # seqnames范围链| gene_id de_log2FC # # < Rle > < IRanges > < Rle > | <人物> <数字> # # [1]chr1 196641878 - 196641878 + | ENSG00000000971.15 4.98711 # # [2] chr6 46119993 - 46119993 + | ENSG00000001561.6 1.92722 # # [3] chr4 17567192 - 17567192 + | ENSG00000002549.12 2.93373 # # [4] chr4 17567192 - 17567192 + | ENSG00000002549.12 2.93373 # # [5] chr4 17567192 - 17567192 + | ENSG00000002549.12 2.93373  ## ... ... ... ... . ... ...# # [27762] chr1 39071316 - 39071316 + | ENSG00000127603.25 -0.385054 # # [27763] chr1 39071316 - 39071316 + | ENSG00000127603.25 -0.385054 # # [27764] chr8 41567187 - 41567187 + | ENSG00000158669.11 0.155922 # # [27765] chr8 41567187 - 41567187 + | ENSG00000158669.11 0.155922 # # [27766] chr8 41567187 - 41567187 + | ENSG00000158669.11 0.155922 # # de_padj重新取样起源peak_id da_log2FC da_padj # # <数字> <整数> <因素> <人物> <数字> <数字> # # 1.37057 e-13 0 de ATAC_peak_21236 -0.546582 [1][3] 2.01310e-11 0 de ATAC_peak_193578 0.222371 3.00939e-11 ## [4] 2.01310e-11 0 de ATAC_peak_193579 -0.281615 7.99889e-05 ## [5] 2.01310e-11 0 de ATAC_peak_193580 0.673705 7.60043e-15 ## ... ... ... ... ... ... ...## [27762] 2.65539e-03 10 not_de ATAC_peak_5357 -1.058236 3.69052e-16 ## [27763] 2.65539e-03 10 not_de ATAC_peak_5358 -1.314112 6.44280e-26 ## [27764] 2.96375e-17 10 not_de ATAC_peak_263397 0.364738 2.08835e-08 ## [27766] 2.96375e-17 10 not_de ATAC_peak_263399 0.317387 1.20088e-08 ## ------- # seqinfo:来自hg38基因组的25个序列(1个循环)

现在我们可以问，相对于“其他”非DE基因，有多少DA峰在DE基因附近?一个基因可以出现不止一次genes_olap_peaks，因为多个峰值可能会重叠在一个基因上，或者因为我们对同一个基因进行了不止一次的重新采样，或者这两种情况的结合。

对于每个基因(即染色体、起始、结束和链的组合)和“起源”(DE vs非DE)，我们可以计算每个基因的不同峰数和基于LFC的最大峰。这是通过reduce_ranges_directed，它允许聚合结果为农庄对象通过合并邻近的基因组区域。使用有向后缀表示我们在维护链信息。在这种情况下，我们只是通过上面提到的组合并范围(基因)。我们还必须考虑我们在计算是否有任何峰值时执行的重采样数量，以确保我们不会重复计算同一个峰值:

gene_peak_max_lfc < -genes_olap_peaks% > %group_by(gene_id起源)% > %reduce_ranges_directed（peak_count =总和（！is.na(da_padj))/n_distinct(重新取样),peak_max_lfc =马克斯（腹肌(da_log2FC))）

然后，我们可以过滤基因，如果它们有任何峰值，并使用箱线图比较非DE和DE基因之间的峰值折叠变化:

图书馆(ggplot2)gene_peak_max_lfc% > %过滤器(peak_count>0）% > %as.data.frame（）% > %ggplot（aes(起源、peak_max_lfc))+geom_boxplot（）

图4.1:与至少有一个DA峰的非DE基因相比，DE基因DA峰的最大lfc的箱图。

一般而言，DE基因相对于非DE基因具有更大的最大DA折叠变化。

接下来，我们研究DA LFC上的阈值如何修改我们在DE或非DE基因附近观察到的DA峰的富集。首先，我们想知道当我们改变峰值LFC的阈值时，DE基因和非DE基因内的峰值数量是如何变化的。作为一个例子，我们可以通过任意选择LFC阈值1或2来计算，如下所示:

origin_peak_lfc < -genes_olap_peaks% > %group_by(来源)% > %总结（peak_count =总和（！is.na(da_padj))/n_distinct(重新取样),lfc1_peak_count =总和（腹肌(da_log2FC)>1，na.rm =真正的）/n_distinct(重新取样),lfc2_peak_count =总和（腹肌(da_log2FC)>2，na.rm =真正的）/n_distinct(重新取样)）origin_peak_lfc

##数据帧2行4列## origin peak_count lfc1_peak_count lfc2_peak_count ## <因子> <数字> <数字> <数字> ## 1 not_de 2391.8 369.5 32.5 ## 2 de 3416.0 1097.0 234.0

在这里，我们看到DE基因倾向于在它们附近有更多的DA峰，并且随着我们增加DA LFC阈值，DA峰的数量减少(正如预期的那样)。我们现在展示了如何计算DE基因与非DE基因的峰值计数的比例，因此我们可以看到这个比例对于各种DA LFC阈值是如何变化的。

对于所有变量，除了起源列中，我们将第一行的值除以第二行，这将是DE基因相对于其他基因的峰的富集。方法将汇总表从长格式重新塑造为宽格式tidyr包中。首先，我们对结果进行pivotpeak_count列的名称-值对，然后再次旋转以将值放入起源列。然后我们创建一个相对丰富的新列:

origin_peak_lfc% > %as.data.frame（）% > %tidyr：：pivot_longer（关口=-起源)% > %tidyr：：pivot_wider（names_from =起源、values_from =值)% > %变异（浓缩=德/not_de)

## # A tibble: 3 × 4 name not_de de enrichment ##     ## 1 peak_count 2392。3416 1.43 ## 2 lfc1_peak_count 370。3 lfc2_peak_count 32.5 234 7.2

上表显示，LFC阈值越大，相对富集程度越高。

由于基因对峰值的一对多映射，当我们增加DA LFC上的阈值时，我们不知道我们是否有相同数量的DE基因参与或更少。我们可以通过分组和聚合两次来检查在不同阈值处DA峰重叠的基因数量。首先，在每个基因、起源和重采样组中计算满足阈值的峰值数量。其次，在原点列中，我们计算满足DA LFC阈值的峰值总数和具有超过零峰值的基因数量(再次对重采样的数量进行平均)。

genes_olap_peaks% > %group_by(gene_id, origin, resample)% > %reduce_ranges_directed（lfc1 =总和（腹肌(da_log2FC)>1，na.rm =真正的)，lfc2 =总和（腹肌(da_log2FC)>2，na.rm =真正的））% > %group_by(来源)% > %总结（lfc1_gene_count =总和(lfc1>0）/n_distinct(重新取样),lfc1_peak_count =总和(lfc1)/n_distinct(重新取样),lfc2_gene_count =总和(lfc2>0）/n_distinct(重新取样),lfc2_peak_count =总和(lfc2)/n_distinct(重新取样)）

##数据帧与2行5列##起源lfc1_gene_count lfc2_gene_count lfc2_peak_count ## <因子> <数字> <数字> <数字> <数字> ## 1 not_de 271.2 369.5 30.3 32.5 ## 2 de 515.0 1097.0 151.0 234.0

为许多阈值执行此操作非常麻烦，并且会创建大量重复的代码。相反，我们创建了一个名为count_above_threshold中每个元素的变量绝对值的和阈值向量。

count_if_above_threshold < -函数(var，阈值){拉普兰人(阈值,函数(.)总和（腹肌(var)>．，na。rm =真正的）)｝

上述函数将计算任意阈值的计数，因此我们可以将其应用于可能感兴趣的LFC阈值。中LFC绝对值的范围为基础，我们选择了包含100个阈值的网格da_peaks农庄对象:

阈值< -da_peaks% > %变异（abs_lfc =腹肌(da_log2FC))% > %与（seq（最小值(abs_lfc),马克斯(abs_lfc),长度。了=One hundred.））

每个阈值的峰值计数作为一个名为价值．首先,农庄对象按照基因、起源和重采样列的数量进行了分组。然后我们对这些列进行汇总，因此每行将包含一个基因、起源和重新采样的所有阈值的峰值计数。我们还维护了另一个包含阈值的列表列。

genes_peak_all_thresholds < -genes_olap_peaks% > %group_by(gene_id, origin, resample)% > %reduce_ranges_directed（值=count_if_above_threshold(da_log2FC阈值),阈值=列表(阈值)）genes_peak_all_thresholds

与8239年# #农庄对象范围和5元数据列:# # seqnames范围链| gene_id起源# # < Rle > < IRanges > < Rle > | <人物> <因素> # # [1]chr1 196641878 - 196641878 + | ENSG00000000971.15 de # # [2] chr6 46119993 - 46119993 + | ENSG00000001561.6 de # # [3] chr4 17567192 - 17567192 + | ENSG00000002549.12 de # # [4] chr7 150790403 - 150790403 + | ENSG00000002933.8 de # # [5] chr4 15768275 - 15768275 + | ENSG00000004468.12 de  ## ... ... ... ... . ... ...## [8235] chr17 43569620-43589619 - | ENSG00000175832.12 not_de ## [8236] chr17 18250534-18270533 + | ENSG00000177427.12 not_de ## [8237] chr20 63885182-63905181 + | ENSG00000101152.10 not_de ## [8238] chr1 39071316-39091315 + | ENSG00000127603.25 not_de ## [8239] chr8 41567187-41587186 + | ENSG00000158669.11 not_de ## resample value threshold ##    ##[1] 0 1,1,1，…0.0658243, 0.1184840, 0.1711436,……##[2] 0 1,1,1，…0.0658243, 0.1184840, 0.1711436,……##[3] 0 6,6,6，…0.0658243, 0.1184840, 0.1711436,……##[4] 0 4,4,4，…0.0658243, 0.1184840, 0.1711436,……##[5] 0 11,11,11，… 0.0658243,0.1184840,0.1711436,... ## ... ... ... ... ## [8235] 10 1,1,1,... 0.0658243,0.1184840,0.1711436,... ## [8236] 10 3,3,2,... 0.0658243,0.1184840,0.1711436,... ## [8237] 10 5,5,5,... 0.0658243,0.1184840,0.1711436,... ## [8238] 10 3,3,3,... 0.0658243,0.1184840,0.1711436,... ## [8239] 10 3,3,3,... 0.0658243,0.1184840,0.1711436,... ## ------- ## seqinfo: 25 sequences (1 circular) from hg38 genome

现在我们可以将这些列表列展开为一个长农庄对象使用expand_ranges ()函数。该函数将取消列表价值而且阈值列并加长结果农庄对象。为了计算每个阈值的峰值和基因计数，我们应用与之前相同的总结:

origin_peak_all_thresholds < -genes_peak_all_thresholds% > %expand_ranges（）% > %group_by(起源、阈值)% > %总结（gene_count =总和(值>0）/n_distinct(重新取样),peak_count =总和(值)/n_distinct(重新取样)）origin_peak_all_thresholds

## DataFrame 200行4列## origin threshold gene_count peak_count ##     ## 1 not_de 0.0658243 708.0 2391.4 ## 2 not_de 0.1184840 698.8 2320.6 ## 3 not_de 0.1711436 686.2 2178.6 ## 4 not_de 0.2238033 672.4 1989.4 ## 5 not_de 0.2764629 650.4 1785.8 ## ... ... ... ... ...## 196 de 5.06849 22 ## 197 de 5.12115 00 ## 198 de 5.17381 00 ## 199 de 5.22647 00 ## 200 de 5.27913 00

同样，我们可以用与以前相同的方式计算lfc中的相对富集，即将结果旋转到长形式，然后返回到宽形式，以计算富集。我们将DE基因相对于其他基因的峰值富集变化可视化为折线图:

origin_threshold_counts < -origin_peak_all_thresholds% > %as.data.frame（）% > %tidyr：：pivot_longer（关口=-c(起源、阈值),names_to =c（“类型”，“var”)，names_sep =“_”，values_to =“数”）% > %选择（-var)origin_threshold_counts% > %过滤器(类型= =“高峰”）% > %tidyr：：pivot_wider（names_from =起源、values_from =数)% > %变异（浓缩=德/not_de)% > %ggplot（aes（x =阈值,y =浓缩)+geom_line（）+实验室（x =“logFC阈值”，y =“相对浓缩”）

##警告:删除了4行包含缺失值(geom_path)。

一个折线图显示了当DA LFC绝对阈值增加时，DE基因与非DE基因之间DA峰值的相对富集程度如何变化。

图4.2:折线图显示了当DA LFC绝对阈值增加时，DE基因与非DE基因之间DA峰值的相对富集程度的变化。

我们计算了DE基因附近DA峰值的总和，以增加LFC阈值对可及性变化的影响。当我们增加阈值时，总峰数下降(每个基因的DA峰的平均值也是如此)。这也有可能是在DA峰值附近的DE基因数量减少了。我们可以用一个图表来研究这个问题，该图表总结了上述丰富图表背后的许多方面。

origin_threshold_counts% > %ggplot（aes（x =阈值,y =数+1，颜色=起源、线型=类型)+geom_line（）+scale_y_log10（）

$折线图显示了当DA LFC绝对阈值增加时基因和峰值计数的变化。线条的颜色取决于它们是否代表DE基因。注意x轴是在$\log_{10}$刻度上。$

图4.3:折线图显示了当DA LFC绝对阈值增加时基因和峰值计数的变化。线条的颜色取决于它们是否代表DE基因。注意x轴在a上\ (\ log_ {10} \)规模。

5讨论

我们已经证明了plyranges而且tximeta(在Bioconductor和tidyverse生态系统)，我们可以流畅地迭代生物数据科学工作流程:从导入，到建模，再到数据集成。

在这个分析中，还有几个进一步的步骤会很有趣;例如，我们可以修改TSS周围的窗口大小，以了解它如何影响富集，并改变DE基因和DA峰值集的FDR截断值。除了自举集的均值之外，我们还可以计算方差，因此在富集线周围画一个区间。

最后，我们的工作流程说明了使用Bioconductor提供的适当数据抽象的好处，例如SummarizedExperiment而且农庄．这些抽象为用户提供了实验数据的心理模型，是构建我们在这里展示的模块化和迭代分析的构建模块。因此，我们已经能够互操作许多分离的R包(来自Bioconductor和tidyverse)，以构建一个无缝的端到端工作流，这对于单个单片工具来说太专业化了。

6软件可用性

工作流材料，包括本文，可以完全按照在Github存储库中找到的说明复制sa-lee / fluentGenomics．此外，工作流的开发版本和所有下游依赖项都可以使用BiocManager通过运行:

#开发版本从GithubBiocManager：：安装（“sa-lee / fluentGenomics”）#版本可从Bioconductor获得BiocManager：：安装（“fluentGenomics”）

本文和分析是用R(R核心团队2019）使用rmarkdown(Allaire等。2019）,knitr(谢2019，2015）包。

6.1会话信息

sessionInfo（）

## 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。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## 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_phone = 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 plyranges_1.17.0 ## [3] DESeq2_1.37.7 GenomicFeatures_1.49.8 ## [7] Biobase_2.57.3 genome ranges_1. 49.1 ## [9] GenomeInfoDb_1.33.17 IRanges_2.31.2 ## [11] S4Vectors_0.35.4 BiocGenerics_0.43.4 ## [13] MatrixGenerics_1.9.1 matrixStats_0.62.0 ## [15] readr_2.1.3 dplyic_1 .0.10 ## [17] tximeta_1.15.4 fluentGenomics_1.11.0 ## ##通过命名空间加载(并且没有附加):# # # # [1] colorspace_2.0-3 rjson_0.2.21 [3] ellipsis_0.3.2 XVector_0.37.1 # # [5] farver_2.1.1 bit64_4.0.5 # # [7] interactiveDisplayBase_1.35.1 fansi_1.0.3 # # [9] xml2_1.3.3 codetools_0.2-18 # # [11] splines_4.2.1 tximport_1.25.1 # # [13] cachem_1.0.6 geneplotter_1.75.0 # # [15] knitr_1.40 jsonlite_1.8.3 # # [17] Rsamtools_2.13.4 annotate_1.75.0 # # [19] dbplyr_2.2.1 png_0.1-7 # # [21] shiny_1.7.3 BiocManager_1.30.19 # # [23] compiler_4.2.1 httr_1.4.4 # # [25] assertthat_0.2.1 Matrix_1.5-1 # # [27][39] jquerylib_0.1.4 vctrs_0.5.0 ## [41] Biostrings_2.65.6 rtracklayer_1.57.0 ## [43] xfun_0.34 string_1 .4.1 ## [45] mime_0.12 lifecycle_1.0.3 ## [47] restfulr_0.0.15 ensembldb_2.21.5 ## [49] XML_3.99-0.12 AnnotationHub_3.5.2 ## [51] scales_1.2.1 zlibbioc_1.43.0 ## [53] vroom_1.6.0 ## [35] glue_1.6.2 GenomeInfoDbData_1.2.9 ##hms_1.1.2 ## [55] promises_1.2.0.1 ProtGenerics_1.29.1 ## [57] parallel_4.2.1 RColorBrewer_1.1-3 ## [59] AnnotationFilter_1.21.0 yaml_2.3.6 ## [61] curl_4.3.3 memoise_2.0.1 ## [63] sass_0.4.2 biomaRt_2.53.5 ## [65] stringi_1.7.8 RSQLite_2.2.18 ## [67] BiocVersion_3.16.0 highr_0.9 ## [69] genefilter_1.79.0 BiocIO_1.7.1 ## [71] filelock_1.0.2 BiocParallel_1.31.15 ## [75] bitops_1.0-7 evaluate_0.17 ## [77] lattice_0.20-45 archive_1.1.5 ## [79] purrr_0.3.5labeling_0.4.2 # # [81] GenomicAlignments_1.33.1 bit_4.0.4 # # [83] tidyselect_1.2.0 magrittr_2.0.3 # # [85] bookdown_0.29 R6_2.5.1 # # [87] generics_0.1.3 DelayedArray_0.23.2 # # [89] DBI_1.1.3 pillar_1.8.1 # # [91] withr_2.5.0 survival_3.4-0 # # [93] KEGGREST_1.37.3 rcurl_1.98 - 1.9 # # [95] tibble_3.1.8 crayon_1.5.2 # # [97] utf8_1.2.2 BiocFileCache_2.5.2 # # [99] tzdb_0.3.0 rmarkdown_2.17 # # [101] progress_1.2.2 locfit_1.5 - 9.6 # # [103] grid_4.2.1 blob_1.2.3 # # [105] digest_0.6.30 xtable_1.8-4 # # [107] tidyr_1.2.1 httpuv_1.6.6 ## [109] munsell_0.5.0 bslib_0.4.0

6.2作者的贡献

所有作者都为工作流的编写和开发做出了贡献。

6.3相互竞争的利益

作者宣称他们之间没有利益冲突。

6.4资金

SL由澳大利亚政府研究培训计划(RTP)奖学金和CSL有限公司的顶级奖学金支持。

MIL的贡献由NIH拨款R01 HG009937支持。

我确认资助者在研究设计、数据收集和分析、发表决定或手稿准备中没有任何作用。

6.5确认

作者要感谢Bioconductor 2019和BiocAsia 2019会议的所有参与者，他们参加了本工作流程论文的早期版本并提供了反馈。

参考文献

Alasoo, Kaur, Daniel Gaffney, 2017。“处理巨噬细胞RNA-seq和ATAC-seq实验的读取计数。”Zenodo。https://doi.org/10.5281/zenodo.1188300．

Alasoo, K, J Rodrigues, S Mukhopadhyay, AJ Knights, AL Mann, K Kundu, HIPSCI-Consortium, C Hale, Dougan G和DJ Gaffney。2018。“染色质和基因表达的共同遗传效应表明增强子启动在免疫反应中发挥作用。”自然遗传学50: 424 - 31所示。https://doi.org/10.1038/s41588-018-0046-7．

Allaire, JJ，谢义辉，Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang和Richard Iannone。Rmarkdown: R动态文档．https://github.com/rstudio/rmarkdown．

Frankish, A, gencode - consortium, P Flicek, 2018。"人类和小鼠基因组的GENCODE参考注释"核酸研究．https://doi.org/10.1093/nar/gky955．

吴志金，吴志金。2012。“使用条件分位数归一化去除RNA-seq数据中的技术可变性。”生物统计学13(2): 204-16。https://doi.org/10.1093/biostatistics/kxr054．

Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo，等。2015。“利用Bioconductor组织高通量基因组分析。”自然方法12(2): 115-21。https://doi.org/10.1038/nmeth.3252．

Köster，约翰内斯和斯文·拉赫曼，2012。“snakmake -一个可扩展的生物信息学工作流引擎。”生物信息学28(19): 2520-2。https://doi.org/10.1093/bioinformatics/bts480．

劳，W慈善，Alhamdoosh妈妈，苏珊安，董学艺，田鲁一，戈登·K·史密斯，马修·E·里奇。2018。“使用Limma、Glimma和edgeR, rna序列分析就像1-2-3一样简单。”F1000研究5(1408): 1408。https://doi.org/10.12688/f1000research.9005.3．

劳伦斯，迈克尔，沃尔夫冈·胡贝尔，Hervé Pagès，帕特里克·阿博尤恩，马克·卡尔森，罗伯特·绅士，马丁·摩根和文森特·凯里。2013。计算和注释基因组范围的软件公共科学图书馆第一版。医学杂志。9.https://doi.org/10.1371/journal.pcbi.1003118．

李、斯图尔特、黛安·库克和迈克尔·劳伦斯，2019年。《Plyranges:基因组数据转换的语法》基因组生物学20(1): 4。https://doi.org/10.1186/s13059-018-1597-8．

爱你的，迈克尔一世，西蒙·安德斯，弗拉迪斯拉夫·金，沃尔夫冈·胡贝尔，2016。RNA-Seq工作流程:基因水平的探索性分析和差异表达F1000研究4(1070): 1070。https://doi.org/10.12688/f1000research.7035.2．

《爱》，迈克尔一世，沃尔夫冈·胡贝尔，西蒙·安德斯，2014。“用DESeq2对RNA-seq数据的折叠变化和离散度进行适度估计。”基因组生物学15(12): 550。https://doi.org/10.1186/s13059-014-0550-8．

爱你的迈克尔，夏洛特·森森，彼得·F·希基，丽莎·K·约翰逊，泰莎·皮尔斯，洛里·谢博德，马丁·摩根和罗伯·帕特罗。Tximeta: RNA-seq中来源识别的参考序列校验和bioRxiv九月，777888。https://doi.org/10.1101/777888．

帕特罗，R, G Duggal, MI Love, RA Irizarry和C Kingsford, 2017。“鲑鱼提供了快速和偏见感知的转录表达量化。”自然方法14: 417 - 19所示。https://doi.org/10.1038/nmeth.4197．

R核心团队，2019年。R:统计计算的语言和环境．维也纳，奥地利:R统计计算基金会。https://www.R-project.org/．

谢博德、洛莉和马丁·摩根，2019年。BiocFileCache:跨会话管理文件．

戈登·史密斯，2004。微阵列实验中评估差异表达的线性模型和经验贝叶斯方法遗传学和分子生物学中的统计学应用3(1)。

Charlotte Soneson, Michael I. Love, Mark Robinson, 2015。“RNA-seq的差异分析:转录水平的估计提高了基因水平的推断。”F1000Research4(1521)。https://doi.org/10.12688/f1000research.7563.1．

维克汉姆，哈德利，玛拉·阿弗里克，詹妮弗·布莱恩，温斯顿·张，露西·达戈斯蒂诺·麦高恩，罗曼François，加勒特·格罗蒙德等。2019。“欢迎来到整洁宇宙。”开源软件杂志4(43): 1686。https://doi.org/10.21105/joss.01686．

谢一辉，2015。动态文档与R和针织．第二版。博卡拉顿，佛罗里达州:查普曼;大厅/ CRC。https://yihui.name/knitr/．

——2019.Knitr: R中动态报告生成的通用包．https://yihui.name/knitr/．

精通plyranges和tximeta基因组学

斯图尔特·李，迈克尔·劳伦斯，迈克尔·勒夫

1简介