一 引言

1.1編寫目的

進行該測試以及撰寫此報告有以下幾個目的

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

2.深入分析各種軟件的各項參數。

3.深入分析各個軟件對應各種數據的情況。

1.2背景

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

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

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

1.3用戶群

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

其他讀者：項目及銷售相關人員。

1.4 數據對象：

1.5 測試階段

本次測試一共分兩個部分：

1. 各個軟件參數摸索；

2. 調整以后參數對后續(xù)結果的影響分析。

1.6測試工具

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’

二 測試概要

主要測試內容如下：

1. EBSeq參數測試

2. DEGSeq參數測試

3. DESeq參數測試

4. edgeR參數測試

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

2.1工作計劃進展

測試內容	計劃開始時間	實際開始時間	計劃完成時間	實際完成時間	工作完成情況
EBSeq參數測試	2013年10月25日	2013年10月25日	2013年10月27日	2013年10月27日	順利
DEGseq參數測試	2013年10月28日	2013年10月28日	2013年10月29日	2013年10月29日	順利
DESeq參數測試	2013年11月2日	2013年11月2日	2013年11月3日	2013年11月3日	順利
EdgeR參數測試	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日	順利
撰寫報告	2013年11月6日	2013年11月6日			順利

表1.工作進度表

2.2測試執(zhí)行

此次測試嚴格按照項目計劃和測試計劃執(zhí)行，按時完成了測試計劃規(guī)定的測試對象的測試。針對測試計劃制定規(guī)定的測試策略，依據測試計劃和測試用例，將網絡數據以及我們觀測的關鍵參數進行了完整的測試。

2.3測試用例

2.3.1功能性

1.測試如何提高差異基因篩選的方法與參數

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

三 測試環(huán)境

3.1軟硬件環(huán)境

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

四 測試結果

4.1 使用默認參數進行測試

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

使用包	篩選結果
DEGSeq	90
DESeq	12
EBSeq	113
EdgeR	8

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

表3 默認參數篩選出的差異基因數目

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

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

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

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

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

4.3 DESeq默認參數與優(yōu)化參數比較結果

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

使用默認參數

下調GO的結果，一共有103個。前30個結果截圖如下

圖1 下調GO篩選P-value結果圖

上調GO的結果，一共有48個。前30個截圖如下：

圖2 上調GO篩選P-value結果圖

按照FDR值小于0.05進行篩選，得出結果如下:

使用默認參數

下調GO的結果，一共20個，截圖如下

圖3 下調GO篩選fdr結果圖

上調GO的結果，

沒有找到結果；前30個結果的截圖如下：

圖4 上調GO篩選fdr結果圖

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

下調GO的結果，一共有103個。前30個結果截圖如下

圖5 下調GO篩選P-value結果圖

上調GO的結果，一共有48個。前30個截圖如下

圖6 上調GO篩選P-value結果圖

根據fdr值小于0.05進行篩選

下調GO的結果，一共20個，截圖如下

圖7 下調GO篩選fdr結果圖

上調GO的結果，

沒有找到結果；前30個結果的截圖如下：

圖8 上調GO篩選fdr結果圖

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

接下來我們分析了Pathway的結果

使用默認參數：按照P-value小于0.05進行篩選的結果如下

下調代謝通路：一共篩選出19個結果，截圖如下

圖9 下調Pathway篩選P-value結果圖

上調代謝通路：一共篩選出50個結果，前30個截圖如下：

圖10 上調Pathway篩選P-value結果圖

使用默認參數：按照fdr小于0.05進行篩選，結果如下

下調代謝通路：一共篩選出3果，截圖如下

圖11 下調Pathway篩選fdr結果圖

上調代謝通路：一共篩選出3結果，截圖如下：

圖12上調Pathway篩選fdr結果圖

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

按照P-value小于0.05進行篩選的結果如下

下調代謝通路：一共篩選出19個結果，截圖如下

圖13下調Pathway篩選P-value結果圖

上調代謝通路：一共篩選出50個結果，前30個截圖如下：

圖14上調Pathway篩選P-value結果圖

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

下調代謝通路：一共篩選出3果，截圖如下

圖15下調Pathway篩選fdr結果圖

上調代謝通路：一共篩選出3結果，截圖如下

圖16上調Pathway篩選fdr結果圖

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

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

4.4 EBSeq默認參數與優(yōu)化參數比較結果

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

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

之后我們做了GO分析，分析結果如下

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

使用默認參數

篩選下調GO結果，一共篩選出111個數據。前30個數據截圖如下：

圖16下調GO篩選P-value結果圖

篩選上調結果，一共篩選出47個數據。前30個數據截圖如下

圖17上調GO篩選P-value結果圖

篩選fdr小于0.05的結果如下

篩選下調GO結果，一共得到23個結果，截圖如下

圖18下調GO篩選fdr結果圖

篩選上調GO結果，沒有篩選到結果，截取前30個結果如下

圖19上調GO篩選fdr結果圖

對于優(yōu)化以后的參數

我們篩選p-value小于0.05

篩選下調GO結果，一共得到111個結果，前30個結果如下

圖20下調GO篩選P-value結果圖

篩選上調GO結果，一共得到47個結果，前30個結果如下

圖21上調GO篩選P-value結果圖

接下來我們分析了Pathway的結果

使用默認參數

按照P-value小于0.05進行篩選的結果如下

下調代謝通路，一共篩選出19個結果，截圖如下

圖22下調GO篩選P-value結果圖

上調代謝通路，一共篩選出43個結果，前30個結果截圖如下

圖23上調GO篩選P-value結果圖

我們使用fdr小于0.05進行篩選

下調代謝通路，一共篩選出3個結果，截圖如下

圖24下調GO篩選fdr結果圖

對于上調通路，一共篩選出2個結果，截圖如下：

圖25上調GO篩選fdr結果圖

對于優(yōu)化以后的參數

我們篩選p-value小于0.05

篩選下調GO結果，一共得到19個結果，截圖如下

圖26下調GO篩選P-value結果圖

篩選上調GO結果，一共得到43個結果，前30個結果截圖如下

圖27上調GO篩選fdr結果圖

篩選fdr小于0.05的結果如下

篩選下調GO結果，一共得到3個結果，截圖如下

圖28下GO篩選fdr結果圖

篩選上調GO的結果，一共得到2個結果，截圖如下

圖29上調GO篩選fdr結果圖

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

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

4.5 EdgeR等其他軟件的測試

正如上述所說，我們在測試edgeR以及DEGSeq的過程中，發(fā)現了一些能夠改進的方法， 但是由于科學必須遵守嚴謹的作風。所以我們并沒有盲目的采用這些做法。

五．測試結論與討論

5.1 DESeq參數改進及解釋

對DESeq的方法而言，程序默認的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  )

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

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

5.2 DESeq腳本中其余幾個函數的參數嘗試

總結DESeq中的幾個參數，在estimateDispersions函數中還有一個我們測試過的參數，fitType，這個參數的作用是設置用于計算dispersion-mean relation關系的方法，一共有兩個參數，parametric，以及local，前者是使用公式dispersion 與asymptDisp， extraPois之間的公式進行計算，而local則是指使用locfit包來進行計算。我們在實驗的過程中，發(fā)現這兩個值的選取對于結果并沒有影響，所以我們并沒有把這個參數的計算放入結果的討論中。

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

nbinomTest，是整個函數中計算差異結果的核心。我們查閱了函數手冊，這個函數并沒有能夠改變靈敏度的參數。我們在實驗的過程中也沒找到對應的參數，所以并沒有進行優(yōu)化。

5.3 EBSeq參數改進及解釋

程序默認的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)化以后的參數如下：

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

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

縱觀整個EBSeq的腳本，函數手冊及其用戶幫助，我們并沒有找到與DESeq，edgeR中都具有的參數dispersion，我們認為這可能是導致EBSeq結果提升不如DESeq明顯的主要因素。

5.4EBSeq腳本中其余幾個函數的參數嘗試

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

5.5edgeR參數改進及解釋

雖然在實驗的過程中，我們并沒有很合理的改善edgeR篩選差異基因的靈敏度的能力，但是我們還是對這個軟件進行了探索。

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中，對于單組數據的程序如下：

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

在我們上述貼出來的代碼中，一共分為3個部分，分別是計算多個樣本重復的兩種方法（存在檢驗（目前平臺使用的方法）以及線性擬合實驗（LRT））的方法，在我們的測試中，這兩種方法對于結果沒有明顯影響。所以我們認為使用哪一種方法都是可行的。在用戶手冊中，作者建議在多個因素（maulti factor）的情況下使用LRT的方法。

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

5.6 DESeq以及DESeq2

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

1.使用SummarizedExperiment用于儲存輸入文件，即時計算結果，以及結果。

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.所有的估計都是基于通用的線性模型（generalized linear model），包括在兩個條件的情況（之前的版本中使用的存在檢驗方法）

5.加入了Wald檢驗方法，同時之前的釋然率檢驗的方法也同樣能夠使用。

6.允許輸入標準化因子（normalized factors）

5.7 各個軟件的輸入文件

軟件	輸入數據	是否能夠輸入標準化數據
DESeq	Counts	否
DESeq2	Counts/	否
EBSeq	Counts	否
DEGSeq	RPKM	是
edgeR	Counts	否

表5 各個軟件輸入文本以及是否能夠輸入標準化數據

以上信息我們是通過查閱相關用戶手冊獲得，并未在實驗中進行驗證。

5.8 DESeq涉及多個條件下的method選擇

在實際生產中，我們發(fā)現，當涉及多個條件的情況是，DESeq中選擇per-condition的方法并不能進行計算，因此我們必須選擇這種情況下能夠使用的方法。因此我們使用了實際生產的數據進行分析得出的結果如下：

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

六．測試總結

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

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

3.由于edgeR中沒有能夠提供一個比較科學的獲取dispersion的方法，所以我們認為并沒有很好的改變結果。但是值得注意的一點是，在技術手冊中提到，當重復是屬于技術重復的時候，是可以使用bcv=0.01這個值的。