bioconductor 芯片表达聚类是哪个包
发布网友
发布时间：2022-04-22 10:23
共1个回答
热心网友
时间：2023-10-10 15:08
原笔记地址： http://x2yline.gitee.io/microarray_data_analysis_note/
芯片数据的结构如下
对于Sample Information，我们至少需要知道样本的组别，在一些更加复杂的使用设计中，我们可能还需要知道更多关于样本的信息如：
在R的数据框类型数据中，每一行代表一个实验样本，每一列代表每一个我们感兴趣的因素（包括处理因素），dataframe可以同时处理文本因子类型和数值类型的值因而非常实用
通常我们在excel中创建一个储存样本信息的表格【每一行代表一个样本，每一列代表一个因素】，并存储为csv或tab分隔的文本文件，再用R读取为dataframe【read.csv，read.delim，read.table函数】
Bioconctor的 Biobase 包通过一个叫做 AnnotatedDataFrame 对象来存储样本信息，这个对象包括：
如我们有8个样本的芯片数据，每个芯片都有200个基因探针：
接下来我们需要创建 AnnotatedDataFrame 对象来描述这8个样本的表型：
可以使用pData查看 AnnotatedDataFrame 对象的样本表型信息
在bioconctor中， ExpressionSet 对象可以把表达矩阵和样本信息整合在一起
可以直接查看 ExpressionSet 的一些基本信息(基因数，样本数等)
我们创建好了一个 ExpressionSet 对象，但是这个对象却没有对基因的注释及相关信息，bioconctor的一个瑕疵就是它允许用户把表达数据与基因信息完全分离开（但大多数人对此都欣然接受）
处理表达数据（exprs）和样本信息（phenoData）外 ExpressoinSet 对象还允许添加其他的一些可选信息
对于Affymetrix公司的芯片数据，bioconctor提供了一个特殊的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' 线段的斜率会比其他样品大，该数据显示每个位置的探针降解速度都不太一致（斜率呈分段变化），但样品之间的降解一致 [1] [2] ，通常从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是独立的
最后校正过后的X值=E(X|S = s)，其中s为该探针的信号强度，E(X|S = s)表示在检测到的信号强度为s的情况下X的期望值，具体推导可以看大牛的博士论文： http://bmbolstad.com/Dissertation/Bolstad_2004_Dissertation.pdf ，暂时没有时间来研究，先搁在这里。
还有一个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算法，所有芯片的数据都得到了很好的利用，但是如果芯片数目比较少时可信度较差
[图片上传失败...(image-b9fc36-1525698429039)]
[图片上传失败...(image-b1b7e2-1525698429039)]
其模型最终PM-MM值为探针亲和力$\varphi_{j}$与表达量$\theta_{i}$的乘积，i代表样本，j代表探针，$v_{i}$代表非特异性杂交的基线水平，文章地址为： http://www.pnas.org/content/98/1/31
RMA 算法是2003年由Rafael A. Irizarry（Bekerley）等人提出的，前面已经大致介绍了背景校正的方法，标准化使用的是一种叫做median polish的方法（以后再介绍），该方法完全丢弃了MM值（他们给出的理由是因为太多MM信号值大于PM值的情况，因此考虑MM的值会引入很大的变异）【他们的想法或许是正确的】
像dChip一样，该算法的summarization方法是建立在一个统计学模型之上的：
[图片上传失败...(image-a12c9f-1525698429039)]
其中芯片编号为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)从而进行质控
模型拟合的残差可以用image函数来可视化显示，可以很清楚地判断是否有异常残差值的出现，网址 http://www.plmimagegallery.bmbolstad.com/ 收集了很多特征性的图片，而Dilution这个数据集基本正常，通常该图片没有瑕疵的情况很少，小的瑕疵一般都可以接受，通过RMA这个稳健的算法基本可以把误差消除掉，而当污点比较大时，我们需要再看其他的质控数据决定是否应该舍弃该芯片。
下面这个就是有问题的图片，除原始图片不太清楚外，其他图片都明显显示有人为误差，即一个“环”，就需要看下一步的QC数据是否有较大异常决定去留
还可以用箱式图查看M值（相对于表达量为中等的那个芯片）分布也叫做Relative Log Expression的分布。 如果有芯片的中位数不是0或跨度明显大于其他芯片，则可能存在误差，这个两个图都是一模一样的，不过名称不同而已
接下来就是看Normalized Unscaled Standard Error (NUSE)值的分布是否异常，与M值（RLE值）的思想差不多，NUSE是每个基因的标准差与该基因中所有芯片中的标准差中位数的比值，因此芯片具有可比性的条件是其所有基因的标准差几乎相同，即NUSE的中位数为1，且分布大致相当。
