bioconductor（Bioconductor模块是什么）

本文目录

• Bioconductor模块是什么
• Lecture 7 Affymetrix，R与bioconductor（芯片数据预处理及质量控制）
• bioconductor入门第二弹--GRanges概述

Bioconductor模块是什么

Bioconductor就是一个基于R语言的、面向基因组信息分析的应用软件集合。Bioconductor的应用功能是以包的集成形式呈现在用户面前，它提供的软件包中包括各种基因组数据分析和注释工具，其中大多数工具是针对DNA微阵列或基因芯片数据的处理、分析、注释及可视化的。同时，Bioconductor。还提供许多与DNA微阵列相关的数据包。

Lecture 7 Affymetrix，R与bioconductor（芯片数据预处理及质量控制）

***隐藏网址***

芯片数据的结构如下

对于 Sample Information ，我们至少需要知道样本的组别，在一些更加复杂的使用设计中，我们可能还需要知道更多关于样本的信息如：

在R的数据框类型数据中，每一行代表一个实验样本，每一列代表每一个我们感兴趣的因素（包括处理因素），dataframe可以同时处理文本因子类型和数值类型的值因而非常实用

通常我们在excel中创建一个储存样本信息的表格【每一行代表一个样本，每一列代表一个因素】，并存储为csv或tab分隔的文本文件，再用R读取为dataframe【read.csv，read.delim，read.table函数】

Bioconductor的 Biobase 包通过一个叫做 AnnotatedDataFrame 对象来存储样本信息，这个对象包括：

如我们有8个样本的芯片数据，每个芯片都有200个基因探针：

接下来我们需要创建 AnnotatedDataFrame 对象来描述这8个样本的表型：

可以使用pData查看 AnnotatedDataFrame 对象的样本表型信息

在bioconductor中， ExpressionSet 对象可以把表达矩阵和样本信息整合在一起

可以直接查看 ExpressionSet 的一些基本信息(基因数，样本数等)

我们创建好了一个 ExpressionSet 对象，但是这个对象却没有对基因的注释及相关信息，bioconductor的一个瑕疵就是它允许用户把表达数据与基因信息 完全分离开 （但大多数人对此都欣然接受）

处理表达数据（exprs）和样本信息（phenoData）外 ExpressoinSet 对象还允许添加其他的一些可选信息

对于Affymetrix公司的芯片数据，bioconductor提供了一个 特殊的ExpressionSet数据类型 ，即 AffyBatch ， affydata 包提供了一个示例数据 Dilution ， Dilution 即为一个AffyBatch对象

首选我们可以看一下芯片的扫描图像，如果出现了异常亮区或异常暗区，则说明质量可能存在问题（被刮了一下）如下图

也可以用箱式图及密度分布图查看每个芯片的表达量分布（注意表达量自动进行log2处理），如果出现了某个芯片与其他芯片的差异很大，说明该芯片质量可能存在问题

affy 包也提供提取特定探针数据的函数

可看出一个序列（一个基因或RNA对应一个probset）对应16个探针（通常一个probset对应11~20种探针，每种探针会重复数百次）

对该基因的探针数据（未经过log2处理）进行可视化：

同时affy包也能查看RNA的降解情况，因为RNA的降解都是从5`向3`开始的，这16（通常每个probset即每个RNA对应这11~20个不同位置的探针）个探针的靶序列也是从RNA的5`向3`方向逐渐靠近的

以杂交的信号为纵坐标，探针的位置作为横坐标作图（其中纵坐标进行了缩放和移动以便对不同样本的斜率进行比较），可以推测5`端的探针序列对应的RNA相对与3`来说密度较低，信号也越低，因此曲线整体呈上升趋势

比斜率更重要的是芯片之间的比较，理想状况下各芯片的曲线是平行的，而降解得很多的样品其5’ 到 3’ 线段的斜率会比其他样品大，该数据显示每个位置的探针降解速度都不太一致（斜率呈分段变化），但样品之间的降解一致 ，通常从5‘到3’端的第5个探针认为是比较稳定的

下图为存在异常降解芯片 脱颖而出 的数据：

affy包提供了 ReadAffy 函数使创建AffyBatch对像变得简单，该函数还提供了GUI（graphical user interface）界面方便用户的使用，然而我们需要对分析过程进行文档化，为了分析的可重复性我们应该避免使用GUI界面

由CEL文件创建AffyBatch对象

构建样本信息的数据AnnotatedDataFrame对象

其中样本的信息文件（data_file_list.txt）格式如下：

使用函数Biobase包的read.AnnotatedDataFrame读取注释文件，设置第二列为行名：

也可以再构建一个description，即MIAME对象

另一种方法为之间使用read.affybatch函数之间构建AffyBatch对象

这四个步骤为：

MAS 5.0 的背景校正方法为把每个芯片（每个CEL文件）分为16个区，在每个区中最最暗的2%作为每个区的背景水平，在对该芯片的每个探针用这16个背景值的加权平均进行校正，每个区的权重取决于该探针与芯片中这16个区域中心的距离。

RMA算法只对PM进行背景校正，MM的值不做校正（放弃MM，因为几乎有1/3的MM值是大于对应的PM的），该模型的假设为PM的值是由两个部分组成的即S=X+Y，其中S为实际检查到的信号，X代表有探针杂交产生的信号服从指数分布Exp(α)，Y代表背景信号服从正态分布N(μ, δ^2)，Y只取大于0的部分，并且X与Y是独立的

***隐藏网址***

还有一个gcrma的算法，考虑到不同序列gc含量与序列亲和力的关系进而估计不同背景信号大小，需要安装gcrma包 bgc.gcrma = gcrma(Data) ，该命令包括了背景校正，标准化和summarization这几个步骤。

在该示例数据中，RMA算法估计的背景值比MAS 5.0估计的背景值大，而且RMA算法的背景值比较恒定

差值大部分为20，可以由芯片的信号值分布看一看20究竟有多大

可以看出20越占大部分信号值的1/10，不同算法的差值还是比较大的

暂时先不管normalization和PM correction这一步，summarization的目的就是把所有同一probset（即同一基因mRNA）的PM和MM探针对转化为一个数值，这个数值是对目标基因表达量的最佳估计，其方法有：

使用不同的summarization方法对probeset进行定量， expresso 函数把所有预处理步骤都整合在了一起

MA plot的作用是为了展示两个值几乎处处相等的变量（x和y）之间的关系，为了展示两个变量之间的变化关系，大多数人的思维都是把x与y分别作为横轴和纵轴进行绘图，如果y=x，则该图呈45度角的直线（如下图中左边图的蓝色直线），可以通过查看点形成的直线偏离预期直线的多少来衡量系统偏差，然而该图存在以下几个缺点：

MA plot的处理方法是把该直线顺时针旋转45度，把参考对角线变为直线，具体做法是把(x+y)/2作为横轴，(y-x)作为纵轴，则参考直线变为一条水平线，如下方右图，这样可以很清楚的在视觉上展示两个相等的变量之间偏离参考值的大小，即存在的系统误差的大小

affy包的函数 mva.pairs 可以快速地做出MA plot（也可以使用MAplot函数），接下来我们使用MA plot来比较不同summarize方法之间的差异， mva.pairs 函数的输入参数就是各变量形成的一个dataframe

上图中的差异也有可能是由于我们没有在summarization之前进行校正，接下来我们使用RMA方法进行背景校正，使用百分位数标准化方法，并且只用PM值进行数据定量汇总

不同算法的组合如下表（摘抄自：NAR文献 Comparison of algorithms for the analysis of Affymetrix microarray data as evaluated by co-expression of genes in known operons ）：

MAS4.0代表一个完整的数据处理算法，对背景校正，标准化，PM校正，probeSet定量4个完整步骤都有其特定的计算过程：

对 avgdiff 的评价：

mas5.0的4个步骤中，除了标准化为constant外，其余均名为mas

PM校正方法名为mas：即把PM-CT作为校正的PM，CT为change threshold是一个调整后的MM，以保证CT值小于PM值，CT值算法为

summarization方法也称为mas：即Signal = anti-log of Tukeybiweight(log (PM –CT))

点评：mas 仍然没有用到不同芯片的数据即建立一个模型

开发了mas 5.0后，affymetrix公司不再投入大量精力开发算法，但提供了一套可以检验算法的数据The Affy Latin Square Experiment.用这个数据可以证明mas 5.0比mas 4.0（avgdiff）表现更优越

Li-Wong算法【由现北大研究员 李程 与现斯坦福教授 王永雄 开发】包括三个小部分，不包括背景校正过程，标准化为invariantset（以后再了解），PM校正方法为只使用PM值或PP和MM值都使用（可自定义），summarization为一个统计学模型称为Li-Wong算法，所有芯片的数据都得到了很好的利用，但是如果芯片数目比较少时可信度较差

***隐藏网址***

RMA 算法是2003年由Rafael A. Irizarry（Bekerley）等人提出的，前面已经大致介绍了背景校正的方法，标准化使用的是一种叫做median polish的方法（以后再介绍），该方法完全丢弃了MM值（他们给出的理由是因为太多MM信号值大于PM值的情况，因此考虑MM的值会引入很大的变异）【他们的想法或许是正确的】

像dChip一样，该算法的summarization方法是建立在一个统计学模型之上的：

其中芯片编号为i，j代表探针

模型的参数由多个芯片的数据进行拟合得到，与dChip不同的是$\epsilon$的值是经过对数处理的，这充分考虑到了值越大的探针其变异也越大的特点

Li Zhang提出了Position-Dependent Nearest Neighbor(PDNN) 模型( Nat Biotech, 2003; 21:818 ).不像dChip和RMA，PDNN模型参数的估计只需要单个芯片的数据（很大程度上是由于参数数目比较少）。该模型考虑序列对杂交的影响，把信号强度值与序列联系起来。

该模型的参数有：

该模型的实现需要单独安装一个包 affypdnn ，在载入这个包之后，会自动载入相关函数，并更新可选的模型

需要注意的是PDNN模型没有按照标准的 expresso 4个步骤处理数据，该模型直接从百分位数校正开始，可以只使用PM或PM和MM同时进行拟合，背景的估计整合到了能量和表达参数里面。

需要使用 expressopdnn 这个 expresso 的变体实现PDNN模型的运算。

第一步是看图形，有没有异常的区域

接着可以看看表达值的分布情况，看有没有异常的芯片

首先我们可以看看每个芯片的背景值是否有可比性，差别非常大的异常值

芯片标准化时，标准的Affymetrix流程为乘上一个校正因子（scaling factors）使各芯片的中位数相同，Affymetrix要求可比较的芯片之间的校正因子不能超过3

接下来是percent of present值（即有多少基因按照 Detection calls 算法是确实存在表达且被检测到的）是否有可比性，如果有离群值或低于30，高于60的值说明可能存在质量问题

Affymetrix设计了一些参考管家基因（通常是GAPDH或beta-actin）的3‘，middle和5’的探针，其中actin的3‘/5’的值应该是小于3，GAPDH的3’/5’值小于1（5’端可能存在部分降解，可与RNA降解图结合起来看），否则可能存在问题

也可以之间对qc结果进行作图查看

affyPLM 是一个与RMA模型相似的probe-level model (PLM)从而进行质控

***隐藏网址*** 当污点比较大时，我们需要再看其他的质控数据决定是否应该舍弃该芯片 。

下面这个就是有问题的图片，除原始图片不太清楚外，其他图片都明显显示有人为误差，即一个“环”，就需要看下一步的QC数据是否有较大异常决定去留

还可以用箱式图查看M值（相对于表达量为中等的那个芯片）分布也叫做Relative Log Expression的分布。 如果有芯片的中位数不是0或跨度明显大于其他芯片，则可能存在误差，这个两个图都是一模一样的，不过名称不同而已

接下来就是看Normalized Unscaled Standard Error (NUSE)值的分布是否异常，与M值（RLE值）的思想差不多，NUSE是每个基因的标准差与该基因中所有芯片中的标准差中位数的比值，因此芯片具有可比性的条件是其所有基因的标准差几乎相同，即NUSE的中位数为1，且分布大致相当。

bioconductor入门第二弹--GRanges概述

Granges是genomic ranges的简称，它是bioconductor中的R包存储genomic intervals的数据结构，快速且有效，是想要深入分析数据所必须掌握的。它的重要性类似于bed注释文件，这么一比较这玩意的重要性就显而易见了

做基因组学分析时，最重要的就是intervals或者是intervals的集合。那么什么是intervals呢？它是基因组上的一些实体，我也想不出来一个合适的名词去翻译它，你可以粗鲁的称呼它为对象。

虽然翻译不出来，但是我举几个例子，一定能让你明白它究竟是个什么玩意

如下图，这是UCSC browser上随便挑选出的一段序列的的截图，我们能从上面看到DNA clusters，SNPs ，promoters， genes，这些东西都可以称之为intervals，甚至一个测序read被map到基因组上的位置也可以称之为一个intervals

许多基因组工作涉及到各种intervals的各种相互关系，如哪些promoters含有SNP位点，某个TF结合位点又与哪些promoter重合，哪些基因被测序reads覆盖。说了这么多intervals，那么它和我们要介绍的GRanges有什么关系呢？

上图是一个GRanges的结构，其中的A1，A2和A3是它的三个intervals，所以我们可以理解，intervals是GRanges对象的横向最小单位。

从纵向上看，GRanges对象由三列组成，第一列是染色体名称，第二列是intervals所在的具体位置（这个位置信息是IRanges对象所提供的），第三列则是该位置所在链的方向

除此之外，我们还能看到左下角是染色体的信息，比如seqlengths告诉了我们chr1的长度，若未定义则为NA

我们可以使用R来处理GRanges对象，但除此之外我们还可以使用bedtools处理intervals之间的相互关系。比如，我最常用的是intersectBed，它可以很快的帮助我知道我所需要的信号处于哪些位置。

IRanges是包含intervals的向量， IRanges和GRanges类似，但是GRanges具有染色体编号和正负链信息，IRanges为GRanges提供了基础功能，可以进行intervals的简单代数运算

IRanges有三列，起始，终止以及宽度。当我们给定了任意两列信息就可以由此推断剩下一列的信息。

下面使用IRanges函数构建IRanges对象

这个个Iranges对象包含三个不同的intervals，每一个intervals都称之为一个range(范围区间)，这里的ir1有三个ranges

start(), end(), width()这三个方法可以获得IRanges对象的起始列，终止列，以及宽度那一列的信息。IRanges对象的每个ranges也可以有自己的名字，如

因为是向量，因此IRanges只有一个维度，不能使用dim()方法，但是我们可以使用length()获取其长度。我们可以像操作向量一样，对IRanges对象进行操作，如

甚至可以利用c()将两个IRanges对象合并

有时候在阅读R包的资料时，你也许会看到NormalIRanges对象，它是由使用reduce函数对IRanges对象处理得到。一个NormalIRanges，他其中的ranges（就是每一 行 数据）有以下几个特点：

用文字可能比较抽象，用一张图应该能帮助大家更好的理解。上面的表示我们的其实Iranges对象，而下面的则是我们使用reduce函数处理后所得到的normal Ranges对象。

下面是reduce()函数的对象：

disjoin函数和reduce函数恰好相反，找出完全不重合的intervals，将重合的intervals进行分割。

还有一些方法可以用于批量转换IRanges对象的坐标位置，比如：

除此之外，还有一些方法对于上述三个方法做了进一步的拓展，由于篇幅原因就不在这里介绍，大家感兴趣可以自己搜搜看

reseize()函数可能是一个处理数据很有用的方法，它有一个fix参数，可以固定原始ranges中不变的列，具体用法如下:

既然IRanges是一个向量，那么肯定可以对IRanges进行集合操作，我们通常使用unnion取并集，intersect取交集：

这些函数，还有一些变种，punion(), pintersect(), psetdiff(), pgap()。它们与R基础函数中的pmax很相似，很少被用到，因此不做讲解

对测序数据进行下游分析时，经常需要寻找overlap，使用findOverlaps()函数寻找两个IRanges对象中的overlap。它的返回值是一个Hists对象，它可以描述两个Iranges之间重叠关系，是一个两列的矩阵，由IRanges的索引组成。可以由queryHist()和subjectHits()得到其索引信息。

findOverlaps的参数很多，是否有最小重叠，是否完全重叠等。我自己也记不住，所以就不在这里介绍了

countOverlaps()函数可以返回overlap的数目，该功能速度更快，占用内存更小，因为不需要记录ranges重叠的具体信息，仅仅考虑重叠的个数

nearest()函数用于寻找两个IRanges中，第一个IRanges对象中的ranges最接近第二个IRanges对象中的哪个ranges。

如果你用过Y叔的R包ChIPseeker，你可能会对这个功能比较熟悉，它将peak注释到基因的一个方法就是通过找该peak与哪个基因最接近

类似的函数还有precede，follow等，实在是太多了，而且不常用，就不一一介绍了……

相信看到这里你们对于IRanges对象已经有了一个比较深的认识，对于GRanges也有了大致的了解。下一期不出意外，我会详细写写GRanges对象的操作，希望能带你进一步认识GRanges对象