高通量蛋白質研究中mzXML文件的探索

06-27

高通量蛋白質研究中mzXML文件的探索

質譜分析是蛋白質組學中常見的一種技術手段,質譜數據處理具有多種策略,多種處理步驟,多種處理軟體的特點.

本文主要探索兩塊內容

通過描述mzXML產生的前世今生來說明一下這種格式產生的原因和這種格式的作用
詳細說明mzXML的數據格式,通過JMZReader的源代碼了解mzXML格式在存儲,解析,搜索等方面的方法和問題

先說說典型的蛋白質鑒定流程如下

我大致的將質譜數據的分析檢測分為四個階段

儀器數據階段

這個階段的數據是指不同的質譜儀器製造商在機器上掃描完以後得到的最原始的數據格式,常見的格式有raw,wiff,t2d,baf,,d,dat,dta,pkl,mgf,MS2等(比如說我們實驗室的儀器產生的格式有raw和wiff兩種),這類格式往往沒有開源文檔手冊可以查詢,如果希望去解析這類格式的文件需要和儀器的製造商單獨聯繫.因此,為了更好的服務各個研究機構,一些蛋白質組學的權威制定組織,例如:

HUPO-PSI(Human Proteome Organization-Proteomics Standards Initiative)

系統生物學研究所(Institute for System Biology,ISB)

歐洲生物信息學研究所(European Bioinformatics Institute,EBI)等

均參與了質譜的標準數據的定製.通過一個維護一個標準的數據格式來統一後續的科學研究,減少重複處理數據帶來的額外工作量.同時為了能夠讓各種原始文件格式轉化為標準格式,科學界也產生了很多相關的轉化工具,而目前為止比較常用的就是msConvert.

那麼再來說一下當年制定標準格式的這件事情到底進展的怎麼樣呢?可以說很曲折

標準原始數據階段

2004年至2008年,用於存儲和交換原始數據的標準主要有mzData和mzXML. 2008年的時候,PSI發布了mzML.目標是可以取代mzData和mzXML.EBI組織還開發了PrideXML這種格式,用於增強搜索結果和實驗信息的存儲.可以說這個階段雖然是全球生物組織一起在制定標準的格式,但是也有些很多互相不兼容和競爭的存在.最終的結果是由ISB開發定製的數據標準mzXML得到了更多的認可.相對的使用也更為廣泛一些(體現在絕大數的質譜數據處理軟體必然都支持mzXML格式)

肽段分析水平階段

通常有PrideXML, pepXML, mzIdentML

蛋白質分析水平階段

通常有mzIdentML和protXML

由於本文主要的探索mzXML格式的相關內容,而在肽段分析水平階段和蛋白質分析水平階段這個格式的使用並不廣泛,因此本文不會過多的涉及肽段分析水平階段和蛋白質分析水平階段的探索.

那麼關於mzXML的格式內容具體是什麼呢?先貼一段給大夥看看

<mzXML xmlns="http://sashimi.sourceforge.net/schema_revision/mzXML_3.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://sashimi.sourceforge.net/schema_revision/mzXML_3.0 http://sashimi.sourceforge.net/schema_revision/mzXML_3.0/mzXML_idx_3.0.xsd"><msRun scanCount="3113" startTime="PT5.0005S" endTime="PT119.982S"> <parentFile fileName="file://PROTEOMICSDEV01/Inetpub/wwwroot/ISB/data/class/day1/runsearch/ICAT_test4.RAW" fileType="RAWData" fileSha1="71be39fb2700ab2f3c8b2234b91274968b6899b1"/> <msInstrument> <msManufacturer category="msManufacturer" value="ThermoFinnigan"/> <msModel category="msModel" value="LCQ Deca"/> <msIonisation category="msIonisation" value="ESI"/> <msMassAnalyzer category="msMassAnalyzer" value="Ion Trap"/> <msDetector category="msDetector" value="EMT"/> <software type="acquisition" name="Xcalibur" version="1.3 alpha 8"/> </msInstrument> <dataProcessing centroided="0"> <software type="conversion" name="ReAdW.exe" version="1.2"/> </dataProcessing> <scan num="19" msLevel="1" peaksCount="1313" polarity="+" scanType="Full" retentionTime="PT353.43S" lowMz="400.39" highMz="1795.56" basePeakMz="445.347" basePeakIntensity="120053" totIonCurrent="1.66755e+007"> <peaks precision="32" byteOrder="network" contentType="m/z-int" compressionType="zlib" compressedLen="782">eJylzskRUEEIBNAJjdAmNELzoBx89VtG5ULRK31+Tp/f5/7aFfjibnDv2ZOrz1779SW9eTPmdMA33unA+4e89/Q0u9hbj3nbP7Pt17f9kfrtae7R6fMf8ZmU61Zvnri95iSdPVuue+M33fT5R7rFU86B3/KTLukT7vwt/jqpP+FJN7uCvsO96ba/kj7xd8G37bz6mjv5a+lLPvP7fM/gBW5/6ttyb+DFky5NBV9x34An3XZv+NbXbHPE9aV8+xMur87cZqe8G3TbH/bq68Crc1KPvLnNnXwzxX7t9bb/VT97+jvwKa/Ak978RtfozNGnzjv11PkedeKppwIub37qu9yp/3Dbn3LtTbj98v6x5aV8ffrv+Z4tT7859rz6K/DFbc6MutRjzgHXJ7/dW06f79Hnlp+pkH/h021/seUddemPtNV525/41D9Ti++13zz1he6GveU64g2eegdvduLtu+Bzpxz/UJfytz82n33i9ovbm/jUZ/4NeAVcnzr/sDf1zKiXT75a8D7fIz+3/cl/4M3xH3P0zb7sLS/pks/8u/DiqT/5nPlH/zaj64DLu1OO/bXgadIf4vJb733k0x/uyUn/iaf/3PZPjj779aVJPZu/4dUV9/CbT968hM/W7x+J9x+3Pf7h6P9XXS93ymvuTVds/eq3fveMPdtfh/v1L/kKuPp0b/kzqWfr1ZfyOuhe/9Jvftqpx7yk17fd6a8Ze7dJuc3edOmfpGt08uKJ93bss9epJTf9kXzFnf7QfwOuz3/sF5e3f7Y9db5n6+1FVwv+v33y9iXdTLOTP+Wbm3Zzm1Potlz5lJd8qe+A69v+cr/me9vrXWz94m55/3Kb673h9hz4pK/z5zF37rvcKeeC2z93s9XdgMunHnl75FNf6m92cSd/0t2Ay2/bvw63Ov8RT/1pzNdXAXf7z4x/vfaYZ27SmedOfv957ZupoDP3wM/WX2z7HHXJ51/2Jp254vqLO/1RbHO2O+ENbs+ma+6z8OLmq0s7/bP9Ze/wFz5tdSlXnz1Ows37AU1Iy5A= </peaks> <scan num="20" msLevel="2" peaksCount="43" polarity="+" scanType="Full" retentionTime="PT356.68S" collisionEnergy="35" lowMz="223.089" highMz="531.078" basePeakMz="428.905" basePeakIntensity="301045" totIonCurrent="764637"> <precursorMz precursorIntensity="120053">445.3466797</precursorMz> <peaks precision="32" byteOrder="network" contentType="m/z-int" compressionType="zlib" compressedLen="55">eJw7wAABDQzYAbr4ARzyB9DoBhw0OnDAYQ66uAOaPDqNbh6x7kC3DxeA6UO3D10cxgcA3f0RQQ== </peaks> </scan> </scan></msRun><index name="scan"> <offset id="19">1209</offset> <offset id="20">2577</offset></index><indexOffset>31713417</indexOffset><sha1>e83e234ed25a2e675ad8a3fbcd56f16585a237c2</sha1></mzXML>

上面的mzXML文件是截取的實際文件中的一小段實例. 而實際上做一次質譜實驗產生的原始數據(以wiff格式為例),大小在1.5G左右,而經過msConvert轉化為mzXML格式以後的大小大概在30G左右.這僅僅只是一次實驗的數據量而已.

我們看一下這個文件的格式內容,先給出mzXML的標準數據格式的官方說明文檔

http://www.peptideatlas.org/tmp/mzML1.1.0.html?

www.peptideatlas.org

這個是mzXML 1.1.0版本的格式詳細說明,其實目前的mzXML的版本號已經到3.0了,不過貌似主體核心的結構沒有變化,加上這個標準的文檔維護者應該有點力不從心.因此後續也就沒有相關的mzXML的文檔的更新,大家看看這個1.1.0版本的文檔也就差不多了.

首先要講一下這個mzXML文檔裡面大概包含的內容.

在Swaths-MS(一種DIA的質譜分析手段)的實驗中,每一個mzXML的文件基本上都是10G左右的大小.通過傳統的DOM解析或者試圖將這個文件一把載入到內存中來進行搜索和數據處理的話都不太實際.當然作為科學實驗室你當然可以買一台20G內存的伺服器來搞.

mzXML包含了一次實驗中所有的分子片段的質譜圖,並且包含了一些實驗的基本數據.其核心數據是峰譜圖,峰譜數據也是佔了整個文件最大頭的一塊數據.大家看到的<peaks>這個標籤內的一段很長的數據就是峰譜數據,其本身是荷質比與強度的Key-Value鍵值對

關於其他的一些諸如msManufacturer這類的標籤我就不說了,大家直接看上面給到的官方說明文檔就可以了解到.我要說的是如何解析這個文件.這個也是mzXML在格式上比較特殊的方法.當然並不是說mzXML這種格式的存儲方式就是最優的,相信大家看完我的解析以後會對蛋白質質譜數據的存儲方法有很多的想法.

首先一次試驗得到的掃描<scan>標籤的記錄大概是100000條以上.行數在數百萬行的量級.因此如果我想查詢其中的某一條scan記錄的話,用一次性載入到內存中的方式顯然是不可行的.這個和使用的計算機語言無關,和計算機硬體有關係.

我們發現mzXML文件的末尾有一大塊叫index的標籤

<index name="scan"> <offset id="19">1209</offset> <offset id="20">2577</offset></index>

這個其實就是mzXML對於方便搜索的一個優化表達. index下標記的是哪一個scan對應在文件中的起始位置給標記出來.比如說第一條就表示id是19的scan是在整個文件的第1209個位元組處開始.然後需要使用一種隨機文件訪問的技術來直接讀取某一個大文件的某一個位置.這樣就不需要將文件完整的讀入到內存中再做搜索了.我們稱這塊數據叫做索引區

而最後還有一個

<indexOffset>31713417</indexOffset>

這個標籤標明了索引區的起始位置(因為你會發現起始索引區也有十萬行以上).而這個indexOffset我們可以通過直接從文件尾部讀取,進行逆向掃描來查找.

public void parseIndexOffset() throws IOException { RandomAccessFile rf = getRandomAccessFile(); long position = rf.length() - 1; rf.seek(position); byte words[] = new byte[1]; for (int i = 0; i < 1000; i++) { rf.read(words); if (words[0] == ) { //detect the line break symbol String line = rf.readLine(); if (line != null && line.contains("<indexOffset>")) { line = line.trim().replace("<indexOffset>", "").replace("</indexOffset>", ""); setIndexOffset(Long.valueOf(line)); return; } } rf.seek(--position); } }

用Java寫了一個實現,可以快速的獲取到尾部的索引塊的坐標值.

拿到了索引值以後,我們可以通過索引值來順序讀取並且解析所有的scan標籤.

Scan標籤內部最重要的一個數據就是peaks標籤.每一個scan可以包含多個scan.而每一個scan的根目錄下一般都只有一個peaks標籤.那麼peaks裡面的那麼長一段代碼是什麼東西呢?

質譜數據最重要的三個數據是 RT時間,m/z質荷比,離子丰度

而每一個scan標籤下都是某一個時刻的離子質荷比與丰度的Key-Value對(有可能有幾十萬對),

<peaks precision="32" byteOrder="network" contentType="m/z-int" compressionType="zlib" compressedLen="55">eJw7wAABDQzYAbr4ARzyB9DoBhw0OnDAYQ66uAOaPDqNbh6x7kC3DxeA6UO3D10cxgcA3f0RQQ==</peaks>

這個peaks標籤裡面會指定壓縮的float數據精度是多少(32位或者64位),壓縮方式是什麼(一般是zlib).只要按照這種規格解碼就可以查看到真實的離子質荷比與丰度的相關數據了.

這就是mzXML的數據格式的相關內容.看完以後大家肯定會有很多疑問.我作為一個軟體工程師,在初次接觸到這種樣子的格式以後也感覺非常的無語.30G大小的單個文件用以保存數據,甚至在XML中加入索引的方式來優化搜索能力.那為啥不直接存在資料庫呢?(事實上有人這麼做過,而且效果也挺好的),因此我也自己寫了一堆代碼並且把所有的數據全部格式化以後存儲到Mongodb資料庫中,最後的結果顯示數據最終壓縮到了原來的85%左右,雖然壓縮比不高,同時因為需要寫入資料庫導致最終的解析速度比較慢,但是一旦存入資料庫以後,做相關的SQL查詢就快了很多.還是有很多好處的.

從官方的JMZReader的開源代碼中,我觀察到JMZReader是先通過暴力方法遍歷所有的索引位置(而不是使用mzXML自帶的indexOffset),JMZReader通過一個叫xxindex的小工具來建立索引,由於演算法的暴力遍歷,因此速度不快.不過xxindex建立的索引相對mzXML自帶的索引來說要強大許多.在索引建立完畢以後的優勢還是不錯的.