# 基于 R 的计算基因组学 本书为 [Computational Genomics with R](http://compgenomr.github.io/book/) 的中文翻译版本，目前正在稳定更新中。 **本书固定访问地址**： ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200710222705.png) ## 中文译者序 初次接触到这本书是在 2019 年底，那时读了几章感觉整体写的不错，很适合 R 语言大致入门后的进阶学习。当时提交了几个原书中可以改进的地方，作者一直迟迟没有反应，在我以为要弃更的时候，2 个月后突然合并了。 本来是想着把整个书仔细过一篇，顺便再提交一些改进建议，然而因为各种事情太多就放下了，以至于我都删除了自己 folk 的版本。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200710223149.png) 直到最近再看这本书的时候，发现作者在致谢部分竟然也把我的名字写上去了，瞬间有点惭愧。 在和果子聊天的过程中，他跟我提到自己也把这本书推荐给了线下课程的学员，但是真正看了的人可能寥寥无几。于是我就动了出一个持续更新的中文翻译版的念想，也算是值得被在致谢中提到一下。 目前计划每周更新一部分，翻译到中后期再出一个完整的中文本链接在网上发布供大家集中阅读。如果你的英文读写还过得去，那么一可以考虑直接阅读这本书的原版内容，二可以考虑和我联系一起进行翻译。 在翻译的过程中，我会尽量忠实于作者的原作，如果有些地方我认为有更好的实现方法，会作为补充放在正文后供你参考。如果你正在找这样一本书，希望它对你有所帮助，如果你偷懒想看看中文的，希望这个翻译版对你有所帮助。 [思考问题的熊](https://kaopubear.top) 2020 年 7 月 # 前言 ## 前言 本书的目的是为基因组学的数据分析提供基础知识。我们根据每年提供的计算基因组学课程开发了本书，受众来自物理、生物、医学、数学、计算机科学或其他领域的交叉学科。我们希望这本书可以成为计算基因组学专业学生的起点，也可以成为基因组学中各具体领域进一步数据分析的指南。这就是为什么我们试图涵盖从编程到基本基因组生物学的各种主题。由于这个领域是跨学科的，因此不同背景的人可以从不同的地方开始。生物学家可能会跳过基因组生物学的部分直接从 R 编程开始，而计算机科学家可能会想从基因组生物学开始。同样，当你需要在没有经验的情况下进行某种类型的分析时，也可能会想参考这本书。 ## 这本书是给谁的？ 这本书包含计算基因组学的实践和理论等方面。生物学和医学产生的数据比以往任何时候都多。因此，我们需要更多同时具有数据分析技能并理解计算基因组学的人。由于计算基因组学是跨学科的，这本书旨在为生物学家，医学科学家，计算机科学家和来自其他背景的人提供一些帮助。本书的受众包括: * 产生数据并热衷于自己分析数据的生物学家和医学科学家。 * 正式开始研究或使用计算基因组学的学生和研究人员，虽然他们没有广泛的领域特定知识，但至少对定量领域有初级理解: 数学、统计 * 经验丰富的研究人员正在寻找快速操作方法，以便开始与计算基因组学相关的特定数据分析任务。 ## 你能从中得到什么？ 该资源描述了这些技能，并提供了帮助读者分析自己基因组数据的方法。 阅读后: * 如果你不熟悉 R，你将会掌握 R 的基础知识，并深入了解 R 在计算基因组学中的特殊用途。 * 你将了解基因组间隔和如何对它们的操作，例如overlap * 你将能够使用 R 及其庞大的包进行序列分析: 例如计算基因组给定片段的 GC 含量或找到转录因子结合位点 * 你将熟悉基因组学中使用的可视化技术，如热图、元基因图和基因组track可视化 * 你将熟悉监督和无监督学习技术，这些技术对数据建模和高维数据的探索性分析很重要 * 你将熟悉不同的高通量测序数据集的分析，主要使用基于 R 的工具。 ## 本书结构 这本书的设计理念是数据分析方法和概念理论的理解同样重要。这也是为什么我们首先尝试给出概念的解释，然后尝试给出数学公式的基本内容以获得更详细的理解。我们也会显示代码并对特定数据分析任务的代码进行解释。此外，我们还为那些希望对数据分析相关方法或概念有更深入理论理解的读者提供额外的参考，如书籍、网站、视频讲座和科学论文。 「基因组学导论」一章介绍了基因组生物学和基因组学的基本概念。理解这些概念对计算基因组学很重要。 「基因组数据分析的 R 介绍」一章除了提供我们在基因组数据分析中观察到的常见数据分析格式外，也包括本书所必需的基本 R 技能。基因组统计，无监督机器学习的探索性数据分析和监督机器学习的预测建模章节介绍了分析高维基因组数据时可能需要的必要技能。 「基因组区间操作和基因组计算」一章介绍了处理基因组区间的基本工具以及它们在基因组上彼此之间的关系。此外，本章还介绍了用于基因组数据可视化的多种方法。本章中介绍的技能是处理下游基因组数据所需的关键技能，这些数据可通过公共数据库 (如 Ensembl 和 UCSC 浏览器) 获得。 其它章节涉及高通量测序数据的具体分析和不同类型的数据集。「高通量测序数据的质检处理和比对」一章介绍了需要对测序数据进行的质检以及进一步处理它们的不同方法。「RNA-seq 分析」、「ChIP-seq 分析」和「BS-seq 分析」章节涉及流行的高通量测序技术的分析技术。最后一章「多组学分析」讨论了整合多组学数据集的方法。 总而言之，这本书是计算基因组学的综合指南。有些部分是为了广泛的跨学科受众和完整性，并不是所有部分都会对这些广泛受众的所有读者同样有用。 ## 软件信息和约定 包名称以粗体文本 (例如，**methylKit**)表示，行内代码和文件名以打字机字体格式化。函数名后跟括号 (例如，`genomation::ScoreMatrix()`)。双冒号运算符 `::` 表示从包访问对象。 ### 赋值运算符约定 传统上， `<-` 是首选赋值运算符。然而，在整本书中，我们交替使用 `=` 和 `<-` 作为赋值运算符。 ### 运行书籍代码所需的包 这本书主要是关于使用 R 包来分析基因组数据，因此如果你想重现本书中的分析，你需要使用安装在每章中安装相关的包`install.packages`或 `BiocManager::install` 功能。在每章中，当我们使用各个包中所需的函数时，我们用 `library ()` 或 `require()`函数加载必要的包。通过查看这些调用，你可以看到该代码块或章节需要哪些包。如果你需要安装本书的所有包依赖项，你可以运行以下命令，在等待时喝杯茶。 ```r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c('qvalue','plot3D','ggplot2','pheatmap', 'cluster', 'NbClust', 'fastICA', 'NMF', 'Rtsne', 'mosaic', 'knitr', 'genomation', 'ggbio', 'Gviz', 'DESeq2', 'RUVSeq', 'gProfileR', 'ggfortify', 'corrplot', 'gage', 'EDASeq', 'citr', 'formatR', 'svglite', 'Rqc', 'ShortRead', 'QuasR', 'methylKit','FactoMineR', 'iClusterPlus', 'enrichR','caret','xgboost','glmnet', 'DALEX','kernlab','pROC','nnet','RANN', 'ranger','GenomeInfoDb', 'GenomicRanges', 'GenomicAlignments', 'ComplexHeatmap', 'circlize', 'rtracklayer', 'BSgenome.Hsapiens.UCSC.hg38', 'tidyr', 'AnnotationHub', 'GenomicFeatures', 'normr', 'MotifDb', 'TFBSTools', 'motifRG', 'JASPAR2018' )) ``` ## 书内数据 我们依靠来自不同的 R 和bioonducter 包的数据。对于那些不附带这些包的数据集，我们创建了自己的包组合数据。你可以通过 ```r devtools::install_github("compgenomr/compGenomRData") ``` 我们使用 `system.file()` 函数来获取文件的路径。我们注意到许多没有经验的用户对这个功能感到困惑。此函数只输出与数据包一起安装的文件的完整路径。 # 第一章：基因组学简介 更新中 本章的目的是为读者提供理解基因组学所需的一些基础知识。需要说明，这绝不是对这一学科的完整概述，而只是一个简单的总结，它将帮助非生物学相关专业的读者理解计算基因组学中反复出现的生物学概念。熟知基因组生物学和全基因组定量分析的读者可以自由跳过这一章或大致浏览一遍。 # 1.1 基因 DNA 和中心法则 本书一个将会反复出现的核心概念是「基因」。在解释基因之前我们需要先介绍其它几个对理解「基因」很重要的概念。人体由数十亿个细胞组成，这些细胞各自有不同功能。例如，在肝脏中一些细胞能够产生分解毒素的酶，在心脏中有专门的肌肉细胞使心脏跳动。这些不同种类的细胞都来自一个单细胞胚胎，所有制造不同种类细胞的指令都包含在这个细胞中，随着细胞每一次分裂指令都会被传送到新的细胞中。这些指令可以编成 DNA 分子，一种由反复出现的核苷酸组成的聚合物。DNA 分子中的四种核苷酸，腺嘌呤、鸟嘌呤、胞嘧啶和胸腺嘧啶（A、C、G 和 T）以特定的序列储存着生命信息，DNA 两个互补链以双螺旋的形式组织起来。 ## 1.1.1 什么是基因组 一个生物体的完整 DNA 序列包含了所有遗传信息，称为基因组。基因组包含了构建和维持生物体的所有信息，其有不同的大小和结构。我们的基因组不仅仅是一段裸露的 DNA，在真核细胞中，DNA 被蛋白质（组蛋白）包裹形成核小体这样的高级结构，组成染色质和染色体（见图 1.1）。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200711214635.png) 图 1.1: 动物的染色体结构 不同的生物体有不同数量的染色体，在某些物种（如大多数原核生物）中，DNA 是以环状形式储存的。不同物种的基因组大小和染色体数量也是不同，人类的基因组有 46 条染色体，超过 30 亿个碱基对，而小麦的基因组有 42 条染色体，170 亿个碱基对。生物的基因组序列可以通过测序技术获得，通过测序获得的基因组中 DNA 序列片段称为读数(read)。利用彼此重叠的读数将最原始的 DNA 片段拼接成较大的片段，从而获得更长的基因组序列。最新的测序技术使基因组测序变得更便宜且用时更短，输出更多、更长和更准确的读数。 1999-2000 年第一个人类基因组的预估成本是 3 亿美元，如今一个高质量的人类基因组用 1500 美元就可以获得。由于成本下降，研究人员和临床医生可以产生更多的数据。这推动了数据存储成本的上升，也推动了对分析基因组数据人才的需求。当然，这也是写这本书的动机之一。 ## 1.1.2 什么是基因 在基因组中，有一些特定区域包含着为遗传信息物理产物编码的精确信息，基因组中含有这种信息的区域就是我们所谓的「基因」。不过基因的准确定义仍在发展之中，根据经典的分子生物学教材，基因是对应于单个蛋白质或单个催化和结构 RNA 分子的一段 DNA 序列（Alberts 等，2002）。关于基因的一个现代定义是：包括编码一个功能转录本所必需的所有序列元素的区域 (Eilbeck 等人，2005)。不过无论定义如何变化，大家认同基因是所有生物体遗传的基本单位。 所有细胞在大多数时候都在以相同方式使用它们的遗传信息；DNA 的复制则是为了将信息传递给新的细胞。基因被激活后会在细胞核中（真核生物）转录成信使 RNA（mRNA），随后 mRNA（如果基因是蛋白质编码的）在细胞质中被翻译成蛋白质。这实质上是遗传信息在 DNA、RNA 和蛋白质之间的传递过程，该过程也被称为分子生物学的「中心法则」（总结见图 1.2）。蛋白质是生命的基本元素。所有活细胞的生长、修复、功能和结构都依赖于它们。 这就是为什么基因是基因组生物学的核心概念，因为一个基因可以编码蛋白质和其他功能分子的信息。基因如何被控制和激活决定了生物体的一切，从细胞特征到免疫反应，细胞发育和对某些刺激的行为都是由基因和其编码的功能分子活性所决定的。肝细胞之所以成为肝细胞，就是因为某些基因被激活并产生对应的功能产物，从而帮助肝细胞完成任务。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200711214644.png) 图 1.2：中心法则 复制转录翻译 ## 1.1.3 基因如何被控制？转录和转录后调控 为了回答这个问题，我们必须对「中心法则」引入的转录概念进行深入挖掘。信息传递过程中的第一步--DNA 到 RNA--称为转录，由 RNA 聚合酶完成。RNA 聚合酶依赖的转录起始是因为 DNA 序列中一个特定区域--核心启动子的存在而实现。核心启动子是 DNA 序列中可以促进转录的区域，位于转录起始位点上游。在真核生物中，被称为通用转录因子的蛋白质能识别并与核心启动子结合进而形成一个转录起始前复合物。RNA 聚合酶识别这些复合物并启动 RNA 合成，聚合酶沿着模板 DNA 前进并产生 RNA 拷贝(Hager, McNally, and Misteli 2009)。mRNA 产生后通常由剪接体进行剪接，被称为「内含子」的部分被移除，被称为「外显子」的部分则被留下。随后剩余的 mRNA 被翻译成蛋白质，最终成熟转录本包括哪些外显子也是可以被调控的，这使得蛋白质具有结构和功能的多样性（见图 1.3）。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200711214654.png) 图 1.3：转录后可以通过剪切产生不同的转录本，进而产生不同的蛋白质亚型，因为产生蛋白质所需的信息在转录本中编码。同一基因的不同转录本可以产生不同的蛋白质亚型。 与蛋白质编码基因相反，非编码 RNA(ncRNAs)基因在转录后经过加工即发挥功能，不进入翻译过程，因此被称为非编码 RNA，某些 ncRNA 也可以进行剪切但仍不会翻译。 ncRNA 和其他 RNA 一般可以分子内形成互补碱基，使它们具有额外的复杂性。这种基于自身互补的结构称为 RNA 二级结构，通常是许多 ncRNA 发挥功能所必需的。 综上所述，从转录起始到产生功能产物的一系列过程称为基因表达，而基因表达的量化和调控则是基因组生物学的基础研究内容。 ## 1.1.4 基因是什么样的 在进入下一话题前，我们最好先了解一下基因是如何被可视化的。作为一个对计算基因组学感兴趣的人，你会经常在电脑屏幕上看到一个基因，而它在电脑上的呈现方式就等同于你听到「基因」这个词时头脑中想到样子。在线数据库中基因会以字母序列的形式出现，或者用一系列彼此链接的方框来展示外显子和内含子的结构，其中也可能包括了转录的方向（见图 1.4）。当然，你遇到更多的是后者，所以当你想到基因时你的脑海中很可能会出现后者的样子。 我们提到 DNA 有两条链，基因其实可以位于其中任何一条链上，转录的方向也取决于此。在下图你可以看到内含子上的箭头（连接方框的线）表示基因的方向。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200711214707.png) 图 1.4: A) UCSC 浏览器中基因的表示方式。方框表示外显子，线表示内含子。B) NCBI GenBank 数据库中显示的 FATE1 基因部分序列。 ## 参考文献 Alberts, B., D. Bray, J. Lewis, M. Raff, K. Roberts, and J.D. Watson. 2002. *Molecular Biology of the Cell*. 4th ed. Garland. Eilbeck, Karen, Suzanna E Lewis, Christopher J Mungall, Mark Yandell, Lincoln Stein, Richard Durbin, and Michael Ashburner. 2005. “The Sequence Ontology: A Tool for the Unification of Genome Annotations.” *Genome Biology* 6 (5): R44. Hager, Gordon L, James G McNally, and Tom Misteli. 2009. “Transcription Dynamics.” *Molecular Cell* 35 (6): 741–53. # 第二章：基于基因组数据的 R 介绍 更新中 计算基因组学的目的是从更高维度的基因组学数据中提供生物学解释和见解。总体而言，它和任何其他类型的数据分析都类似，但是做计算基因组学需要该领域特定的知识和工具。 随着高通量实验技术的兴起，数据分析能力也成为研究者们追求的一项技能。本章的目的是首先让读者熟悉数据分析步骤，然后在基因组数据分析的背景下提供 R 编程的基础知识。R 是一种开源免费的统计编程语言，在研究人员和数据挖掘人员中很受欢迎，可以用于构建软件和进行数据分析。虽然有很多 R 编程教程可以学习，但我们的目标是在基因组学的背景中进行介绍。当你尝试用 R 分析基因组数据时，书中提到的这些例子都来自于现实工作。为了分析基因组数据而学习这种编程语言时需要根据基因组学的实际背景来对学习材料进行筛选。 # 2.1 基因组数据分析步骤 无论分析何种类型数据，数据分析都有一个共同的模式。我们将讨论这种一般模式以及如何将其应用于基因组学问题。数据分析步骤通常包括数据收集、质量检查和清理、数据处理、数据建模、数据可视化和报告几个部分。 尽管人们期望以线性方式完成这些步骤，但用不同的参数或工具重复这些步骤都是正常的。实际上，数据分析需要一遍又一遍地经历同样的步骤，以便能够: a) 回答（开始没有意识到的）其它相关问题，b) 处理后期分析中意识到的数据质量问题，以及，c) 处理在分析中加入的新数据集。 接下来，我们将在基因组数据分析的背景下简要解释这些步骤。 ## 2.1.1 数据收集 数据收集是指为需要分析的问题提供任何来源、实验或调查得到的相关数据。在基因组学中，数据收集是由第一章介绍的高通量分析完成的。我们也可以使用公开可用的数据集和在第一章中提到的那些专业数据库。你应该收集多少数据和什么类型的数据取决于你试图回答的问题以及你正在研究的技术和生物多样性。 ## 2.1.2 质量检查和清理 一般来说，数据分析几乎总是从处理不完善的数据开始。有噪音的缺失值或测量值是很常见的，数据质量检查和清理的目的在于识别数据中存在的问题并将其从数据集中清理出去。 高通量基因组学数据在产生的时候通常有技术原因引入的偏差，以测序为例，测序得到的 reads 具有不同的碱基质量。在 reads 的末尾可能就有很多质量不好的碱基。识别那些低质量碱基并删除它们可以提高比对步骤的效果。 ## 2.1.3 数据处理 这一步指的是将数据处理成适合探索性分析和建模的格式。通常我们拿到手的数据不会是可以直接进行后续分析的格式。你可能需要通过转换 (如 log 转换、标准化等) 将其调整为其他格式，或者用一些预定义条件从原始数据集中提取子集。就基因组学数据而言这些处理包括多个步骤。进行完质量检查之后，处理过程包括将 reada 与基因组比对，并对感兴趣的基因或区域进行定量。这只是计算有多少 reads 覆盖到了你感兴趣的区域，如果你的实验方案是 RNA 测序，这个数量通过后续一些标准化的方法可以让你知道每个基因表达量是多少。 ## 2.1.4 探索性数据分析和建模 这个阶段通常采用已处理或半处理过的数据并应用机器学习或统计方法对数据进行探索性分析。比较典型的内容例如我们需要看到变量之间的关系或者基于变量看到样本之间的关系。具体而言，我们可能会观察处理组是否如实验设计预期，是否有异常值或任何其他的异常。在此步骤之后，你可能需要进行额外的清理或重新处理以解决掉各种在在这一步发现的异常。 另一个相关步骤是建模，通常指的是基于你测量的其他变量来对你感兴趣的变量进行建模。在基因组学的背景下，可能你试图通过从组织样本中测量的基因表达来预测患者的疾病状态，具体方法可以是回归或任何其他机器学习方法，这个过程通常被我们称为预测建模。 统计建模也是这个步骤的一部分，如假设检验。首先我们有一个预期，然后尝试确认预期也与模型有关。一个很好的例子就是差异基因表达分析，比较某种条件下的两个数据集，如条件 A 和条件 B 的表达值，我们假设条件 A 和条件 B 具有相似的表达值然后进行检验。你将在第三章中看到更多相关信息。 ## 2.1.5 可视化和报告 可视化是所有步骤或多或少都需要的内容，但是在最后阶段你的分析结果最终需要以数字、表格和文本的方式进行呈现，这就是你的分析报告。在基因组学中，我们会使用常见的数据可视化方法以及由基因组数据分析开发或推广的一些特定可视化方法。你会在第三章看到很多流行的可视化内容。 ## 2.1.6 为什么使用 R 进行基因组学？ 凭借天生的统计分析能力、绘图优势和丰富的扩展包，R 是分析基因组数据的最佳语言之一。高维基因组数据集通常适合用核心 R 包和函数进行分析，最重要的是 bioconductor 和 CRAN 有一系列专门的工具来进行基因组学特异性分析。以下是可以使用 R 完成的计算基因组学任务列表。 ### 2.1.6.1 数据清理和处理 大多数数据清理任务，例如删除不完整的列和值、重组和转换数据都可以使用 R 实现。此外，在 R 包的帮助下还可以连接到各种格式的数据库，如 mySQL，mongoDB 等，并使用数据库特定工具查询和获取数据到 R 环境中。 除此之外，基因组数据的特定处理和质量检查可以通过 R/bioconductor 实现。例如，reads 比对质量检查甚至比对也可以通过 R 包实现。 ### 2.1.6.2 一般数据分析和探索 大多数基因组数据集也适用于一般数据分析工具的应用。在某些情况下，你可能需要对数据进行预处理以使其处于适合应用某类工具的状态。以下是通过 R 可以做的部分事情。 * 无监督数据分析: 聚类 (k-means，分层)，矩阵分解 (PCA，ICA 等) * 监督数据分析: 广义线性模型、支持向量机、随机森林 ### 2.1.6.3 特定基因组学数据分析方法 R/Bioconducto 让你可以实现许多其他生物信息学特定的算法。以下是可以做的部分事情。 * 序列分析: 给定 DNA 序列的 TF 结合 motif、 GC 含量和 CpG 计数等 * 差异表达 (芯片和测序数据) * 基因集/通路分析: 基因集中富集了什么样的基因 * 基因组区间操作，如与转录起始位点重叠的 CpG 岛，以及基于位置重叠的过滤 * 与外显子重叠的 reads 数和计算每个基因的 reads 数 ### 2.1.6.4 可视化 可视化是包括计算基因组学在内的所有数据分析技术的重要组成部分。同样，你可以在 R 中使用基本可视化技术，也可以在特定包的帮助下使用基因组相关的特定技术。这里是部分可以用 R 做的事情。 * 基本图: 直方图，散点图，柱状图，箱线图，热图 * 基于全基因组的 ideograms 和 circos 图提供了整个基因组不同特征的可视化。 * 基因组特征区间的 meta-profiles ，如所有启动子上的 read 富集展示 * 基因组中特定位点定量分析的可视化展示 # 2.2 从 R 起步 下载并安装 [R](http://cran.r-project.org/) 和 [RStudio](http://www.rstudio.com/) 。虽然 Rstudio 并非必须但它是一个很好的工具，如果你刚刚开始学习 R 非常建议你使用 Rstudio 进行学习，你需要特定的数据集来运行本文档中的代码。下载 data.zip 并将其解压到你选择的目录，文件夹名称应该是‘data’，你的 R 工作目录应该在 data 文件夹上一层级，也就是在你的 R 控制台中输入 `dir("data")` 时，应该能够看到数据文件夹的内容。 你可以通过 `setwd()` 命令更改工作目录，使用 `getwd()` 命令获取当前工作目录。在 RStudio 中也可以单击顶部菜单并通过可视化的操作来更改工作目录的位置。 > 译者注： 目前实际并没有 data.zip 这个数据，本书中使用到的数据可以通过 devtools::install\_github("compgenomr/compGenomRData") 进行安装。 > > 同时我也把这个数据包上传到了网盘，你可以下载后在本地通过 install.packages("\~/Desktop/compGenomRData.tar.gz", repos=NULL,type="source") 进行安装 > > 网盘链接: 密码: `24mh` > > 备用链接： 密码：`kaopubaer` ## 2.2.1 安装包 R 包可以理解为基础 R 的附加内容，可以帮你实现基础 R 中不直接支持的任务。正是通过这些扩展包才让 R 成为适合计算基因组学的工具。[Bioconductor](http://bioconductor.org/) 项目是计算生物学相关软件包的专用库，同时 R 的主包存储库 CRAN 也有计算生物学相关的包。除此以外，[R-Forge](http://bioconductor.org/)，[GitHub](http://bioconductor.org/) 和 [googlecode](http://bioconductor.org/) 也可能托管了部分 R 包。 你可以用 `install.packages()` 安装 CRAN 包(需要说明， `#` 是 R 中的注释字符)。 ```r # 从 CRAN 安装名为 "randomForests" 的 R 包 install.packages("randomForests") ``` 你可以使用特定的安装方法来安装 bioconductor 包。 ```r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("rtracklayer") ``` 可以使用 devtools 的 `install_github()` 函数从 github 安装软件包。 ```r library(devtools) install_github("hadley/stringr") ``` > 译者注：从 GitHub 或者 Bitbucket 等地方安装 R 包，还可以借助 [remote 包](https://remotes.r-lib.org/#:~:text=remotes%20Install%20R%20Packages%20from%20remote%20or%20local,lightweight%20replacement%20of%20the%20install_%2A%20functions%20in%20devtools.)进行安装。Remote 本身可以理解为 devtools 一系列 install\_\* 函数的精简版本。remote 同时支持下载 bioconductor 包。 安装软件包的另一种方法是从源代码进行本地安装。 ```r # 下载源文件 download.file("http://goo.gl/3pvHYI", destfile="methylKit_0.5.7.tar.gz") # 从源文件安装包 install.packages("methylKit_0.5.7.tar.gz", repos=NULL,type="source") # 删除源文件 unlink("methylKit_0.5.7.tar.gz") ``` 你还可以更新 CRAN 和 Bioconductor 包。 ```r # 升级 CRAN 包 update.packages() # 升级 bioconductor 包 BiocManager::install(update = T) ``` ## 2.2.2 在自定义位置安装包 如果你在服务器或集群上使用 R，那就不太可能拥有管理员权限来安装包。在这种情况下可以通过告诉 R 在哪里寻找额外的包来自定义位置。 打开你家目录中的 renvironon 文件，并添加以下行: ``` R_LIBS=~/Rlibs ``` 这句命令将告诉 R 家目录的 'Rlibs' 目录是查找包和安装包的第一位置选择，接下来你应该去创建那个目录然后启动一个新的 R 会话并开始安装包。这之后的 R 包将被安装到你具有读写访问权限的本地目录中。 > 译者注：关于 R 的相关配置和在 macOS 中的安装，可以参考译者早前写过的两篇文章 [R 安装升级后的若干规定动作](https://kaopubear.top/blog/2018-07-09-chineseuser/) 和 [macOS 10.15 安装 R 包](https://kaopubear.top/blog/2019-10-29-macos15user/) 。 ## 2.2.3 获取有关函数和包的帮助信息 你可以通过 `help()` 和 `help.search()` 函数获得有关函数的帮助文档。也可以用 `ls()` 函数列出一个包中的函数 ```r library(MASS) ls("package:MASS") # functions in the package ls() # objects in your R enviroment # get help on hist() function ?hist help("hist") # search the word "hist" in help pages help.search("hist") ??hist ``` ### 2.2.3.1 需要更多帮助？ 此外，你还可以检查包的 vignette 以获得帮助和对函数更深入实际的理解。所有 Bionconductor 包都有 vignette 来引导你完成示例分析。学会搜索也总是有帮助的，有许多博客和网页上都有关于 R 的帖子，Stackoverflow 和 R-blogger 通常是优秀和可靠的信息来源。 # 2.3 R 中的计算 R 可以当作普通计算器来使用，也有人说 R 其实就是一个非常复杂的计算器。这里有一些例子，请记住 `#` 是注释字符，注释给出了操作的细节含义。 ```r 2 + 3 * 5 # 注意运算符的优先级. log(10) # 以e为底的自然对数 5^2 # 5的平方 3/2 # 除法 sqrt(16) # 开平方 abs(3-7) # 绝对值 pi # 数字 exp(2) # 指数函数 # 这是注释 ``` # 2.4 R 的数据结构 R 具有多种数据结构，了解 R 中常见的数据结构以及如何使用它们是至关重要的。 ### 2.4.1 向量 向量是一个核心的数据结构。它可以理解为一个相同类型 (数字，字符或逻辑) 的元素列表。稍后你将看到表的每一列都将表示为向量。R 可以轻松直观地处理向量，可以用 `c()` 函数创建向量，但这不是唯一的方法。对向量的操作将作用于向量的所有元素。 ```r x<-c(1,3,2,10,5) #创建含有5个元素的向量 x = c(1,3,2,10,5) x ``` ``` ## [1] 1 3 2 10 5 ``` ```r y<-1:5 #创建一个包含5个连续整数的向量 y+2 #加法运算 ``` ``` ## [1] 3 4 5 6 7 ``` ```r 2*y #乘法运算 ``` ``` ## [1] 2 4 6 8 10 ``` ```r y^2 #对每个元素平方运算 ``` ``` ## [1] 1 4 9 16 25 ``` ```r 2^y #对数字2进行对应元素的次方运算 ``` ``` ## [1] 2 4 8 16 32 ``` ```r y #y 本身并不会被改变 ``` ``` ## [1] 1 2 3 4 5 ``` ```r y<-y*2 y #此时y发生了变化 ``` ``` ## [1] 2 4 6 8 10 ``` ```r r1<-rep(1,3) # 创造一个长度为3的向量 length(r1) #向量长度 ``` ``` ## [1] 3 ``` ```r class(r1) # 向量类型 ``` ``` ## [1] "numeric" ``` ```r a<-1 # 实际这是一个长度为1的向量 ``` ### 2.4.2 矩阵 矩阵是指由行和列组成的数字数组。你可以将其视为向量的一个叠加版本，其中每行或每列都是向量。创建矩阵的最简单方法之一是使用 `cbind()` 组合相等长度的向量，其含义是 “通过列合并”。 ```r x<-c(1,2,3,4) y<-c(4,5,6,7) m1<-cbind(x,y);m1 ``` ``` ## x y ## [1,] 1 4 ## [2,] 2 5 ## [3,] 3 6 ## [4,] 4 7 ``` ```r t(m1) #转置 ``` ``` ## [,1] [,2] [,3] [,4] ## x 1 2 3 4 ## y 4 5 6 7 ``` ```r dim(m1) # 返回对象维度 ``` ``` ## [1] 4 2 ``` 你也可以直接列出元素并指定矩阵: ```r m2<-matrix(c(1,3,2,5,-1,2,2,3,9),nrow=3) m2 ``` ``` ## [,1] [,2] [,3] ## [1,] 1 5 2 ## [2,] 3 -1 3 ## [3,] 2 2 9 ``` 矩阵和另一个数据结构数据框都是是表格型的数据结构，你可以按需提取行和列提供给子集。图 2.1 显示了它们是如何工作的。 ![](https://kaopubear-1254299507.cos.ap-shanghai.myqcloud.com/picgo/20200718190453.png) 图 2.1: 从矩阵中提取子集 ### 2.4.3 数据框（Data Frames） 数据框比框矩阵更加通用，因为不同列可以是不同的数据类型 (数字，字符，因子等)。可以通过 `data.frame()` 函数构造数据框。接下来我们说明了如何从基因组区段或者坐标构建数据框。 ```r chr <- c("chr1", "chr1", "chr2", "chr2") strand <- c("-","-","+","+") start<- c(200,4000,100,400) end<-c(250,410,200,450) mydata <- data.frame(chr,start,end,strand) #修改列名 names(mydata) <- c("chr","start","end","strand") mydata ``` ``` ## chr start end strand ## 1 chr1 200 250 - ## 2 chr1 4000 410 - ## 3 chr2 100 200 + ## 4 chr2 400 450 + ``` ```r # 另一种方法 mydata <- data.frame(chr=chr,start=start,end=end,strand=strand) mydata ``` ``` ## chr start end strand ## 1 chr1 200 250 - ## 2 chr1 4000 410 - ## 3 chr2 100 200 + ## 4 chr2 400 450 + ``` 有多种方法可以提取数据框的元素，你可以使用列数或列名来提取某些列，也可以使用行号提取某些行，还可以使用逻辑参数来提取数据，例如提取列中值大于某个阈值的所有行。 ```r mydata[,2:4] # 提取2-4列 ``` ``` ## start end strand ## 1 200 250 - ## 2 4000 410 - ## 3 100 200 + ## 4 400 450 + ``` ```r mydata[,c("chr","start")] # 提取chr和start两列 ``` ``` ## chr start ## 1 chr1 200 ## 2 chr1 4000 ## 3 chr2 100 ## 4 chr2 400 ``` ```r mydata$start # 数据框中的start变量 ``` ``` ## [1] 200 4000 100 400 ``` ```r mydata[c(1,3),] # 提取第一和第三行 ``` ``` ## chr start end strand ## 1 chr1 200 250 - ## 3 chr2 100 200 + ``` ```r mydata[mydata$start>400,] # 提取所有start大于400的行 ``` ``` ## chr start end strand ## 2 chr1 4000 410 - ``` ### 2.4.4 列表 列表可以理解为对象 (组件) 的有序集合。列表允许你收集各种 (可能不相关的) 对象。 ```r # 具有4个成分的列表事例 # 字符串，数值向量，矩阵和标量 w <- list(name="Fred", mynumbers=c(1,2,3), mymatrix=matrix(1:4,ncol=2), age=5.3) w ``` ``` ## $name ## [1] "Fred" ## ## $mynumbers ## [1] 1 2 3 ## ## $mymatrix ## [,1] [,2] ## [1,] 1 3 ## [2,] 2 4 ## ## $age ## [1] 5.3 ``` 您可以使用 `[]` 用列表中的位置或名称提取列表的元素。 ```r w[[3]] # 列表的第三个元素 ``` ``` ## [,1] [,2] ## [1,] 1 3 ## [2,] 2 4 ``` ```r w[["mynumbers"]] # 名字为 mynumbers 的元素 ``` ``` ## [1] 1 2 3 ``` ```r w$age ``` ``` ## [1] 5.3 ``` ### 2.4.5 因子 因子用于存储分类数据，它们对于统计建模很重要，因为分类变量在统计模型中与连续变量会被区别对待。这确保了在统计模型中可以正确地处理分类数据。 ```r features=c("promoter","exon","intron") f.feat=factor(features) ``` 需要注意的一点是，当你使用`read.table()` 来读取数据框或者使用 `data.frame()` 来创造数据框时，字符列默认被存储为因子，如果想要修改这个默认设置可以在这两个函数中设定 `stringsasfactor = FALSE` 。 > 译者注：R 4.0 版本开始，默认设置已经为 `stringsasfactor = FALSE`