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基础设施(劳伦斯et al。2013；Huber et al。2015）和tidyverse设计原则维克汉姆等人(2019）．

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

1．1实验数据

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

的巨噬细胞包包含24个RNA-seq样本的RNA-seq量化，RNA-seq样本的子集由Alasoo等人(2018）．配对端读数用大马哈鱼(Patro et al。2017），使用Gencode 29人类参考转录本(法兰克语，gencode - consortium和flick2018）．有关量化的更多细节和使用的确切代码，请参阅巨噬细胞包中。包还包含Snakemake文件，用于分发大马哈鱼集群上的量化作业(Koster和Rahmann2012）．

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

下面的代码将加载到缓存数据文件的路径，如果它不存在，则创建缓存并生成一个SummarizedExperiment使用的BiocFileCache包(牧羊人和摩根2019）．

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

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

对我们如何得到它的精确描述SummarizedExperiment对象可在节中找到2．2．

2导入数据作为SummarizedExperiment

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

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

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

路径相对于当前的R会话。但是，在本例中，我们将文件分布在巨噬细胞包中。可以使用。来定位相关目录和相关文件执行．

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

关于实验的信息包含在coldata.csv文件。我们利用dplyr而且readr包(作为tidyverse)将该文件读入R(韦翰et al。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):复制，ng_ul_mean, rna_auto ## ##ℹ使用' spec() '检索该数据的完整列规范。##ℹ指定列类型或设置' show_col_types = FALSE '来关闭此消息。

coldata

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

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

现在我们有一个总结实验设计和量化位置的表格。以下代码行为分析人员做了大量工作:导入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-04-21 # #发现匹配TxDb通过“AnnotationHub”从缓存加载# # # #加载所需的包:GenomicFeatures # #加载所需的包:AnnotationDbi # # # #附加包:“AnnotationDbi”掩盖了# # # #以下对象从“包:dplyr”:## ##选择## ##生成抄本范围

se

##类:rangedsummarizeexperiment ## dim: 205870 24 ##元数据(6):tximetaInfo quantInfo…txomeInfo txdbInfo ## assays(3): counts abundance length ## 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和Robinson (2015）．

gse < -summarizeToGene(se)

##加载现有的TxDb已创建:2012-04-29 12:03:19

从数据库获取转录到基因的映射

##生成基因范围

# #总结丰富

# #总结计数

# #总结长度

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

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

的SummarizedExperiment对象中包含ATAC-seq峰值，可以从以下制表符分隔的文件中创建Alasoo和Gaffney (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(汉森，伊里扎里，吴2012）和log2转变。加载cqn-规范化矩阵的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使用以下代码块(爱，休伯和安德斯2014）．设计公式表明我们希望控制供体基线(行)，并检测在这种情况下基因表达的差异。有关Bioconductor中DE工作流的更全面的讨论，请参见爱等。2016）而且法等。2018）．

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

##使用tximeta的计数和平均文本长度

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

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

dds < -DESeq(dds)

##估计大小因素

##使用'avgTxLength' from assays(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]see 'cooksCutoff'参数的结果## [2]see 'independentFiltering'参数的结果

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

图3.1:可视化DESeq2结果是一个“MA图”。有调整过的基因假定值0.01以下为红色。

现在我们将结果输出为农庄对象，并由于约定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]0.0859708 0.00000 1.0000000 1.000000 ## [3] 1.2846179 0.1969067 1.44545 0.1483329 1.000000 ## [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 ## [17803] -0.0339296 0.188005 0 1 1 ## [17804] 0.3123477 0.522700 0 1 1 ## [17805] 0.9945187 1.582373 0 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)山峰

## seqnames ranges string | peak_id da_log2FC ##    |   ## [1] chr1 9979-10668 * | ATAC_peak_1 0.266185 ## [2] chr1 10939-11473 * | ATAC_peak_2 0.322177 ## [3] chr1 155039 -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包(劳伦斯et al。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]chr7 5898725 - 5898725 + | ENSG00000122674.11 0.2031642 # # [2] chr12 22625075 - 22625075 + | ENSG00000139163.15 0.2575908 # # [3] chr1 43974487 - 43974487 + | ENSG00000117410.13 0.1148944 # # [4] chr19 55635016 - 55635016 + | ENSG00000213015.8 -0.0892649 # # [5] chr16 21402237 - 21448567 | ENSG00000169246.16 0.1878874  ## ... ... ... ... . ... ...# # [745] chrX 15319451 - 15335580 | ENSG00000165195.15 0.22848693 # # [746] chr5 134194334 - 134226142 | ENSG00000113575.9 0.06052272 # # [747] chr12 990509 - 1495933 + | ENSG00000082805.19 0.00646525 # # [748] chr11 61130256 - 61161617 | ENSG00000167987.10 -0.09338839 # # [749] chr17 17810399 - 17837002 | ENSG00000072310.16 0.30044350 # # de_padj # # <数字> 7.06653 e-07 # # # # [1] [2] 4.69348 e-06 e-13 # # 4.53819 # # [3] [4] 2.61090 e-08 # # 6.06624 e-07 [5]  ## ... ...## [745] 2.09461e-07 ## [746] 8.24597e-24 ## [747] 8.07219e-09 ## [748] 1.09401e-12 ## [749] 2.45385e-04 ## ------- ## 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 ## ------- ## seqinfo:来自hg38基因组的25个序列(1个圆形)

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

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

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

在距离的5英尺处锚定会固定结束的负链范围，并固定开始正链范围。

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

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

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

我们现在准备计算RNA-seq基因(我们的DE集和引导集)和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]1.15274e-04 ## [2] 3.17478e-05 0 de ATAC_peak_231183 1.453297 9.73225e-17 ## [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] .65539e-03 10 not_de ATAC_peak_5357 -1.058236 3.69052e-16 ## [27763] .65539e-03 10 not_de ATAC_peak_5358 -1.314112 6.44280e-26 ## [27764] 2.96375e-17 10 not_de ATAC_peak_263396 -0.904080 8.19577e-13 ## [27765] 2.96375e-17 10 not_de ATAC_peak_263397 0.364738 2.08835e-08 ## [27766] 2.96375e-17 10 not_de ATAC_peak_263397 0.317387 1.20088e-08 ## ------- ## seqinfo: 25个序列(1个循环)来自hg38基因组

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

对于每个基因(即染色体、开始、结束和链的组合)和“起源”(DE vs not-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包中。首先我们对结果进行透视peak_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。1097 2.97 ## 3 lfc2_peak_count 32.5 234 7.2

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

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

genes_olap_peaks% > %group_by(gene_id、起源、重新取样)% > %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(重新取样)）

## DataFrame with 2行5列## origin lfc1_gene_count lfc2_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、起源、重新取样)% > %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 ##重采样值阈值##    ##[1] 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

## 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

同样，我们可以用同样的方式计算低碳碳化合物中的相对富集，即将结果旋转到长形式，然后再返回到宽形式来计算富集。我们将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”）#版本可从BioconductorBiocManager：：安装（“fluentGenomics”）

本文采用R(R核心团队2019）使用rmarkdown(阿莱尔et al。2019）,knitr(谢2019，2015）包。

6.1会话信息

sessionInfo（）

## R版本4.2.0 RC (22-04-19 r82224) ##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.4 LTS ## ##矩阵产品:default ## BLAS: /home/biocbuild/bbs-3.15-bio /R/lib/libRblas. ##因此## LAPACK: /home/biocbuild/bbs-3.15-bio /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_TELEPHONE= c# [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ## [8] base ## ##其他附加的包:# [1] ggplot_2 .3.5 plyranges_1.16.0 ## [3] DESeq2_1.36.0 genomic icfeatures_1 .48.0 ## [5] AnnotationDbi_1.58.0 summarizedexperiment_1 .48.0 ## [7] Biobase_2.56.0 genomic icranges_1 .48.0 ## [9] GenomeInfoDb_1.32.1 IRanges_2.30.0 ## [11] S4Vectors_0.34.0 biocgenerics_0.2.0 ## [13] MatrixGenerics_1.8.0 matrixStats_0.62.0 ## [15] readr_2.1.2 dplyr_1.0.9 ## [17] tximeta_1.14.0 fluentgenomics ics_1.8.0 ## ##通过命名空间加载(并没有连接):# [9] xml2_1. 24.0 splines_4.2.0 ## [11] tximport_1.24.0 cachem_1.0.6 ## [13] geneplotter_1.74.0 knitr_1.39 ## [15] jsonlite_1.8.0 Rsamtools_2.12.0 ## [17] annotate_1.74.0 dbplyr_2.1.1 ## [19] png_0.1-7 shiny_1.7.1 ## [21] BiocManager_1.30.17 compiler_4.2.0 ## [23] httr_1.4.2 assertthat_0.2.1 ## [25] Matrix_1.4-1 fastmap_1.1.0 ## [27]## [35] GenomeInfoDbData_1.2.8 rappdirs_0.3.3 ## [37] Rcpp_1.0.8.3 jquerylib_0.1.4 ## [39] vctrs_0.4.1 Biostrings_2.64.0 ## [41] rtracklayer_1.56.0 xfun_0.30 ## [43] string_1 .4.0 mime_0.12 ## [45] lifecycle_1.0.1 restfulr_0.0.13 ## [47] integrationhub_2 .20.0 xml_3 .4.0 scales_1.2.0 ## [49] zlibbioc_1.42.0 vroom_1.5.7 ## [53] hms_1.1.1 . ## [35] genomeinfodbdata_1.1.1 ## [39] vctrs_0.4.1promises_1.2.0.1 # # [55] ProtGenerics_1.28.0 parallel_4.2.0 # # [57] RColorBrewer_1.1-3 AnnotationFilter_1.20.0 # # [59] yaml_2.3.5 curl_4.3.2 # # [61] memoise_2.0.1 sass_0.4.1 # # [63] biomaRt_2.52.0 stringi_1.7.6 # # [65] RSQLite_2.2.12 BiocVersion_3.15.2 # # [67] highr_0.9 genefilter_1.78.0 # # [69] BiocIO_1.6.0 filelock_1.0.2 # # [71] BiocParallel_1.30.0 rlang_1.0.2 # # [73] pkgconfig_2.0.3 bitops_1.0-7 # # [75] evaluate_0.15 lattice_0.20-45 # # [77] purrr_0.3.4 labeling_0.4.2 # # [79]基因组icalignments_1 .32.0 bit_4.0.4 ## [81] tidyselect_1.1.2 magrittr_2.0.3 ## [83] bookdown_0.26 R6_2.5.1 ## [85] generics_0.26 r6_2 .22.0 ## [87] DBI_1.1.2 pillar_1.7.0 ## [89] withr_2.5.0 survival_1 .3-1 ## [91] KEGGREST_1.36.0 rcurl_1 .3.0 ## [95] utf8_1.2.2 biocfilecache_1 .4.0 ## [97] tzdb_0.3.0 rmarkdown_2.14 ## [99] progress_1.2.2 locfit_1.5-9.5 ## [101] grid_4.2.0 blob_1.2.3 ## [105] tidyr_1.2.0 httpuv_1.6.5 ## [107] munsell_0.5.0 bslib_0.3.1

6.2作者的贡献

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

6.3相互竞争的利益

作者声明他们没有竞争利益。

6.4资金

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

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

我确认资助方在研究设计、数据收集和分析、出版决定或稿件准备方面没有任何作用。

6.5确认

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

参考文献

考尔和丹尼尔·加夫尼，2017。“巨噬细胞RNA-seq和ATAC-seq实验的读取计数处理。”Zenodo。https://doi.org/10.5281/zenodo.1188300．

阿拉索，K, J罗德里格斯，S Mukhopadhyay, AJ Knights, AL Mann, K Kundu, hipscii - consortium, C Hale, Dougan G和DJ Gaffney。2018。“染色质和基因表达的共同遗传效应表明增强子启动在免疫反应中起着作用。”自然遗传学50: 424 - 31所示。https://doi.org/10.1038/s41588-018-0046-7．

阿莱尔、JJ、谢一辉、乔纳森·麦克弗森、哈维尔·卢拉西、凯文·乌西、阿伦·阿特金斯、哈德利·威克姆、乔·程、温斯顿·张和理查德·伊诺内。2019。Rmarkdown: R的动态文档．https://github.com/rstudio/rmarkdown．

francish, 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, Johannes和Sven Rahmann. 2012。“snakmake -一个可扩展的生物信息学工作流引擎。”生物信息学28日(19):2520 - 2。https://doi.org/10.1093/bioinformatics/bts480．

Law, Charity W, Monther Alhamdoosh, Su shan，董学义，田鲁一，Gordon K Smyth, Matthew E Ritchie. 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·约翰逊，N特莎·皮尔斯，洛莉·谢博德，马丁·摩根和罗布·帕特罗。2019。Tximeta: RNA-seq中物源识别的参考序列校验和。bioRxiv, 777888年9月。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和Knitr．第二版。博卡拉顿，佛罗里达州:查普曼;大厅/ CRC。https://yihui.name/knitr/．

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

流利的基因组学与plyrange和tximeta

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

1简介