一 引言

1.1編寫目的

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

1.目前我們公司使用的幾種經(jīng)典算法，對(duì)于差異基因的篩選過于嚴(yán)格，導(dǎo)致結(jié)果后續(xù)分析的過程中遇到數(shù)據(jù)過少的情況。因此我們希望能夠找到更好的篩選出更多可能性，卻又不影響后續(xù)分析的參數(shù)來(lái)改善這種情況。

2.深入分析各種軟件的各項(xiàng)參數(shù)。

3.深入分析各個(gè)軟件對(duì)應(yīng)各種數(shù)據(jù)的情況。

1.2背景

差異基因（diffGene）篩選的定義是用統(tǒng)計(jì)學(xué)的方法對(duì)高通量測(cè)序的結(jié)果進(jìn)行篩選，跳出樣本間有顯著差異的基因。通過分析實(shí)驗(yàn)組，控制組；不同時(shí)間下的各個(gè)組之間的差異基因，可以幫助科學(xué)家找到與現(xiàn)象相關(guān)的基因，用于后續(xù)的實(shí)驗(yàn)驗(yàn)證。通過這種方法，大大加速了實(shí)驗(yàn)研究的進(jìn)度。能夠更加高效的進(jìn)行生物學(xué)研究。

目前在差異基因篩選中，世界上最常用的軟件有：DESeq，EBSeq，DEGSeq，以及其他的一些方法，比如TTest。目前的這些軟件中，通用的優(yōu)化方法有FPKM（RKPM），UpQuantile，median，TMM，TPM。在上述提到的這些軟件中，分別使用到了這些方法中的一種或者幾種方法，我們?cè)谟懻撝袝?huì)詳細(xì)的解釋這些結(jié)果。

鑒于我們目前對(duì)于差異基因的篩選結(jié)果過于嚴(yán)格，這樣使得我們后續(xù)的分析不能很好的展開，所以開展這個(gè)分析項(xiàng)目，以改善目前的這種情況。

1.3用戶群

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

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

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

1.5 測(cè)試階段

本次測(cè)試一共分兩個(gè)部分：

1. 各個(gè)軟件參數(shù)摸索；

2. 調(diào)整以后參數(shù)對(duì)后續(xù)結(jié)果的影響分析。

1.6測(cè)試工具

R Version 3.0.0

EdgeR version 3.4.0

DEGSeq version 1.2.2

DESeq version 1.14.0

DESeq2 version 1.2.5

EBSeq version 1.3.1

1.7 參考資料

edgeR: dierential expression analysisof digital gene expression data User's Guide

Package ‘edgeR’
Dierential expression of RNA-Seq data at the gene level --the DESeq package

Package ‘DESeq’

How to use the DEGseq Package

DEGseq

EBSeq: An R package for differential expression analysis using RNA-seq data

Package ‘EBSeq’

二 測(cè)試概要

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

1. EBSeq參數(shù)測(cè)試

2. DEGSeq參數(shù)測(cè)試

3. DESeq參數(shù)測(cè)試

4. edgeR參數(shù)測(cè)試

5. 后續(xù)GO_Pathway分析驗(yàn)證

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

測(cè)試內(nèi)容	計(jì)劃開始時(shí)間	實(shí)際開始時(shí)間	計(jì)劃完成時(shí)間	實(shí)際完成時(shí)間	工作完成情況
EBSeq參數(shù)測(cè)試	2013年10月25日	2013年10月25日	2013年10月27日	2013年10月27日	順利
DEGseq參數(shù)測(cè)試	2013年10月28日	2013年10月28日	2013年10月29日	2013年10月29日	順利
DESeq參數(shù)測(cè)試	2013年11月2日	2013年11月2日	2013年11月3日	2013年11月3日	順利
EdgeR參數(shù)測(cè)試	2013年11月4日	2013年11月4日	2013年11月5日	2013年11月5日	順利
GO_Pathway分析	2013年11月5日	2013年11月5日	2013年11月6日	2013年11月6日	順利
撰寫報(bào)告	2013年11月6日	2013年11月6日			順利

表1.工作進(jìn)度表

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ù)

2.測(cè)試提升靈敏度以后的參數(shù)對(duì)于后續(xù)分析的影響。

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

3.1軟硬件環(huán)境

表2 軟硬件環(huán)境

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

4.1 使用默認(rèn)參數(shù)進(jìn)行測(cè)試

在這次的測(cè)試中，我們使用的默認(rèn)腳本都是我們公司自己平臺(tái)的腳本，并且也參閱了各個(gè)軟件對(duì)應(yīng)的腳本，確認(rèn)了我們使用的腳本是默認(rèn)腳本。在選擇數(shù)據(jù)方面，我們選擇了客戶潘雪潔的關(guān)于人類的數(shù)據(jù)，這個(gè)數(shù)據(jù)包含了兩個(gè)樣本，每個(gè)樣本具有3次重復(fù)。在使用各個(gè)軟件進(jìn)行處理以后，我們將得出的結(jié)果按照總共6個(gè)樣本總和小于18，小于60，小于600的數(shù)據(jù)全部去除用以后續(xù)分析，進(jìn)行這一步的主要目的是為了降低那些低豐度的表達(dá)量對(duì)于我們結(jié)果的影響。

使用包	篩選結(jié)果
DEGSeq	90
DESeq	12
EBSeq	113
EdgeR	8

*DEGSeq使用FPKM，其余均使用counts

表3 默認(rèn)參數(shù)篩選出的差異基因數(shù)目

4.2參數(shù)優(yōu)化探索

針對(duì)不同的軟件，我們參照其用戶手冊(cè)以及函數(shù)說明進(jìn)行了優(yōu)化，在優(yōu)化的過程中，我們著重考量了不同程序的預(yù)處理過程以及處理數(shù)據(jù)的核心部分的參數(shù)控制。針對(duì)我們優(yōu)化以后的結(jié)果，我們將在討論部分進(jìn)行詳細(xì)介紹。

*DEGSeq使用FPKM，其余均使用counts

表4 優(yōu)化以后篩選出的差異基因數(shù)目

在上述結(jié)果中，我們看到只有DESeq以及EBSeq的結(jié)果有了顯著提升，而其他兩個(gè)軟件的結(jié)果并沒有變化。對(duì)于DEGseq，由于使用的是FPKM的值，我們?cè)诓檎覅?shù)變化的時(shí)候并沒有找到很好的參數(shù)，所以并沒有產(chǎn)生變化。對(duì)于EdgeR，我們?cè)跍y(cè)試的過程中，只發(fā)現(xiàn)了調(diào)節(jié)dipersion這個(gè)參數(shù)能夠很好的改變結(jié)果，但是在查詢這個(gè)參數(shù)的統(tǒng)計(jì)學(xué)意義的過程中，我們發(fā)現(xiàn)這個(gè)參數(shù)是用于反應(yīng)樣本之間的數(shù)值隨機(jī)變化的顯著性，即當(dāng)這個(gè)參數(shù)為設(shè)置為0的時(shí)候就會(huì)認(rèn)為數(shù)值之間的任何數(shù)值差異（即便相差1%）都會(huì)認(rèn)為這兩者在這次統(tǒng)計(jì)上是有差異的，于是我們參照了函數(shù)聲明中求這個(gè)參數(shù)的函數(shù)進(jìn)行估計(jì)，并沒有發(fā)現(xiàn)軟件支持的很好提高靈敏度的方法，同時(shí)我們認(rèn)為人為的定義一個(gè)值在科學(xué)的意義上而言顯得很不嚴(yán)謹(jǐn)。所以我們的放棄了這個(gè)軟件的提升靈敏度的嘗試。

4.3 DESeq默認(rèn)參數(shù)與優(yōu)化參數(shù)比較結(jié)果

在實(shí)驗(yàn)過程中，我們發(fā)現(xiàn)優(yōu)化前后，DESeq的結(jié)果變化選出的基因從12個(gè)提升到了697個(gè)，提升了685個(gè)。是所有變化中最多的，為了確定結(jié)果對(duì)后續(xù)分析沒有影響，我們使用輸出的結(jié)果進(jìn)行了后續(xù)的GO，Pathway分析。在處理之前，為了降低樣本低豐度結(jié)果導(dǎo)致的影響，我們將輸入結(jié)果按照六列樣本數(shù)據(jù)和小于18,60,600的量分為了三次重復(fù)。我們希望通過這個(gè)做法來(lái)知道低豐度對(duì)于結(jié)果的影響。以下是我們?cè)敿?xì)分析的刪除小于18的結(jié)果的GO，Pathway分析的結(jié)果。

使用默認(rèn)參數(shù)

下調(diào)GO的結(jié)果，一共有103個(gè)。前30個(gè)結(jié)果截圖如下

圖1 下調(diào)GO篩選P-value結(jié)果圖

上調(diào)GO的結(jié)果，一共有48個(gè)。前30個(gè)截圖如下：

圖2 上調(diào)GO篩選P-value結(jié)果圖

按照FDR值小于0.05進(jìn)行篩選，得出結(jié)果如下:

使用默認(rèn)參數(shù)

下調(diào)GO的結(jié)果，一共20個(gè)，截圖如下

圖3 下調(diào)GO篩選fdr結(jié)果圖

上調(diào)GO的結(jié)果，

沒有找到結(jié)果；前30個(gè)結(jié)果的截圖如下：

圖4 上調(diào)GO篩選fdr結(jié)果圖

之后我們使用我們優(yōu)化的結(jié)果進(jìn)行GO分析，其結(jié)果如下

下調(diào)GO的結(jié)果，一共有103個(gè)。前30個(gè)結(jié)果截圖如下

圖5 下調(diào)GO篩選P-value結(jié)果圖

上調(diào)GO的結(jié)果，一共有48個(gè)。前30個(gè)截圖如下

圖6 上調(diào)GO篩選P-value結(jié)果圖

根據(jù)fdr值小于0.05進(jìn)行篩選

下調(diào)GO的結(jié)果，一共20個(gè)，截圖如下

圖7 下調(diào)GO篩選fdr結(jié)果圖

上調(diào)GO的結(jié)果，

沒有找到結(jié)果；前30個(gè)結(jié)果的截圖如下：

圖8 上調(diào)GO篩選fdr結(jié)果圖

我們將上述的結(jié)果進(jìn)行了對(duì)照分析，發(fā)現(xiàn)兩者的結(jié)果是一一對(duì)應(yīng)的。這證明了我們提升靈敏度的方法對(duì)于GO分析沒有影響。

接下來(lái)我們分析了Pathway的結(jié)果

使用默認(rèn)參數(shù)：按照P-value小于0.05進(jìn)行篩選的結(jié)果如下

下調(diào)代謝通路：一共篩選出19個(gè)結(jié)果，截圖如下

圖9 下調(diào)Pathway篩選P-value結(jié)果圖

上調(diào)代謝通路：一共篩選出50個(gè)結(jié)果，前30個(gè)截圖如下：

圖10 上調(diào)Pathway篩選P-value結(jié)果圖

使用默認(rèn)參數(shù)：按照fdr小于0.05進(jìn)行篩選，結(jié)果如下

下調(diào)代謝通路：一共篩選出3果，截圖如下

圖11 下調(diào)Pathway篩選fdr結(jié)果圖

上調(diào)代謝通路：一共篩選出3結(jié)果，截圖如下：

圖12上調(diào)Pathway篩選fdr結(jié)果圖

使用優(yōu)化以后的參數(shù)

按照P-value小于0.05進(jìn)行篩選的結(jié)果如下

下調(diào)代謝通路：一共篩選出19個(gè)結(jié)果，截圖如下

圖13下調(diào)Pathway篩選P-value結(jié)果圖

上調(diào)代謝通路：一共篩選出50個(gè)結(jié)果，前30個(gè)截圖如下：

圖14上調(diào)Pathway篩選P-value結(jié)果圖

使用優(yōu)化過后的參數(shù)：按照fdr小于0.05進(jìn)行篩選，結(jié)果如下

下調(diào)代謝通路：一共篩選出3果，截圖如下

圖15下調(diào)Pathway篩選fdr結(jié)果圖

上調(diào)代謝通路：一共篩選出3結(jié)果，截圖如下

圖16上調(diào)Pathway篩選fdr結(jié)果圖

我們將上述的結(jié)果進(jìn)行了對(duì)照分析，發(fā)現(xiàn)兩者的結(jié)果是一一對(duì)應(yīng)的。我們發(fā)現(xiàn)，無(wú)論是否按照排序還是飛排序的情況，我們的各種GO/PathwayID都是一一對(duì)應(yīng)的。這證明了我們提升靈敏度的方法對(duì)于GO分析沒有影響的預(yù)期。

這個(gè)結(jié)果證明我們對(duì)于DESeq使用的優(yōu)化方法是可行的。

4.4 EBSeq默認(rèn)參數(shù)與優(yōu)化參數(shù)比較結(jié)果

在實(shí)驗(yàn)過程中，我們發(fā)現(xiàn)優(yōu)化前后，EBSeq的結(jié)果變化選出的基因從12個(gè)提升到了697個(gè)，提升了685個(gè)。是所有變化中最多的，為了確定結(jié)果對(duì)后續(xù)分析沒有影響，我們使用輸出的結(jié)果進(jìn)行了后續(xù)的GO，Pathway分析。與DESeq中的處理方式相似，在處理之前，為了降低樣本低豐度結(jié)果導(dǎo)致的影響，我們將輸入結(jié)果按照六列樣本數(shù)據(jù)和小于18,60,600的量分為了三次重復(fù)。我們希望通過這個(gè)做法來(lái)知道低豐度對(duì)于結(jié)果的影響。以下是我們?cè)敿?xì)分析的刪除小于18的結(jié)果的GO，Pathway分析的結(jié)果。

可以看出，在EBSeq中，篩選出的基因從113個(gè)提升到了317個(gè)，提升了204個(gè)。

之后我們做了GO分析，分析結(jié)果如下

篩選P-value小于0.05的值如下

使用默認(rèn)參數(shù)

篩選下調(diào)GO結(jié)果，一共篩選出111個(gè)數(shù)據(jù)。前30個(gè)數(shù)據(jù)截圖如下：

圖16下調(diào)GO篩選P-value結(jié)果圖

篩選上調(diào)結(jié)果，一共篩選出47個(gè)數(shù)據(jù)。前30個(gè)數(shù)據(jù)截圖如下

圖17上調(diào)GO篩選P-value結(jié)果圖

篩選fdr小于0.05的結(jié)果如下

篩選下調(diào)GO結(jié)果，一共得到23個(gè)結(jié)果，截圖如下

圖18下調(diào)GO篩選fdr結(jié)果圖

篩選上調(diào)GO結(jié)果，沒有篩選到結(jié)果，截取前30個(gè)結(jié)果如下

圖19上調(diào)GO篩選fdr結(jié)果圖

對(duì)于優(yōu)化以后的參數(shù)

我們篩選p-value小于0.05

篩選下調(diào)GO結(jié)果，一共得到111個(gè)結(jié)果，前30個(gè)結(jié)果如下

圖20下調(diào)GO篩選P-value結(jié)果圖

篩選上調(diào)GO結(jié)果，一共得到47個(gè)結(jié)果，前30個(gè)結(jié)果如下

圖21上調(diào)GO篩選P-value結(jié)果圖

接下來(lái)我們分析了Pathway的結(jié)果

使用默認(rèn)參數(shù)

按照P-value小于0.05進(jìn)行篩選的結(jié)果如下

下調(diào)代謝通路，一共篩選出19個(gè)結(jié)果，截圖如下

圖22下調(diào)GO篩選P-value結(jié)果圖

上調(diào)代謝通路，一共篩選出43個(gè)結(jié)果，前30個(gè)結(jié)果截圖如下

圖23上調(diào)GO篩選P-value結(jié)果圖

我們使用fdr小于0.05進(jìn)行篩選

下調(diào)代謝通路，一共篩選出3個(gè)結(jié)果，截圖如下

圖24下調(diào)GO篩選fdr結(jié)果圖

對(duì)于上調(diào)通路，一共篩選出2個(gè)結(jié)果，截圖如下：

圖25上調(diào)GO篩選fdr結(jié)果圖

對(duì)于優(yōu)化以后的參數(shù)

我們篩選p-value小于0.05

篩選下調(diào)GO結(jié)果，一共得到19個(gè)結(jié)果，截圖如下

圖26下調(diào)GO篩選P-value結(jié)果圖

篩選上調(diào)GO結(jié)果，一共得到43個(gè)結(jié)果，前30個(gè)結(jié)果截圖如下

圖27上調(diào)GO篩選fdr結(jié)果圖

篩選fdr小于0.05的結(jié)果如下

篩選下調(diào)GO結(jié)果，一共得到3個(gè)結(jié)果，截圖如下

圖28下GO篩選fdr結(jié)果圖

篩選上調(diào)GO的結(jié)果，一共得到2個(gè)結(jié)果，截圖如下

圖29上調(diào)GO篩選fdr結(jié)果圖

我們比較使用默認(rèn)參數(shù)與優(yōu)化以后的參數(shù)的變化，發(fā)現(xiàn)這兩者做出的GO，以及Pathway的結(jié)果是一致的。無(wú)論是否按照排序都是出現(xiàn)相同的結(jié)果。

所以我們可以知道，我們對(duì)于EBSeq的處理方法是可行的。

4.5 EdgeR等其他軟件的測(cè)試

正如上述所說，我們?cè)跍y(cè)試edgeR以及DEGSeq的過程中，發(fā)現(xiàn)了一些能夠改進(jìn)的方法， 但是由于科學(xué)必須遵守嚴(yán)謹(jǐn)?shù)淖黠L(fēng)。所以我們并沒有盲目的采用這些做法。

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

5.1 DESeq參數(shù)改進(jìn)及解釋

對(duì)DESeq的方法而言，程序默認(rèn)的R腳本如下：

library(DESeq2) data = read.table(fileName, he=T, sep="\t", row.names=1) conds = factor( c("c", "c", "p", "p", "p", "c") ) colData=data.frame(conds) cds =  DESeqDataSetFromMatrix( data, colData,,formula(~ conds) ) cds = estimateSizeFactors(cds) cds =  estimateDispersions(cds) res = nbinomTest( cds, "p", "c" ) write.table( res, file="/home/novelbio/hxl_temp/DifferenceExpression_result/CPVS_DESeq.xls",sep="\t",row.names=F  )

我們優(yōu)化的腳本如下

filePath = "/home/novelbio/hxl_temp" fileName = "/home/novelbio/hxl_temp/All_Counts600.txt" setwd(filePath) library(DESeq) data = read.table(fileName, he=T, sep="\t", row.names=1) conds = factor( c("c", "c", "p", "p", "p", "c") ) cds = newCountDataSet( data, conds ) cds = estimateSizeFactors(cds) cds =  estimateDispersions(cds,method="per-condition",sharingMode="gene-est-only") res = nbinomTest( cds, "p", "c" ) write.table( res, file="/home/novelbio/hxl_temp/CPVS_DESeq600_optimization.xls",sep="\t",row.names=F  )

在上述命令中，我們的改動(dòng)使用的紅色進(jìn)行標(biāo)注。其解釋如下：estimateDispersions這個(gè)函數(shù)的功能是用于估計(jì)dispersion的參數(shù)，正如我們上述所說，這個(gè)參數(shù)是用于衡量樣本之間數(shù)據(jù)變化對(duì)于很衡量樣本差異的靈敏度的值，當(dāng)這個(gè)值設(shè)為0的時(shí)候，就意味著即使兩個(gè)樣本有一點(diǎn)點(diǎn)的差異，也認(rèn)為這個(gè)數(shù)據(jù)是能夠反映基因表達(dá)的差異的，也就是說，這個(gè)數(shù)據(jù)越小，越能說明靈敏的反映基因差異的顯著性。在這里，我們額外的加入了兩個(gè)參數(shù)，分別是method，這個(gè)參數(shù)是指定計(jì)算dispersion的方法，而我們指定的per-condition是指，針對(duì)不同的樣本重復(fù)計(jì)算出不等的dispersion值，在我們的例子中，我們一個(gè)樣本有三組重復(fù)，就意味著每一次重復(fù)都對(duì)應(yīng)了一個(gè)離差值。而通過查閱函數(shù)手冊(cè)，我們知道計(jì)算離差的方法一共有4種，分別是pooled，pooled-CR，per-condition，blind。其中，第一種pooled的方法最為粗糙， 這個(gè)方法對(duì)于所有的數(shù)據(jù)而言，只估計(jì)了一個(gè)經(jīng)驗(yàn)性的dispersion值，并且將這個(gè)值用于所有的樣本，顯然這個(gè)做法并不嚴(yán)謹(jǐn)，在使用的過程中我們并沒有使用這種做法。第二種pooled-CR的做法是通過Cox-Reid adjusted pro?le likelihood的方法估計(jì)出dispersion值，顯然這種做法的運(yùn)算速度要慢于第一種做法。第三種方法就是我們采用的方法，，對(duì)于不同條件的數(shù)據(jù)分別計(jì)算不同的dispersion值。但是這種方法在用于多種條件分析的情況（并非樣本重復(fù)）時(shí)，并不適用。即當(dāng)我們?cè)谕淮芜\(yùn)行中，不應(yīng)該使用多個(gè)比較（比如同時(shí)分析時(shí)間點(diǎn)以及實(shí)驗(yàn)，對(duì)照組，但是分開是可以使用的）。第四種方法blind的做法是認(rèn)為所有樣本都是同一條件下的不同重復(fù)，忽略樣本的標(biāo)記（labels），并計(jì)算經(jīng)驗(yàn)性的dispersion值。即便是沒有生物學(xué)重復(fù)的情況下也能使用，但是這種方法會(huì)導(dǎo)致樣本失去某些重要的東西（loss of power），顯然這種方法也是不嚴(yán)謹(jǐn)?shù)?。綜合考慮這些方法的原理以及我們?cè)趯?shí)驗(yàn)中的結(jié)果，我們認(rèn)為使用per-condition的方法顯然是更加合理的。我們加入的另一個(gè)參數(shù)sharing-mode的作用是選擇將優(yōu)化以后的數(shù)據(jù)還是沒有優(yōu)化的數(shù)據(jù)進(jìn)行后續(xù)的分析，在這里，一共有三個(gè)參數(shù)，fit-only，maximum，gene-est-only。這三個(gè)參數(shù)的解釋如下：fit-only的作用是只使用優(yōu)化以后的數(shù)據(jù)。函數(shù)手冊(cè)中建議只有當(dāng)具有少量的重復(fù)的時(shí)候才能使用這種方法。Maximum的方法是在優(yōu)化的數(shù)據(jù)以及未優(yōu)化的數(shù)據(jù)中選取較大值。函數(shù)手冊(cè)建議當(dāng)只有3個(gè)重復(fù)以下的時(shí)候使用這個(gè)參數(shù)。而第三個(gè)參數(shù)，gene-est-only在樣本重復(fù)比較多并且dispersion值是可靠的時(shí)候適用，當(dāng)樣本重復(fù)少的時(shí)候，這個(gè)參數(shù)會(huì)比較靈敏的反映出樣本變化情況。我們?cè)跍y(cè)試的過程中發(fā)現(xiàn)，當(dāng)sharing-mode選擇gene-est-only的時(shí)候，得出的結(jié)果能夠很好的提高靈敏度，這與我們進(jìn)行試驗(yàn)的目的是一致的。所以我們選擇了這個(gè)參數(shù)。而method選擇的目的則是為了更好的提高我們實(shí)驗(yàn)過程中的嚴(yán)謹(jǐn)性。

通過后續(xù)的GO，Pathway分析，我們發(fā)現(xiàn)我們優(yōu)化以后以及優(yōu)化以前的的數(shù)據(jù)做出的GO，Pathway分析的結(jié)果是一致的，這直接的證明了我們處理方式是合理的。我們可以看出，無(wú)論是使用P-value還是FDR進(jìn)行篩選，其結(jié)果都是一致的。由于GO,Pathway分析后續(xù)的步驟都是相同，即，使用相同的數(shù)據(jù)其結(jié)果也是相同的，所以我們并沒有做畫圖分析。

5.2 DESeq腳本中其余幾個(gè)函數(shù)的參數(shù)嘗試

總結(jié)DESeq中的幾個(gè)參數(shù)，在estimateDispersions函數(shù)中還有一個(gè)我們測(cè)試過的參數(shù)，fitType，這個(gè)參數(shù)的作用是設(shè)置用于計(jì)算dispersion-mean relation關(guān)系的方法，一共有兩個(gè)參數(shù)，parametric，以及local，前者是使用公式dispersion 與asymptDisp， extraPois之間的公式進(jìn)行計(jì)算，而local則是指使用locfit包來(lái)進(jìn)行計(jì)算。我們?cè)趯?shí)驗(yàn)的過程中，發(fā)現(xiàn)這兩個(gè)值的選取對(duì)于結(jié)果并沒有影響，所以我們并沒有把這個(gè)參數(shù)的計(jì)算放入結(jié)果的討論中。

estimateSizeFactors，在DESeq設(shè)置sizeFactor的方法中，默認(rèn)使用median 的方法。官方手冊(cè)上所說在低count數(shù)的情況下使用shorth的方法我們并沒有在R語(yǔ)言中找到。指定方法所使用的參數(shù)是locfunc。

nbinomTest，是整個(gè)函數(shù)中計(jì)算差異結(jié)果的核心。我們查閱了函數(shù)手冊(cè)，這個(gè)函數(shù)并沒有能夠改變靈敏度的參數(shù)。我們?cè)趯?shí)驗(yàn)的過程中也沒找到對(duì)應(yīng)的參數(shù)，所以并沒有進(jìn)行優(yōu)化。

5.3 EBSeq參數(shù)改進(jìn)及解釋

程序默認(rèn)的R腳本如下：

filePath = "/home/novelbio/soft/NBCnew6/rscript/tmp/" fileName = "/home/novelbio/hxl_temp/All_Counts18.txt" setwd(filePath) library(EBSeq) data = read.table(fileName, he=T, sep="\t", row.names=1) data = as.matrix(data) size = QuantileNorm(data, 0.75) dataRun=data[,c(1, 2, 3, 4, 5, 6)] sizeRun = size[c(1, 2, 3, 4, 5, 6)] EBOutRun<- EBTest(Data=dataRun, Conditions=as.factor(c("c", "c", "p", "p", "p", "c")),sizeFactors=sizeRun, maxround=5) GenePP1=GetPPMat(EBOutRun) FDR<-GenePP1[,1] c<- unlist(EBOutRun$C1Mean) p<- unlist(EBOutRun$C2Mean) LogFC=log2(PostFC(EBOutRun)$RealFC) out<-cbind(c,p, LogFC,FDR ) colnames(out)<-c("c","p","LogFC", "FDR") write.table(out, file="/home/novelbio/hxl_temp/CPVS_EBSeq18_withOutOptimization.xls", row.names=TRUE,col.names=TRUE,sep="\t")

我們優(yōu)化以后的參數(shù)如下：

fileName = "/home/novelbio/hxl_temp/All_Counts600.txt" library(EBSeq) data = read.table(fileName, he=T, sep="\t", row.names=1) data = as.matrix(data) size = QuantileNorm(data,0.75) dataRun=data[,c(1,2,3,4,5,6)] sizeRun = size[c(1,2,3,4,5,6)] EBOutRun<- EBTest(Data=dataRun, Conditions=as.factor(c("C", "C", "P", "P", "P", "C")) ,sizeFactors=sizeRun, maxround=1 ,ApproxVal=0.3) GenePP1=GetPPMat(EBOutRun) FDR<-GenePP1[,1]

C<- unlist(EBOutRun$C1Mean) P<- unlist(EBOutRun$C2Mean) LogFC=log2(PostFC(EBOutRun)$RealFC) out<-cbind(C,P, LogFC,FDR ) colnames(out)<-c("DF","HF","LogFC", "FDR") write.table(out, file="/home/novelbio/hxl_temp/CVSP_EBSeq600_optimization.xls", row.names=TRUE,col.names=TRUE,sep="\t")

在上述命令中，我們的改動(dòng)使用的紅色進(jìn)行標(biāo)注。其解釋如下：maxround是指EBTest迭代的次數(shù)，迭代次數(shù)越多，得出的結(jié)果越少。程序默認(rèn)的迭代次數(shù)是5次。而我們公司平臺(tái)指定的次數(shù)是15次，這樣雖然看似很嚴(yán)謹(jǐn)，但是也會(huì)過濾掉很多有用的信息，顯然，這不是將實(shí)驗(yàn)結(jié)果最大化利用的最合理的方法。在實(shí)驗(yàn)的過程中，我們發(fā)現(xiàn)即便這個(gè)參數(shù)設(shè)置成為1，也不會(huì)對(duì)后續(xù)的分析產(chǎn)生影響。所以我們?cè)谶@里就大膽的設(shè)置成為1次迭代。ApproxVal的英文解釋是The variances of the transcripts with mean < var will be approximated as mean/(1-ApproxVal)，是用于調(diào)整mean的值的一個(gè)參數(shù)，我們?cè)趯?shí)驗(yàn)的過程中發(fā)現(xiàn)，當(dāng)這個(gè)值設(shè)置為0.3的時(shí)候，是最能提升差異基因篩選的結(jié)果的。而在后續(xù)的實(shí)驗(yàn)中，我們發(fā)現(xiàn)這個(gè)結(jié)果的變化并不會(huì)對(duì)結(jié)果產(chǎn)生影響，所以我們也接受了這個(gè)值的選擇。

縱觀整個(gè)EBSeq的腳本，函數(shù)手冊(cè)及其用戶幫助，我們并沒有找到與DESeq，edgeR中都具有的參數(shù)dispersion，我們認(rèn)為這可能是導(dǎo)致EBSeq結(jié)果提升不如DESeq明顯的主要因素。

5.4EBSeq腳本中其余幾個(gè)函數(shù)的參數(shù)嘗試

在EBSeq中，運(yùn)行的代碼可以大致分為兩個(gè)部分，數(shù)據(jù)預(yù)處理以及差異基因篩選，在數(shù)據(jù)預(yù)處理中，我們發(fā)現(xiàn)計(jì)算size這個(gè)變量的函數(shù)一共有如下幾個(gè)：quantileNorm，這個(gè)函數(shù)是使用百分位點(diǎn)的方法進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化；RankNorm，這個(gè)函數(shù)是利用rank標(biāo)準(zhǔn)化方法進(jìn)行標(biāo)準(zhǔn)化。另外一種標(biāo)準(zhǔn)化方法是MedianNorm，這個(gè)方法使用Anders et. al., 2010.提出的中位數(shù)標(biāo)準(zhǔn)化方法進(jìn)行標(biāo)準(zhǔn)化。在我們測(cè)試的過程中，針對(duì)這一組數(shù)據(jù)，我們發(fā)現(xiàn)不同的標(biāo)準(zhǔn)化方法得出的結(jié)果是幾乎一樣的。針對(duì)具體方法的具體使用對(duì)象，我們并沒有找到各自針對(duì)的情況。在函數(shù)EBTest中，還有幾個(gè)參數(shù)，比如pool，函數(shù)手冊(cè)建議在沒有樣本重復(fù)的時(shí)候設(shè)置這個(gè)參數(shù)為TRUE，而PoolLower，PoolUpper兩個(gè)參數(shù)這是指定只用于分析的數(shù)據(jù)集，即只取在這兩者范圍內(nèi)的數(shù)據(jù)集用于分析。其余的參數(shù)我們并沒有發(fā)現(xiàn)能夠改善分析步奏的嚴(yán)謹(jǐn)性或者提高差異基因篩選靈敏度能力的參數(shù)，所以并沒有進(jìn)行指定。對(duì)于這些參數(shù)，我們認(rèn)為使用軟件默認(rèn)的參數(shù)就不錯(cuò)了。

5.5edgeR參數(shù)改進(jìn)及解釋

雖然在實(shí)驗(yàn)的過程中，我們并沒有很合理的改善edgeR篩選差異基因的靈敏度的能力，但是我們還是對(duì)這個(gè)軟件進(jìn)行了探索。

filePath = "/home/novelbio/hxl_temp" fileName = "/home/novelbio/hxl_temp/diff.txt" setwd(filePath) data2 = read.table(fileName, he=T, sep="\t", row.names=1) library(edgeR) bcv = 0 plus = 0 data = read.table(fileName, he=T, sep="\t", row.names=1) data = data + plus group = factor( c("c", "c","p", "p", "p", "c" ) ) #######################一般方法##################### Normfactors=calcNormFactors(data,method="RLE") p = DGEList(counts=data,group=group) p = estimateCommonDisp(p) p= estimateTagwiseDisp(p) result = exactTest(p, pair=c("c","p" ) ,dispersion=p$tagwise.dispersion); result = exactTest(p, pair=c("c","p" ) ); resultFinal=result$table; result2=cbind(resultFinal,fdr=p.adjust(resultFinal[,4])); write.table(result2, file="/home/novelbio/hxl_temp/CPVS_EdgeR.xls",sep="\t") ###############LRT_Method################# design=model.matrix(~group) m = DGEList(counts=data,group=group) m=estimateGLMCommonDisp(m,design) m <- estimateGLMTrendedDisp(m,design) m<- estimateGLMTagwiseDisp(m,design) fit <- glmFit(m,design) lrt <- glmLRT(fit,coef=2) lrtFinal=lrt$table; lrt2=cbind(lrtFinal,fdr=p.adjust(lrtFinal[,4])); write.table(lrt2, file="/home/novelbio/hxl_temp/CPVS_EdgeR.xls",sep="\t")

在EdgeR中，對(duì)于單組數(shù)據(jù)的程序如下：

filePath = "/home/novelbio/hxl_temp" fileName = "/home/novelbio/hxl_temp/diff.txt" setwd(filePath) data2 = read.table(fileName, he=T, sep="\t", row.names=1) library(edgeR) bcv = 0 plus = 0 data = read.table(fileName, he=T, sep="\t", row.names=1) data = data + plus group = factor( c("c", "c","p", "p", "p", "c" ) ) #######################一般方法##################### Normfactors=calcNormFactors(data,method="RLE") p = DGEList(counts=data,group=group) p = estimateCommonDisp(p) result = exactTest(p, pair=c("c","p" ) ,dispersion=p$common.dispersion); result = exactTest(p, pair=c("c","p" ) ); resultFinal=result$table; result2=cbind(resultFinal,fdr=p.adjust(resultFinal[,4])); write.table(result2, file="/home/novelbio/hxl_temp/CPVS_EdgeR.xls",sep="\t")

我們通過查詢edgeR的函數(shù)手冊(cè)，發(fā)現(xiàn)這個(gè)結(jié)果中也有關(guān)于dispersion的函數(shù)，通過查閱函數(shù)手冊(cè)，發(fā)現(xiàn)有幾種可以指定這個(gè)值的方法，最簡(jiǎn)單的也是我們現(xiàn)在使用的，是先指定一個(gè)bcv（這個(gè)值是dispersion的平方根），我們平臺(tái)目前指定的是0.3，然后在exactTest函數(shù)中指定dispersion=bcv^2。在函數(shù)手冊(cè)中，有幾個(gè)關(guān)于這個(gè)值的經(jīng)驗(yàn)數(shù)值，比如對(duì)于人類數(shù)據(jù)這個(gè)值去0.4，對(duì)于相似的模式生物（identical model organisms）可以取0.1，對(duì)于重復(fù)樣本可以取0.01（technical replicates）.在我們看來(lái)，使用這種方法并不利于批量化作業(yè)，而且從科學(xué)的角度而言也顯得不嚴(yán)謹(jǐn)。其余估計(jì)dispersion的函數(shù)有estimateTagwiseDisp()這個(gè)函數(shù)與DESeq中的per-condition是類似的，都是按照列來(lái)估計(jì)dispersion。而estimateCommonDisp()這個(gè)函數(shù)則和DESeq中的pooled的方法類似，也是對(duì)于所有數(shù)據(jù)估計(jì)一個(gè)dispersion值，我們認(rèn)為這樣的做法并不嚴(yán)謹(jǐn)，因此也就沒有采用。

在我們上述貼出來(lái)的代碼中，一共分為3個(gè)部分，分別是計(jì)算多個(gè)樣本重復(fù)的兩種方法（存在檢驗(yàn)（目前平臺(tái)使用的方法）以及線性擬合實(shí)驗(yàn)（LRT））的方法，在我們的測(cè)試中，這兩種方法對(duì)于結(jié)果沒有明顯影響。所以我們認(rèn)為使用哪一種方法都是可行的。在用戶手冊(cè)中，作者建議在多個(gè)因素（maulti factor）的情況下使用LRT的方法。

在我們的測(cè)試過程中，我們發(fā)現(xiàn)，當(dāng)bcv設(shè)置為0的時(shí)候，結(jié)果能夠提升到1200個(gè)左右，因此，我們認(rèn)為在滿足條件（technical replicates）的時(shí)候，設(shè)置bcv=0.01的條件是可以提升這個(gè)結(jié)果的。

5.6 DESeq以及DESeq2

在測(cè)試DESeq的過程中，我們發(fā)現(xiàn)這個(gè)軟件推出了新版本，DESeq2。查閱相關(guān)手冊(cè)，我們發(fā)現(xiàn)DESeq2相比DESeq最大的改變是能夠直接讀取序列文件，然后自己計(jì)算counts數(shù)。其余并沒有發(fā)現(xiàn)變化。查閱DESeq2的用戶手冊(cè)，我們發(fā)現(xiàn)如下變化：

1.使用SummarizedExperiment用于儲(chǔ)存輸入文件，即時(shí)計(jì)算結(jié)果，以及結(jié)果。

2. Maximum a posteriori estimation of GLM coecients incorporating a zero-mean normal prior withvariance estimated from data (equivalent to Tikhonov/ridge regularization). This adjustment haslittle eect on genes with high counts, yet it helps to moderate the otherwise large spread in log2fold changes for genes with low counts (e. g. single digits per condition). 3. Maximum a posteriori estimation of dispersion replaces the sharingMode options fit-only or maximum of the previous version of the package.

4.所有的估計(jì)都是基于通用的線性模型（generalized linear model），包括在兩個(gè)條件的情況（之前的版本中使用的存在檢驗(yàn)方法）

5.加入了Wald檢驗(yàn)方法，同時(shí)之前的釋然率檢驗(yàn)的方法也同樣能夠使用。

6.允許輸入標(biāo)準(zhǔn)化因子（normalized factors）

5.7 各個(gè)軟件的輸入文件

軟件	輸入數(shù)據(jù)	是否能夠輸入標(biāo)準(zhǔn)化數(shù)據(jù)
DESeq	Counts	否
DESeq2	Counts/	否
EBSeq	Counts	否
DEGSeq	RPKM	是
edgeR	Counts	否

表5 各個(gè)軟件輸入文本以及是否能夠輸入標(biāo)準(zhǔn)化數(shù)據(jù)

以上信息我們是通過查閱相關(guān)用戶手冊(cè)獲得，并未在實(shí)驗(yàn)中進(jìn)行驗(yàn)證。

5.8 DESeq涉及多個(gè)條件下的method選擇

在實(shí)際生產(chǎn)中，我們發(fā)現(xiàn)，當(dāng)涉及多個(gè)條件的情況是，DESeq中選擇per-condition的方法并不能進(jìn)行計(jì)算，因此我們必須選擇這種情況下能夠使用的方法。因此我們使用了實(shí)際生產(chǎn)的數(shù)據(jù)進(jìn)行分析得出的結(jié)果如下：

從上述結(jié)果可以看出，效果最好的是pooled-CR，而Blind的method對(duì)應(yīng)也報(bào)錯(cuò)（與per-condition的一樣），而pooled的結(jié)果和pooled-CR的結(jié)果都比默認(rèn)參數(shù)的效果好，因此我們可以使用這兩種方法，之后我們比較了對(duì)應(yīng)兩種條件下pooled和per-condition的效果，發(fā)現(xiàn)per-condition的效果（11903）比pooled-CR（8848）的效果要好。因此對(duì)應(yīng)兩種條件下的選擇不變。

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

1.在測(cè)試過程中，我們成功在不影響后續(xù)結(jié)果的情況下提升了差異基因靈敏度的方法，并且證明了我們方法的正確性

2.對(duì)于edgeR以及DESEq而言，能夠很好的控制其靈敏度的參數(shù)是控制dispersion的值。對(duì)于EBSeq而言，則是控制迭代的次數(shù)以及控制ApproxVal的值。

3.由于edgeR中沒有能夠提供一個(gè)比較科學(xué)的獲取dispersion的方法，所以我們認(rèn)為并沒有很好的改變結(jié)果。但是值得注意的一點(diǎn)是，在技術(shù)手冊(cè)中提到，當(dāng)重復(fù)是屬于技術(shù)重復(fù)的時(shí)候，是可以使用bcv=0.01這個(gè)值的。