专注于Agilent microRNA 芯片数据的处理R包-AgiMicroRna

希望所有学员都可以站在生信技能树的舞台上发光发热！确实没想到如此小众的R包也可以有详细的笔记教程:

下面是EIM伟随机投稿

1. R包简介

R包作者：Pedro Lopez-Romero 最后一次更新：October 27，2020

AgiMicroRna主要用于Agilent microRNA array数据的处理、质量评估和差异表达分析。该包基于limma结构生成处理数据，以便使用limma的线性模型进行差异表达的统计推理。

AgiMicroRna将Agilent Feature Extraction (AFE)图像分析软件导出的扫描数据读入R。

AgiMicroRna主要目的是产生处理后的数据，并使用limma软件包进行分析。

该包基于miRBase数据库：

In addition, HTML files that contains links of the declared significant microRNAs to the Sanger miRBase http://microrna.sanger.ac.uk/ are given.

1.1 示例数据特征

The data come from human mesenchymal stem cells obtained from bone marrow.100 ng of each RNA sample were hybridized onto Agilent Human microRNA Microarray v2.0 (G4470B,Agilent Technologies).

The Human microRNA microarray v2.0 contains 723 human and 76 human viral microRNAs, each of them replicated 16 times. There are 362 microRNAs interrogated by 2 different oligonucleotides, 45 microRNAs by 3 and 390 microRNAs interrogated by 4 different oligonucleotides. Only 2 microRNAs are interrogated by the same oligonucleotide. The array contains also a set of positive and negative controls that are replicated a different number of times.

人骨髓间充质干细胞的microRNA，

共有三种治疗效果A,B,C，每个实验条件2个重复；

将两种处理(MSC_B和MSC_C)与对照MSC_A进行比较。

1.2 示例数据处理路线

• 准备数据&数据可视化(数据检查)

• be normalized between arrays；

• The signal is background corrected using the exponential + normal convolution model；
• the background signal is normalized between arrays；
• the total gene signal is estimated from a linear model that takes into account the probe effect. The estimates of the model are obtained using a robust methods such as the median polish.
• Obtaining the Total microRNA Gene Signal processed by AFE；
• normalization between arrays.
• 通过AFE算法计算TGS(Total Gene Signal)
• 利用RMA算法评估 the gene signal

• 质控并输出至 ExpressionSet

• After obtaining the normalized total gene signal, some of the genes are eliminated from the analysis using some of the quality flags that AFE attaches to each feature.

• 差异分析：利用 limma 的线性模型特征完成了差异表达分析。

• M value, moderated t and F statistics, p values and FDR, etc

http://microrna.sanger.ac.uk/

2. R包所需数据：Target File

需要一个目标文件，以便将每个数据文件分配给指定的实验组。【便于之后导入数据并配对】

目标文件是一个由用户创建的以tab分隔的文本格式文件。以下列必须出现在目标文件中。

• 第一列***FileName***（必须），包括图像数据文件的名称。

• 第二列***Treatment***（必须），包括治疗效果。

• 第三列 GErep（必须），它以数字代码表示治疗效果，从1到n, n是治疗效果的级别数。

目标文件中的其他列是可选的。它们可能包括其他说明实验条件的解释变量的信息，如年龄、性别和考虑到实验设计的阻塞变量(配对、阻塞设计等)。

这些变量应该包含在目标文件中，以便最终在limma模型中使用。

2.1 示例数据的Targe file

Subject列：1）表示两个重复；2）便于配对分析；

2.2 获取Targe file

示例数据的Targe file

library("AgiMicroRna")
data(targets.micro)
print(targets.micro)
## FileName Treatment GErep Subject
## mscA1 mscA1.txt A 1 1
## mscA2 mscA2.txt A 1 2
## mscB1 mscB1.txt B 2 1
## mscB2 mscB2.txt B 2 2
## mscC1 mscC1.txt C 3 1
## mscC2 mscC2.txt C 3 2


导入自己的数据：readTargets

函数readTargets的作用：

1. 帮助查看txt文件是否含有必须列(FileName、Treatment、GErep)；
2. 设置行名。

library("AgiMicroRna")
targets.micro=readTargets(infile="targets.txt",verbose=TRUE)


3. 数据读取

To read the scanned data files into R we use the readMicroRnaAFE.

【注】readMicroRnaAFE函数：

但要求上述函数中所出现的列在txt内均存在；

该函数可以创建uRNAList类的新对象，即类似limma的RGList类。

1. 列替换：TGS="gTotalGeneSignal", TPS="gTotalProbeSignal", meanS="gMeanSignal", procS="gProcessedSignal"；
2. 数据预处理和差异表达分析所必须列：gTotalGeneSignal, gMeanSignal, gIsGeneDetected, ControlType, ProbeName, and GeneName.
• 其中AFE完整算法包括：Typically it contains the Multiplicatively Detrended BackgroundSubtracted Signal or the BackgroundSubtractedSignal.
• 含义：the TotalProbeSignal times the number of probes per gene；The gTotalProbeSignal is the robust average of all the processed signals for each replicated probe multiplied by the total number of probe replicates. These signals are used by AFE algorithms to estimate the gTotalGeneSignal.
• 用途：在标准化后，用于差异表达分析；
1. gTotalGeneSignal：
2. gMeanSignal：是每个探测器的原始信号值，经AFE算法获得gProcessedSignal。
3. gIsGeneDetected：逻辑值，判断基因是否被the microRNA microarray 检测到；
4. ControlType
5. The ProbeName is an Agilent-assigned identifier for the probe synthesized on the microarray.
6. The GeneName is an identifier for the gene for which the probe provides expression information.

4. 绘图-数据可视化

AgiMicroRna可直接通过gMeanSignal绘制箱线图(boxplot)、密度图(density)、MA图(MA plots)、Relative Log Expression (RLE)和样本层次聚类(hierarchical cluster of samples)等功能，可以评估数据质量，并检查处理步骤的性能。并且该包提供了一步化出图函数：qcPlots

4.1 箱线图(boxplot)&密度图(density)

两者均是针对“log2(dd.micromeanS)=5-6的密度大，数量多。

[ps] 可以将密度图看做箱线图的另一种表示形式；

4.2 MA图

The MA plots represent the fold-change (M) in the y-axis against the average log expression (A) for two given arrays.

• M为Y轴，M = (x - y)；A为X轴，A = (x + y)/2；
• y = apply(uRNAList$meanS, 1, median)
• x = uRNAList$meanS[, i]
• 该图所用到的变量值：
• dd.micro$meanS;
• dd.microGeneName  (miRNA的ID)
• dd.microControlType （feature分类）

The signal values of the reference array are computed as the median spots taken over the whole set of arrays Every kind of feature is identified with different color. 该包自定了一个参考数组(y值)，将其他数组(x值)与该参考数组进行比较。

https://www.jianshu.com/p/cdfac0bfb733

4.2.1 miRNA的ID：

uRNAListGeneName值：

• 

  https://zhuanlan.zhihu.com/p/81384798

• 各个ID转换；https://zhuanlan.zhihu.com/p/136744701

4.2.2 该包的MA图未进行log2处理

https://www.jianshu.com/p/cdfac0bfb733

[ps]：作者已经进行添加了ddaux$G列，可能是想进行log2处理，但导入写函数时并未使用到G列。

4.3 RLE plot

RLE plot 表示经过the Relative Log Expression (RLE)算法处理的boxplot图；

The RLE plot displays for each sample a Boxplot with the Relative Log Expression (RLE).

The RLE is computed for every spot in the array as the difference between the spot and the median of the same spot across all the arrays.（所有数组中同一点的中位数与点之间的差值）

4.4 hierarchical 层次聚类图

qcPlots makes a hierarchical cluster of the samples using the hclust function of the stats package.

• The options for the distance measures are euclidean and pearson.

• 聚类图未按照样本分组聚类的原因：

• 可能是进行层次聚类的基因略少；

• The variables that distinguish the experimental conditions from one another are the differential expressed genes, and that the number of genes may be few relative to the full set of genes of the data set, and hence the cluster analysis will often not reflect the influence of these relevant genes. Therefore if the percentage of differential expression is low, the samples might not be grouped according to their experimental group, since the whole set of genes has very little information regarding the experimental grouping, and the plot will mainly show other grouping features or simply random noise.

4.5 CV图

在探针和基因水平识别重复特征，并计算阵列的变异系数。

The cvArray identifies the replicated non-controlprobes for each feature in the array and computes CV for every microRNA probe set. Then, the median of the CV for each probe set is reported as the array reproducibility.（识别阵列中每个特征的复制的非控制探针，并计算每个microRNA探针集的CV。然后，每个探针组CV的中位数作为阵列重现性报告。）

• To obtain the CV using the cvArray function, we can either choose the MeanSignal or ProcessedSignal. 使用cvArray函数得到CV，我们可以选择 means signal 或 ProcessedSignal。
• A lower CV median indicates a better array reproducibility.
• CV：标准差/均值。

5. 数据标准化(2种方法)

目的：试图补偿芯片之间的系统技术差异，以更清楚地看到样品之间的生物学差异。

AgiMicroRna uses two strategies to obtain a gene signal estimate normalized between arrays.

5.1 AFE算法

• uses the TGS signal processed by the AFE algorithms.
• This TGS can be normalized between arrays using either the quantile (default) or the scale methods.
• If we use the none option the TGS is only log2 transformed. 即：函数tgsNormalization 中有三个标准化参数：“none”(仅做log2转化)，“quantile”，"scale"。
• 如果我们设置 makePLOTpre = TRUE 并且 makePLOTpost = TRUE, 则可以一步出图：
• 密度图（density plot）
• 箱图（boxplot）
• MA图（MA plot）
• RLE图（Relative Log Expression plot (RLE) ）
• 分层聚类图（hierarchical clustering plot）

5.2 RMA算法-affy包

基于affy包的RMA算法：uses the robust multiarray average (RMA) method implemented in the affy package.

RMA通过对该基因的所有探针测量得到对每个基因的表达测量的估计。该算法生成的模型将产生一个考虑到探针效应的基因信号估值。

5.2.1 主要包括三个部分

背景校正

通过拟合normal和exponential 卷积模型到观测强度的失量；

• The normal part represents the background
• the exponential part represents the signal intensities

标准化

使用分位数规范化对 arrays 进行规范化

拟合

用线性模型拟合log2变换后的探针强度（probe intensities）以获得microRNA基因信号的估值。

该模型参数估值使用了median-polish算法；

• Median-Polish算法：用于数据表robust 探索性分析；该方法通常应用于双向表，它通过将行和列标签作为分类因素来拟合数据的附加模型(constant + rows + columns)。

• 该方法由 John Tukey 提出。

  [github.com/borisvish/Median-Polish](https://github.com/borisvish/Median-Polish#:~:text=Median-Polish. Python implementation of Median Polish algorithm for,considering row and column labels as categorical factors.)

5.2.2 rmaMicroRna函数

该函数可以获得每个microRNA信号的RMA估值

ddTGS.rma=rmaMicroRna(dd.micro,normalize=TRUE,background=TRUE)


该函数的具体原理如下：

1. 使用preprocessCore的rma.background.correct函数对信号进行背景校正；
2. 使用limma函数normalizeBetweenArrays用分位数(quantile)方法对信号进行arrays之间的归一化。
3. 对每个基因进行2或4个不同的测量，获得重复探针的中位数。
• 对同一基因的不同探针检测。并且探针测量强度经log2转化，然后通过affy包的rma_c_complete_copy整合成一个单基因测量。
4. rmaMicroRna 将RMA信号估值存储在四个变量中：ddTGS.rma$TGS、ddTGS.rma$TPS, ddTGS.rma$meanS and ddTGS.rma$procS  。

利用RMA处理后的基因信号，可以得到一些图。代码及结果如下：

MMM=ddTGS.rma$meanS
colnames(MMM)=colnames(dd.micro$meanS)
maintitle='TGS.rma'
colorfill='blue'
ddaux=ddTGS.rma
ddaux$meanS=MMM
mvaMicroRna(ddaux,maintitle,verbose=TRUE)


【注】：mvaMicroRna在默认情况下绘制包含在ddaux$meanS中的任何东西，因此，我们需要先创建ddaux对象，然后在ddaux$meanS中存储想使用的矩阵。

在获得全基因信号后，无论是提取AFE产生的TGS，还是使用RMA算法，我们都可以使用AFE算法附加的FLAGS对signals 进行过滤。

结果如下：

rm(ddaux)
RleMicroRna(MMM,"RLE TGS.rma",colorfill)
boxplotMicroRna(MMM,maintitle,colorfill)
plotDensityMicroRna(MMM,maintitle)


6. 质控：过滤探针(probes)

主要利用filterMicroRna函数进行质控；并且质控通常在总基因信号归一化后进行。

• AFE附在每个特征上一个标识，用来识别信号的不同量化错误，可以用来过滤较差质量的microRNAs。
• 该函数返回一个uRNAList，其中包含已过滤的数据。

目的是：

• 去除对照特征；
• 去除未检测到的基因(gIsGeneDetected = 0)；
• 根据阴性对照的表达筛选出表达未达到一定水平的microRNA。

主要过滤标准：

• 通常情况下，根据gIsGeneDetected过滤microRNA。阈值设置为：gisgenedetectflag = 1 即必须保持至少一个实验条件；
• 更高标准：去除gMeanSignal表达接近阴性对照特征的microRNAs。阈值设置为：Mean(negative control expression) + 1.5 x SD(negative control expression)

filterMicroRna函数返回值中包含被过滤的值：

• NOCtrl_FlagIsGeneDetected.txt
• IsGeneDetected Flag for the Non Control Genes, 1 = detected
• IsNOTGeneDetected.txt
• 不被IsGeneDetected Flag检测到的基因

具体示例代码：

ddPROC = filterMicroRna(
    ddNORM,
    dd.micro,
    control=TRUE,
    IsGeneDetected=TRUE,
    wellaboveNEG=FALSE,
    limIsGeneDetected=75,
    limNEG=25,
    makePLOT=FALSE,
    targets.micro,
    verbose=TRUE,
    writeout=FALSE
)


【注】

• 其中 ddNORM 为 通过AFE提取总基因信号估值并通过分位数标准化的数据；
• 也可以将ddNORM替换为RMA产生的ddTGS.rma对象

7. 生成 ExpressionSet 对象

在完成所有处理步骤之后，esetMicroRna创建一个可以用于统计分析的ExpressionSet对象。如果设置makePLOT=TRUE将会展示一些图。(Heatmaps, PCAs)

esetPROC=esetMicroRna(ddPROC,targets.micro,makePLOT=FALSE,verbose=TRUE)


感兴趣的小伙伴也可以阅读文献;

Published: 22 January 2010 Processing of Agilent microRNA array data