失眠多梦脱发心悸？物质鉴定拯救你的茫然！
BioDeep公共免费物质注释平台MetAnno目前已经支持MetaDNA算法，欢迎老师试用：http://msms.biodeep.cn/new_task.php?app=new_task
MetaDNA文献引用以及BioDeep平台物质注释程序算法原理请参见：http://www.metabolomics-shanghai.org/publications.php
物质鉴定的过程实际上就是一个通过碎片相似度来查找meta注释信息的一个过程。
假设我们所认知的物质是由两部分信息构成的：
1. meta.db 物质注释信息：名字，质量，分子式，数据库编号，等
2. lib_ms2 物质结构信息：二级质谱图碎片(mz, intensity)
程序包通过lib_ms2和样本的peak_ms2进行二级比对，然后通过一一对应关系得到lib_ms2所对应的meta.db部分的数据。
上面的步骤即是一个完整的未知代谢物物质注释的过程。
如何评判两个二级质谱图相似
说道相似度的计算，对于大部分的生物信息学的老师可能会首先想到的是Smith–Waterman或者Needleman-Wunsch这类用于序列相似度计算的动态规划的算法，当然这两个算法是可以应用于质谱的二级碎片甚至是三级或者四级碎片上面的相似度的计算的。当然，对于序列相似度老师您也可以使用Pearson或者Spearman相关度来进行评价。但是对于质谱图之间的比较，我们使用的是一个来自于朱正江老师的论文之中的更加简单的算法：质谱图比对公式(SSM)。
质谱图比对公式有些类似于余弦相似度(cosA)的计算，只不过cosA的结果值是位于[-1, 1]这个区间内的，而SSM的结果值则是位于[0, 1]这个区间内，通过SSM相似度可以比余弦相似度更加直观的考察两个质谱图之间的相似度的高低。
为了容易理解SSM算法，我们可以首先从一个二维平面上面的cosA结果值的计算来理解：
我们已经明白，两个向量间的余弦值可以通过使用欧几里得点积公式求出：
a.b = ||a||*||b||cosA
cosA = a.b / (||a||||b||)给定两个属性向量，A和B，其余弦相似性θ由点积和向量长度给出：
cosA = func(a, b) sum(a*b) / sqrt(sum(a^2)) * sqrt(sum(b^2))
由于cosA的结果值是位于[-1, 1]这个开区间内，并且我们将向量限制在第一象限范围内（实际上二级质谱图的碎片的响应强度都是正实数，都位于第一象限内)。故而cosA的结果值在物质鉴定程序之中是位于[0, 1]这个开区间内的。
在这里，为了便于理解，我们假设某一个待鉴定的物质只有两个碎片，并且各个碎片的响应强度都被归一化到了(0, 1]这个半开区间内，并且假设标准品库reference参考是上面的示例图之中的X坐标轴，A与B则分别是样本原始数据之中的两个二级feature的时候，可以计算出下面的cos夹角结果：
对于feature A，其坐标为(0.8, 0.05)，则与reference参考之间的夹角为：
# [1] 0.9999219
对于feature B，其坐标为(0.2, 0.98)，则与reference参考之间的夹角为：
# [1] 0.19996
通过上面的例子，老师您也可以明白，对于两个向量而言，如果二者的夹角越小，则cosA值越接近于1，反之二者的夹角越接近于90度的时候，则cosA值越接近于0，这样子我们就可以评判出两个质谱图之间的相似性程度了。上面的两个简单的cosA计算的例子大致就是SSM算法的基本原理了。当然，质谱图并不都是像例子之中的一样只有两个碎片，都是位于一个二维平面内的，高分辨率下的QE平台产生的HCD模式的质谱图可能有好几百个碎片，则我们只需要同理，将上面的计算过程从两个元素长度的向量替换为n个元素长度的向量就好了。
因为cosA的值可能会位于[-1, 1]这个区间内，所以SSM算法对余弦相似度做了一些小修改，结果值就位于[0, 1]范围区间内，这样子就更加容易表述相似度这个概念：0完全不相似，1表示完全相似。修改后得到的SSM计算公式为：
#' SSM alignment score
#'
#' @param q The query sample spectra matrix
#' @param s The library spectra matrix
#'
#' @return A single score value in range \code{[0, 1]}.
#'
#' @details The \code{q} and \code{s} spectra matrix should contains
#' two column: \code{mz}, \code{into}. And the row number of this
#' two spectra matrix should be the same, and the orders of the mz
#' raw should match each other.
#'
SSM <- function(q, s) {
q <- q[ ,2];
s <- s[ ,2];
score <- sum(q * s) / sqrt(sum(q ^ 2) * sum(s ^ 2)); if (is.nan(score) || is.na(score) || score == Inf || score == -Inf) { 0;
} else {
score;
}
}ExtractMzI
老师您可能会发现，不论是哪种相似度的计算方法，都会要求样本的向量和数据库之中的参考值是等长的，所以在计算相似度之前，我们就还需要完成一个操作：ExtractMzI，即将样本之中的数据和数据库之中的参考之间进行补齐。老师您可以使用下面的一段来自文献的代码进行质谱图碎片的补齐操作：
#'
#' @description Extract m/z and intensity from scan, applied on the extraction
#' of both precursor ions and fragment ions
#'
#' @param scan.all all scans of ms1 or ms2
#' @param mzlib mz of targeted ions in library
#' @param tolerance The given tolerance, usually is tolerance with da 0.3. This tolerance
#' parameter is a list object that from \code{MetaDNA} package.
#'
#' @details \code{scan.all} is a matrix of the ms2 spectra with at least two column:
#' \code{mz}, \code{into}. The \code{mzlib} is the mz vector that extract from#' subject.
#'
#' @return This function returns a spectra matrix that with row subset extract from#' \code{scan.all} spectra matrix by match mz with \code{mzlib}.
#'
ExtractMzI <- function(scan.all, mzlib, tolerance = MetaDNA::tolerance(0.3, "da")) {
scan.all <- rbind(NULL, scan.all);
colnames(scan.all) <- NULL;
scan.all <- as.matrix(scan.all);
assert <- tolerance$assert;
scan.mz <- as.vector(unlist(scan.all[, 1]));
scan.into <- as.vector(unlist(scan.all[, 2]));
into <- sapply(mzlib, function(mz) {
# Get m/z row match in scan.all spectra matrix
mzi <- assert(scan.mz, mz);
if (sum(mzi) == 0) {
# no match
0;
} else {
mz.into <- scan.into[mzi];
# returns the highest intensity value for match
max(mz.into);
}
});
into <- as.numeric(into);
cbind(mzlib, into);
}
补齐的原理其实就是：我们只需要以标准品库之中的质谱图为参考模板，在一定的误差范围内，从query之中选出对应的mz的碎片的响应强度用以构成一个和参考等长的计算向量即可。如果参考模板上面的某一个mz的碎片在query之中是缺失的，则将其相对应的响应强度设置为零即可。
构建标准品库的MetaDB信息
从这篇小文章的最开始我们了解到，物质注释的过程大致就是一个根据相似度找对应的注释信息的过程，那么在进行了上面的相似度比较之后，我们只需要得到对应的注释信息就可以完成一个简单的注释过程了。为了批量的构建物质注释过程，我们需要在本地构建一个物质注释信息库。帕诺米克的BioDeep标准品库之中的注释信息都是来自于NCBI的PubChem数据库之中的compound物质信息，老师您可以从pubchem的ftp服务器下载这些注释信息：
ftp://ftp.ncbi.nlm.nih.gov/pubchem/Compound/CURRENT-Full
在这个数据库之中，老师您可以找到几乎每一种物质的分子式，数据库编号，分子结构信息等用于化学信息学计算所需要的原始数据。
帕诺米克免费公共物质注释平台
最后，我想给老师介绍一下帕诺米克于近期推出的一个免费的公共物质注释平台：免费公共物质注释平台 MetAnno 项目目前已经完工，老师您可以访问 http://msms.biodeep.cn/ 进行使用。
其它文章：