一 引言

1.1編寫目的

進(jìn)行該測(cè)試以及撰寫此報(bào)告有以下幾個(gè)目的

1.通過對(duì)測(cè)試結(jié)果的分析，得到對(duì)軟件質(zhì)量的評(píng)價(jià)；

2.分析在Illumina測(cè)序平臺(tái)下，MapSplice能夠獲得最大junction數(shù)目以及mapping率的參數(shù)；

3.分析在ionproton測(cè)序平臺(tái)下，MapSplice能夠獲得最大junction數(shù)目以及mapping率的參數(shù)；

4.嘗試找到參數(shù)與測(cè)序長度的經(jīng)驗(yàn)性關(guān)系。

1.2背景

MapSplice是一個(gè)RNA-seq數(shù)據(jù)分析工具，其核心程序是bowtie.可以快速的確認(rèn)exon-exon剪切拼接。主要功能和Tophat差異不大。

與Tophat不同的是，MapSplice并沒有針對(duì)某一種測(cè)序平臺(tái)而開發(fā)，所以對(duì)于75bp以下的短序列以及75bp以上的長序列reads都可以使用。目前，全球最大的癌癥研究項(xiàng)目TCGA（The Cancer Genome Atlas）正在主要推崇使用這個(gè)軟件。

Ionproton屬于二代測(cè)序中較新的平臺(tái)，可以認(rèn)為是二點(diǎn)五代測(cè)序平臺(tái)，其測(cè)序長度平均在100個(gè)bp以上。目前我們公司使用的就是這個(gè)平臺(tái)的進(jìn)行二代測(cè)序分析。

鑒于之前使用Tophat進(jìn)行參數(shù)優(yōu)化以后發(fā)現(xiàn)結(jié)果并不是很理想，所以決定跟換軟件進(jìn)行測(cè)試，尋找更好的結(jié)果。因此，提出此次工作內(nèi)容，探索更好的參數(shù)配置，提高mapping率以及junction數(shù)目。

1.3用戶群

主要讀者：公司研發(fā)部，公司管理人員。

其他讀者：項(xiàng)目及銷售相關(guān)人員。

1.4 數(shù)據(jù)對(duì)象：

1.5 測(cè)試階段

軟件測(cè)試

1.6測(cè)試工具

Samtools version：0.1.18；

IGV version：2.3.18；

Awk；

1.7 參考資料

MapSplice userguide

Wang K, Singh D, Zeng Z, et al. MapSplice: accurate mapping of RNA-seq reads for splice junction discovery[J]. Nucleic acids research, 2010, 38(18): e178-e178.

Trapnell C, Pachter L, Salzberg S L. TopHat: discovering splice junctions with RNA-Seq[J]. Bioinformatics, 2009, 25(9): 1105-1111.

二 測(cè)試概要

關(guān)于MapSplice參數(shù)測(cè)試從2013年9月19日開始到2013年9月26日結(jié)束，共持續(xù)7天，一共25個(gè)測(cè)試用例。

主要測(cè)試內(nèi)容如下：

1. 軟件安裝以及依賴性測(cè)試。

2. 文件分割以后查找junction數(shù)目以及不進(jìn)行分割查找junction數(shù)目的差異大小，能否接受，為今后并行化文件回帖提供依據(jù)。

3. Segment參數(shù)進(jìn)行優(yōu)化工作。

4. 針對(duì)Illumina測(cè)序平臺(tái)數(shù)據(jù)以及ionproton測(cè)序平臺(tái)數(shù)據(jù)的mapping能力差異。

5. 簡要測(cè)試MapSplice檢測(cè)融合基因的能力

2.1工作計(jì)劃進(jìn)展

測(cè)試內(nèi)容	計(jì)劃開始時(shí)間	實(shí)際開始時(shí)間	計(jì)劃完成時(shí)間	實(shí)際完成時(shí)間	工作完成情況
軟件安裝	2013年9月19日	2013年9月19日	2013年9月19日	2013年9月23日	本地安裝受阻，服務(wù)器端安裝正常。
軟件依賴性查找	2013年9月24日	2013年9月24日	2013年9月24日	2013年9月24日	順利
不同測(cè)序平臺(tái)回帖能力	2013年9月24日	2013年9月24日	2013年9月24日	2013年9月24日	順利
文件分割與否回帖差異	2013年9月25日	2013年9月25日	2013年9月25日	2013年9月25日	順利
Segment參數(shù)優(yōu)化	2013年9月26日	2013年9月26日	2013年9月26日	2013年9月26日	順利
融合基因檢測(cè)	2013年9月26日	2013年9月26日	2013年9月27日	2013年9月27日	順利

2.2測(cè)試執(zhí)行

此次測(cè)試嚴(yán)格按照項(xiàng)目計(jì)劃和測(cè)試計(jì)劃執(zhí)行，按時(shí)完成了測(cè)試計(jì)劃規(guī)定的測(cè)試對(duì)象的測(cè)試。針對(duì)測(cè)試計(jì)劃制定規(guī)定的測(cè)試策略，依據(jù)測(cè)試計(jì)劃和測(cè)試用例，將網(wǎng)絡(luò)數(shù)據(jù)以及我們觀測(cè)的關(guān)鍵參數(shù)進(jìn)行了完整的測(cè)試。

2.3測(cè)試用例

2.3.1功能性

1.測(cè)試主要實(shí)現(xiàn)，包括較高的mapping率以及較多的junction數(shù)目。

2.測(cè)試junction數(shù)目與文件分割與否的相關(guān)性大小。

三 測(cè)試環(huán)境

3.1軟硬件環(huán)境

四 測(cè)試結(jié)果

4.1 軟件安裝

安裝中，我們使用的軟件版本是MapSplice 2.1.5。在本地進(jìn)行測(cè)試的時(shí)候由于當(dāng)時(shí)未知的軟件依賴關(guān)系，并沒有安裝成功。軟件提示報(bào)錯(cuò)為本地bowtie沒有在系統(tǒng)中找到，于是在本地安裝了與軟件要求對(duì)應(yīng)的bowtie 0.12.7 。本地可以使用bowtie，但是MapSplice仍然報(bào)這個(gè)錯(cuò)誤，于是放棄在本地進(jìn)行安裝。在服務(wù)器的安裝很順利，很快就測(cè)試通過。

4.2文件分割mapping與未分割mapping

進(jìn)行文件分割運(yùn)行的最主要的考慮是為了嘗試能否進(jìn)行分布式的計(jì)算，所以我們?cè)谶@一部分的工作中將文件分割成4份分開進(jìn)行運(yùn)算，然后將這4個(gè)文件運(yùn)行出的junction數(shù)目相加比較與未分割情況下的junction數(shù)目差異。為了得到更加準(zhǔn)確的效果，在本次測(cè)試中，我們使用了3個(gè)測(cè)序深度的ionproton測(cè)序平臺(tái)得出的reads，分別是20萬個(gè)reads，200萬個(gè)reads以及整個(gè)文件（一共33926644個(gè)reads）進(jìn)行分析。文件統(tǒng)一分割為4個(gè)文件。測(cè)試結(jié)果如下：

表 20萬個(gè)reads運(yùn)行所得結(jié)果及差值

表 200萬個(gè)reads運(yùn)行所的結(jié)果及差值

表 所有reads運(yùn)行所有的結(jié)果及差值

通過上述結(jié)果可以知道分開以后與未分開時(shí)相差大約5%以上（占未分開的junction數(shù)目）所以可以認(rèn)為并不是很適合將reads分開以后進(jìn)行mampping。

4.3 segment參數(shù)優(yōu)化

由于在Tophat參數(shù)探索的過程中，我們通過分析發(fā)現(xiàn)在所有參數(shù)中，segment_length是對(duì)junction影響最顯著的參數(shù)，所以我們?cè)趯?duì)MapSplice進(jìn)行分析時(shí)，主要也是分析這個(gè)參數(shù)。在測(cè)試過程中我們發(fā)現(xiàn)當(dāng)這個(gè)參數(shù)大于30的時(shí)候就會(huì)報(bào)錯(cuò)，而軟件參數(shù)中segment_length的下限值限定為18，軟件說明中推薦對(duì)于50bp以上的reads文件，建議使用25這個(gè)長度，而根據(jù)文獻(xiàn)中算法設(shè)計(jì)的思路可以知道，當(dāng)這個(gè)數(shù)據(jù)越大的時(shí)候，整個(gè)junction的敏感度就會(huì)越低，而對(duì)應(yīng)的程序運(yùn)行時(shí)間就會(huì)越短，與之相反，當(dāng)這個(gè)數(shù)據(jù)越小的時(shí)候，整個(gè)junction的敏感度就越高，而對(duì)應(yīng)程序的運(yùn)行時(shí)間就會(huì)越長。在本次測(cè)試中，我們對(duì)這個(gè)數(shù)值從18到28進(jìn)行了抽樣實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如下：

表 segment_length數(shù)值與junction數(shù)目關(guān)系表（以上ratio省略%）

從上述表格中可以看出對(duì)于不同的segment_length而言，junction數(shù)目的百分比確實(shí)是有變化的，總體趨勢(shì)是segment_length越長，junction的數(shù)目就越少，由于RNA-seq回帖率與測(cè)序深度正相關(guān)的關(guān)系，我們可以推測(cè)對(duì)于更多數(shù)目的數(shù)據(jù)而言，這個(gè)數(shù)值會(huì)有提高。在數(shù)據(jù)記錄中，我們同時(shí)也記錄了任務(wù)運(yùn)行的總時(shí)間，與文獻(xiàn)符合的是，segment_length長度越短，運(yùn)行時(shí)間就會(huì)越低，而且我們發(fā)現(xiàn)時(shí)間增長的速度是很夸張的。當(dāng)segment_length為28的時(shí)候，運(yùn)行時(shí)間是12分鐘12秒，而當(dāng)segment_length為18的時(shí)候，運(yùn)行時(shí)間是1小時(shí)7分45秒。可以看出這個(gè)時(shí)間差是很大的。綜考慮我們認(rèn)為如果有需要，取20到22都是不錯(cuò)的選擇。

4.4融合基因檢測(cè)參數(shù)測(cè)試

本實(shí)驗(yàn)中，我們主要檢測(cè)了檢測(cè)融合基因以及檢測(cè)junction之間的關(guān)系。我們的檢測(cè)了在尋找融合基因情況下，junction數(shù)目的變化，全部結(jié)果如下所示：

指定參數(shù)non_canonical_fusion	500000
seg_length	18	19	20	22	28
junction_numbers	39786	40269	40203	39450	未測(cè)試
threads	20	20	23	20	未測(cè)試
time(min)			13:34
按照>10KB都過濾掉
	18	19	20	22	28
junctions_filted	36952	37399	37361	36621	未測(cè)試
junction_numbers	39786	40269	40203	39450	未測(cè)試
conseved_ratio	92.87689	92.87293	92.93088	92.8289
ratio	7.3904	7.4798	7.4722	7.3242
按照<20bp的都過濾掉
	18	19	20	22	28
junctions_filted	39775	40257	40191	39440	未測(cè)試
junction_numbers	39786	40269	40203	39450	未測(cè)試
conseved_ratio	99.97235	99.9702	99.97015	99.97465
ratio	7.955	8.0514	8.0382	7.888
按照<20bp以及>10Kb的都過濾掉
	18	19	20	22	28
junctions_filted	36941	37387	37349	36611	未測(cè)試
junction_numbers	39786	40269	40203	39450	未測(cè)試
conseved_ratio	92.84924	92.84313	92.90103	92.80355
ratio	7.3882	7.4774	7.4698	7.3222
mapping_ratio	41.7	41.8	41.54	40.92	未測(cè)試
real_junction_ratio	17.71751	17.88852	17.98219	17.89394	未測(cè)試

由于之前的測(cè)試，我們考慮的參數(shù)中已經(jīng)放棄了segment_length等于28這個(gè)參數(shù)，所以在這一步中，為了節(jié)約計(jì)算資源，我們并沒有計(jì)算segment_length等于28情況下的測(cè)試數(shù)據(jù)。從上表中可以很明顯的看出當(dāng)檢測(cè)融合基因時(shí)，整體數(shù)據(jù)的mapping率明顯下降。因此導(dǎo)致的real_junction_ratio數(shù)目的提升并不能認(rèn)為可能是真的提升。

五．測(cè)試結(jié)論與討論

5.1平臺(tái)差異

通過查閱已經(jīng)有的資料，我們知道Illumina測(cè)序平臺(tái)和ionproton平臺(tái)最直觀的差別在于后者的平均測(cè)序長度比前者長；在我們測(cè)試的例子中，Illumina的測(cè)序長度在50-97個(gè)bp之間，而ionproton的測(cè)序長度在50到235個(gè)bp之間。從此可以看出兩者的最合適參數(shù)應(yīng)該是有差別的。通過上一次tophat與這一次MapSplice的比較，我們發(fā)現(xiàn)，無論如何提高tophat的參數(shù)，我們都很難接近MapSplice使用默認(rèn)參數(shù)下的junction數(shù)目，所以我們認(rèn)為對(duì)于公司ionproton測(cè)序平臺(tái)，我們使用MapSplice會(huì)更加適合。而在我們的測(cè)試結(jié)果中，對(duì)于Illumina測(cè)序平臺(tái)測(cè)試時(shí)，進(jìn)行單端實(shí)驗(yàn)的結(jié)果如下：

表 ionproton與Illumina單端結(jié)果計(jì)算junction百分?jǐn)?shù)

5.2文件分割測(cè)試

通過這個(gè)測(cè)試的結(jié)果，我們可以看出分割前后運(yùn)行得出的junction數(shù)目差距為5%（相比未分割的情況）以上，并且這個(gè)數(shù)目隨著我們的測(cè)序深度的提高而提高。所以從這個(gè)結(jié)果而言，我們認(rèn)為不適合將文件分割進(jìn)行處理。

5.3segment參數(shù)測(cè)試結(jié)果

在測(cè)試實(shí)驗(yàn)中，我們發(fā)現(xiàn)segment_length參數(shù)從28向18變化的過程中，總體趨勢(shì)是由少變多變少，整體趨勢(shì)圖如下：

圖 segment_length與junction_ratio關(guān)系圖

從上圖中可以看出大約在20到22的時(shí)候是最好的。在官方說明文檔中，作者推薦當(dāng)序列長度大于50的時(shí)候推薦使用參數(shù)25。下表是segment_length與測(cè)試時(shí)間之間的關(guān)系：

圖 segment_length與running_time關(guān)系圖

上圖中，我們計(jì)算時(shí)間是使用實(shí)際計(jì)算總時(shí)間乘以運(yùn)行的CPU數(shù)目。其中在22這個(gè)長度上時(shí)運(yùn)行的CPU數(shù)目是20個(gè)，所以時(shí)間有所波動(dòng)，總體而言來看，在長度為20到28之間時(shí)間變化還是可以接受的，然后當(dāng)長度繼續(xù)下降的時(shí)候，時(shí)間就開始指數(shù)級(jí)的上升的，這一點(diǎn)可以從圖中看出。

所以，綜合取舍junction率以及運(yùn)行時(shí)間，我們認(rèn)為使用默認(rèn)參數(shù)是可以接受的，但是使用20至22也許會(huì)有更好的結(jié)果。

5.4融合基因檢測(cè)參數(shù)設(shè)置

在我們的測(cè)試數(shù)據(jù)中，我們可以很明顯的看出在各個(gè)segment_lengh情況下，mapping率都有下降，相比不做這一步檢測(cè)，mapping率下降了至少30%，我們一開始認(rèn)為是把部分junction的數(shù)據(jù)被認(rèn)為是融合基因，當(dāng)我們檢測(cè)的時(shí)候才發(fā)現(xiàn)實(shí)際情況與我們的預(yù)測(cè)是不符合的。軟件找到的融合基因數(shù)目十分少，并且基本都是跨染色體的。因此我們提出了新的想法，程序在同時(shí)進(jìn)行查找junction以及融合基因的時(shí)候，為了確保計(jì)算時(shí)間不會(huì)超過單查找junction時(shí)的時(shí)間太多，并且由于查找融合基因是比較消耗計(jì)算資源的，所以程序在查找junction的時(shí)候并沒有分配過多的資源，導(dǎo)致了更多的reads沒有被程序mapping上去，因此我們?cè)诖苏J(rèn)為實(shí)際應(yīng)用中，應(yīng)該將查找融合基因以及查找junction分成兩步分開進(jìn)行，如何能夠使得兩步的資源能更加節(jié)省，將是我們接下來的工作。

5.5 測(cè)試中的問題

在測(cè)試過程中我們發(fā)現(xiàn)了一個(gè)有趣的情況，如下圖所示：

圖 不同segment參數(shù)下查找junction數(shù)目的能力

上圖中，首先可以看到在參考基因組中這個(gè)部分是有junction的，中間的四個(gè)條帶從上到下依次是長度為18,19,22,28四個(gè)參數(shù)情況下的對(duì)應(yīng)這個(gè)位置的回帖情況，可以很清楚的看到，在參數(shù)為18,28的時(shí)候是找到了這個(gè)區(qū)域的，但是在中間參數(shù)19,22的情況下，并沒有回帖到這個(gè)位置，由于這個(gè)部分并不是很短，我們可以認(rèn)為這個(gè)部分在染色體上是唯一的，所以排除了這兩個(gè)參數(shù)情況下回帖到其他地方的可能性，確定這個(gè)部分對(duì)應(yīng)的reads并沒有回帖上去。因?yàn)椴⒉皇钱?dāng)segment_length變小或變大的情況下才逐漸出現(xiàn)的，所以可以認(rèn)為是隨機(jī)的結(jié)果，這暗示我們?nèi)绻骋淮蔚慕Y(jié)果不是很理想的情況下可以通過重復(fù)或更改參數(shù)重復(fù)來提高junction數(shù)目。

另外在測(cè)試中，我們找到了支持segment_length越短，查找的敏感性就越明顯的圖像證據(jù)，如下圖所示：

圖 不同segment_length下查找junction的敏感性

上圖中可以看出，在參考基因中，這個(gè)部分是有junction的，而在segment_length為19,22,28的時(shí)候，都沒有找到回帖上，我們認(rèn)為這個(gè)結(jié)果對(duì)于文獻(xiàn)中提到的segment_length越短，敏感性越強(qiáng)這個(gè)說法。

六．測(cè)試總結(jié)

1.由于MapSplice在我們已經(jīng)配置好的服務(wù)器上能夠很流暢的直接使用，所以對(duì)于我們的hdfs而言，我們認(rèn)為可以直接裝配使用。對(duì)于本地的軟件使用的可能需要復(fù)雜的軟件支持，由于在這一步我們花費(fèi)了部分的時(shí)間，所以在此并沒有進(jìn)行詳細(xì)的尋找軟件依賴關(guān)系。

2.綜合考慮junction查找能力以及運(yùn)行時(shí)間，我們認(rèn)為在一般情況下，默認(rèn)參數(shù)就是可以的了。當(dāng)有特殊需求時(shí)，可以考慮使用參數(shù)在20到22內(nèi)的任意值。

3.對(duì)于Illumina的單端數(shù)據(jù)而言，我們認(rèn)為使用Tophat的效果比使用MapSplice的效果好，對(duì)于Illumina的雙端數(shù)據(jù)而言， 對(duì)于ionproton的數(shù)據(jù)而言，我們認(rèn)為使用MapSplice的效果遠(yuǎn)比使用Tophat好，不論是mapping率還是junction數(shù)都顯示使用MapSplice更加合適。

4.我們認(rèn)為查找junction以及查找融合基因這兩個(gè)工作應(yīng)該分開進(jìn)行。由于時(shí)間關(guān)系，我們并沒有查找弱化junction數(shù)目查找情況下，對(duì)融合基因查找的影響。

5.鑒于MapSplice查找junction時(shí)存在一定幾率不能找全所有的junction，所以對(duì)于查找情況不好的數(shù)據(jù)，我們可以通過簡單的重復(fù)運(yùn)行或更改參數(shù)運(yùn)行來嘗試提高這個(gè)數(shù)據(jù)。

6.我們測(cè)試結(jié)果顯示回帖操作并不能通過將源文件分割分別回帖來實(shí)現(xiàn)分布式運(yùn)行。

七．測(cè)試中使用的命令，參數(shù)及說明

測(cè)試的結(jié)果在/media/hdfs/nbCloud/public/test/Illuminaandionproton0906/MapSplice-v2.1.5文件夾中下。

測(cè)試中統(tǒng)計(jì)junction數(shù)目的命令為

awk -F"\t" '{if((($3-$2)>20)&&(($3-$2)<10000)){total+=$5}}END{print total}' ./split4_test_segment_length20_non_canonical/junctions.txt

分割文件使用的perl文件見附件

使用MapSplice命令如下

Python MapSplice1.py  -c /media/hdfs/nbCloud/public/nbcplatform/genome/human/hg19_GRCh37/ChromFa/sep \

-x /media/winE/genome/human/hg19_GRCh37/ChromFa/all/hg19_GRCh37_bowtie_index -o ./split4_test_segment_length25/ -1  ../split4.fq  -p 10 -s 25

Mapsplice參數(shù)說明

其中重要的參數(shù)是粗體表現(xiàn)。

必須參數(shù)：

-c 序列文件的文件夾，注意：文件必須是fasta格式，后綴是.fa文件。

-x       bowtie_index指定的路徑及前綴。注意：只支持bowtie1的索引，并不支持bowtie2的索引。如果沒有設(shè)定這個(gè)選項(xiàng)，或者指定的路徑?jīng)]有對(duì)應(yīng)的索引，則會(huì)在結(jié)果輸出路徑下自動(dòng)建立索引。

-1 FATSA格式或者是FASTQ格式。對(duì)于雙端的回帖，這對(duì)應(yīng)編號(hào)為1的文件。多個(gè)文件用逗號(hào)隔開，文件名之間不能有空格

-2 FATSA格式或者是FASTQ格式。對(duì)于雙端的回帖，這對(duì)應(yīng)編號(hào)為2的文件，并且兩個(gè)文件順序必須一致。多個(gè)文件用逗號(hào)隔開，文件名之間不能有空格

-p/--threads 線程數(shù)目，默認(rèn)是1；

-o/--output 指定Mapsplice輸出文件夾，默認(rèn)是./mapsplice_out/這個(gè)文件夾 沒有寫清楚輸出文件的具體樣式，譬如輸出文件前綴，文件名，輸出文件類型

--qual-scale 輸入文件的打分類型。默認(rèn)是自動(dòng)尋找，可以指定如下：phred33,phred64,solexa64

--bam 默認(rèn)的輸出文件時(shí)SAM格式的文件，通過這個(gè)選項(xiàng)可以指定輸出BAM文件。

--keep-tmp 保存中間文件。

-s/--seglen 指定segment_length,通常默認(rèn)是25，我們測(cè)試的結(jié)果暗示這個(gè)結(jié)果在20到22都是不錯(cuò)的。最小值是18，目前測(cè)試的結(jié)果暗示最大值不要超過30.

--min-map-len 軟件只會(huì)記錄完全匹配或者匹配數(shù)目不小于這個(gè)參數(shù)的序列。默認(rèn)參數(shù)是50.

-k/--max-hits 每個(gè)read的最大匹配數(shù)，大于這個(gè)數(shù)的序列都丟棄掉。默認(rèn)參數(shù)是4.

-i/--min-intron 最小intron長度，默認(rèn)是50

-I/--max-intron 最大intron長度，默認(rèn)是300000

--non-canonical 同樣也搜索非經(jīng)典的junction，我們測(cè)試的結(jié)果是這個(gè)參數(shù)能夠提高junction數(shù)目，但是并不明顯。

-m/--splice-mis 允許第一以及最后一個(gè)部分的最大不匹配數(shù)目。允許范圍是0-2，默認(rèn)參數(shù)是1。

--max-append-mis 允許匹配高出錯(cuò)率片段以及鄰近的低出錯(cuò)片段不匹配的數(shù)目。默認(rèn)參數(shù)是3

--ins 最大插入長度，默認(rèn)是6，范圍是0-10

--del 最大刪除長度，默認(rèn)是6，范圍是0-10

--fusion| --fusion-non-canonical 查找融合（非經(jīng)典）基因

--filtering junctions過濾級(jí)別，取值為1,2.默認(rèn)是2；1代表更高的敏感度。2是標(biāo)準(zhǔn)過濾。