R 中的資料探勘演算法/序列挖掘/DEGSeq
使用 DEGseq
首先，我們必須選擇組織資料的方式。有兩種方法可以做到這一點：使用基於 MA 圖的隨機抽樣模型方法或使用基於 MA 圖的技術重複方法。
基於 MA 圖的隨機抽樣模型方法
MA 圖是一種統計分析工具，已用於檢測和視覺化微陣列資料的強度依賴比率。設 C1 和 C2 表示從兩個樣本獲得的對映到特定基因的讀取次數，其中 Ci ~ 二項式(ni, pi)，I = 1,2，其中 ni 表示對映讀取的總數，pi 表示讀取來自該基因的機率。我們定義 M = log2C1 - log2C2，以及 A = (log2C1 + log2C2)/2。可以證明，在隨機抽樣假設下，給定 A = a（a 是 A 的觀測值）的 M 的條件分佈遵循近似正態分佈。對於 MA 圖上的每個基因，我們對 H0: p1 = p2 與 H1: p1 ≠ p2 進行假設檢驗。然後，可以基於條件正態分佈分配一個 P 值。
基於 MA 圖的技術重複方法
儘管據報道測序平臺具有低背景噪聲，但技術重複仍然可以提供資訊，用於質量控制和估計由於不同機器或平臺造成的變異。基於 MA 圖的另一種方法透過比較資料中的技術重複（如果可用）來估計噪聲水平。在這種方法中，首先在兩個技術重複的 MA 圖上沿 A 軸應用滑動視窗，以估計對應於不同表達水平的隨機變異。透過區域性迴歸對強度依賴噪聲水平進行平滑估計，並在正態分佈假設下轉換為以 A 為條件的 M 的區域性標準差。然後，區域性標準差用於識別兩個樣本之間基因表達的差異。
多重檢驗校正
對於上述方法，針對每個基因計算的 P 值會調整為 Q 值，以進行多重檢驗校正。使用者可以設定 P 值或錯誤發現率 (FDR) 閾值來識別差異表達基因。
選擇資料後，您可以應用 R 包 DEGseq。有五種不同的方法來執行該程式。它們是 DEGexp、DEGseq、readGeneExp、samWrapper 和 getGeneExp。
描述
當用戶已經擁有基因表達值（例如對映到每個基因的讀取次數）時，此函式用於識別差異表達基因。
用法
DEGexp( geneExpFile1, geneCol1=1, expCol1=2, depth1=rep(0, length(expCol1) ), groupLabel1=" geneExpFile2, geneCol2=1, expCol2=2, depth2=rep(0, length(expCol2) ), groupLabel2=" header=TRUE, sep="", method=c("LRT", "CTR", "FET", "MARS", "MATR", "FC"), pValue=1e-3, zScore=4, qValue=1e-3, foldChange=4, thresholdKind=1, outputDir="none", normalMethod=c("none", "loess", "median"), replicate1="none", geneColR1=1, expColR1=2, depthR1=rep(0, length(expColR1) ), replicate2="none", geneColR2=1, expColR2=2, depthR2=rep(0, length(expColR2) ), rawcount=TRUE )
引數
| 引數 | 描述 |
|---|---|
| geneExpFile1 | 包含樣本 1（或 method="CTR" 時的重複 1）重複的基因表達值的檔案。 |
| geneCol1 | geneExpFile1 中的基因 ID 列。 |
| expCol1 | geneExpFile1 中樣本 1（數值向量）重複的表達值列。注意：每列對應於樣本 1 的一個重複。 |
| depth1 | 樣本 1（數值向量）每個重複的唯一對映到基因組的讀取總數，預設值：將對映到所有註釋基因的讀取總數作為每個重複的深度。 |
| groupLabel1 | 圖中組 1 的標籤。 |
| GeneExpFile2 | 包含樣本 2（或 method="CTR" 時的重複 2）重複的基因表達值的檔案。 |
| geneCol2 | geneExpFile2 中的基因 ID 列。 |
| expCol2 | geneExpFile2 中樣本 2（數值向量）重複的表達值列。注意：每列對應於樣本 2 的一個重複。 |
| depth2 | 樣本 2（數值向量）每個重複的唯一對映到基因組的讀取總數，預設值：將對映到所有註釋基因的讀取總數作為每個重複的深度。 |
| groupLabel2 | 圖中組 2 的標籤。 |
| header | 一個邏輯值，指示 geneExpFile1 和 geneExpFile2 是否包含變數名稱作為第一行。 |
| sep | 欄位分隔符字元。如果 sep = ""（read.table 的預設值），則分隔符為空格，即一個或多個空格、製表符、換行符或回車符。 |
| method | 識別差異表達基因的方法。可能的方法包括：• "LRT"：似然比檢驗 (Marioni 等人，2008)，• "CTR"：檢查技術重複之間的變異是否可以用隨機抽樣模型 (Wang 等人，2009) 解釋，• "FET"：Fisher 精確檢驗 (Joshua 等人，2009)，• "MARS"：基於 MA 圖的隨機抽樣模型方法 (Wang 等人，2009)，• "MATR"：基於 MA 圖的技術重複方法 (Wang 等人，2009)，• "FC"：MA 圖上的倍數變化閾值。 |
| pValue | p 值閾值（對於方法：LRT、FET、MARS、MATR）。僅在 thresholdKind=1 時使用。 |
| zScore | z 分數閾值（對於方法：MARS、MATR）。僅在 thresholdKind=2 時使用。 |
| qValue | q 值閾值（對於方法：LRT、FET、MARS、MATR）。僅在 thresholdKind=3 或 thresholdKind=4 時使用。 |
| thresholdKind | 閾值型別。可能的型別包括：• 1：p 值閾值，• 2：z 分數閾值，• 3：q 值閾值 (Benjamini 等人，1995)，• 4：q 值閾值 (Storey 等人，2003)。 |
| foldChange | MA 圖上的倍數變化閾值（對於方法：FC）。 |
| outputDir | 輸出目錄。 |
| normalMethod | 歸一化方法："none"、"loess"、"median"。建議："none"。 |
| replicate1 | 包含重複批次 1（僅在 method="MATR" 時使用）的基因表達值的檔案。注意：重複 1 和重複 2 是樣本的兩個（組）技術重複。 |
| GeneColR1 | 重複批次 1（僅在 method="MATR" 時使用）的表達檔案中的基因 ID 列。 |
| expColR1 | 重複批次 1（數值向量）（僅在 method="MATR" 時使用）的表達檔案中的表達值列。 |
| depthR1 | 重複批次 1（數值向量）中每個重複的唯一對映到基因組的讀取總數，預設值：將對映到所有註釋基因的讀取總數作為每個重複的深度（僅在 method="MATR" 時使用）。 |
| ReplicateLabel1 | 圖中重複批次 1 的標籤（僅在 method="MATR" 時使用）。 |
| replicate2 | 包含重複批次 2（僅在 method="MATR" 時使用）的基因表達值的檔案。注意：重複 1 和重複 2 是樣本的兩個（組）技術重複。 |
| geneColR2 | 重複批次 2（僅在 method="MATR" 時使用）的表達檔案中的基因 ID 列。 |
| expColR2 | 重複批次 2（數值向量）（僅在 method="MATR" 時使用）的表達檔案中的表達值列。 |
| depthR2 | 重複批次 2（數值向量）中每個重複的唯一對映到基因組的讀取總數，預設值：將對映到所有註釋基因的讀取總數作為每個重複的深度（僅在 method="MATR" 時使用）。 |
| ReplicateLabel2 | 圖中重複批次 2 的標籤（僅在 method="MATR" 時使用）。 |
| rawCount | 一個邏輯值，指示基因表達值是基於原始讀取計數還是歸一化值。 |
示例
> library(DEGseq)
> geneExpFile <- system.file("data", "GeneExpExample5000.txt",
+ package = "DEGseq")
> if (geneExpFile == "") {
+ zipFile <- system.file("data", "Rdata.zip", package = "DEGseq")
+ if (zipFile != "") {
+ unzip(zipFile, "GeneExpExample5000.txt", exdir = tempdir())
+ geneExpFile <- file.path(tempdir(), "GeneExpExample5000.txt")
+ }
+ }
> layout(matrix(c(1, 2, 3, 4, 5, 6), 3, 2, byrow = TRUE))
> par(mar = c(2, 2, 2, 2))
> DEGexp(geneExpFile1 = geneExpFile, expCol1 = c(7, 9, 12, 15, 18), groupLabel1 = "kidney", geneExpFile2 = geneExpFile, expCol2 = c(8, 10, 11, 13, 16), groupLabel2 = "liver", method = "MARS")
Please wait...
geneExpFile1: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/GeneExpExample5000.txt
gene id column in geneExpFile1: 1
expression value column(s) in geneExpFile1: 7 9 12 15 18
total number of reads uniquely mapped to genome obtained from sample1: 345504 354981 334557
geneExpFile2: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/GeneExpExample5000.txt
gene id column in geneExpFile2: 1
expression value column(s) in geneExpFile2: 8 10 11 13 16
total number of reads uniquely mapped to genome obtained from sample2: 274430 274486 264999
method to identify differentially expressed genes: MARS
pValue threshold: 0.001
output directory: none
The outputDir is not specified!
Please wait …
Identifying differentially expressed genes ...
Please wait patiently ...
output …
The results can be observed in directory: none
描述
此函式用於識別來自 RNA 測序資料的差異表達基因。它將來自 RNA 測序資料的兩個樣本的唯一對映讀取與基因註釋作為輸入。因此，使用者應提前將讀取（從樣本的測序文庫獲得）對映到相應的基因組。
用法
DEGseq ( mapResultBatch1, mapResultBatch2, fileFormat="bed", readLength=32, strandInfo=FALSE, refFlat, groupLabel1="group1", groupLabel2="group2", method=c("LRT", "CTR", "FET", "MARS", "MATR", "FC"), pValue=1e-3, zScore=4, qValue=1e-3, foldChange=4, thresholdKind=1, outputDir="none", normalMethod=c("none", "loess", "median"), depthKind=1, replicate1="none", replicate2="none", replicateLabel1="replicate1", replicateLabel2="replicate2" )
引數
| 引數 | 描述 |
|---|---|
| mapResultBatch1 | 樣本 1（或 method="CTR" 時的重複 1）技術重複的唯一對映結果檔案。 |
| MapResultBatch2 | 樣本 2（或 method="CTR" 時的重複 2）技術重複的唯一對映結果檔案。 |
| fileFormat | 檔案格式："bed" 或 "eland"。 "bed" 格式示例：chr12 7 38 readID 2 "eland" 格式示例：readID chr12.fa 7 U2 F 注意：欄位分隔符字元為 TAB。並且檔案必須遵循其中一個示例的格式。 |
| ReadLength | 讀取的長度（僅在 fileFormat="eland" 時使用）。strandInfo 在 cDNA 克隆過程中是否保留了鏈資訊。 "TRUE"：保留，"FALSE"：未保留。 |
| RefFlat | UCSC refFlat 格式的基因註釋檔案。 |
| GroupLabel1 | 圖中組 1 的標籤。 |
| GroupLabel2 | 圖中組 2 的標籤。 |
| Method | 識別差異表達基因的方法。可能的方法包括： "LRT"：似然比檢驗，"CTR"：檢查兩個技術重複之間的變異是否可以用隨機抽樣模型解釋，"FET"：Fisher 精確檢驗，"MARS"：基於 MA 圖的隨機抽樣模型方法，"MATR"：基於 MA 圖的技術重複方法，"FC"：MA 圖上的倍數變化閾值。 |
| pValue | p 值閾值（對於方法：LRT、FET、MARS、MATR）。僅在 thresholdKind=1 時使用。 |
| zScore | z 分數閾值（對於方法：MARS、MATR）。僅在 thresholdKind=2 時使用。 |
| qValue | q 值閾值（對於方法：LRT、FET、MARS、MATR）。僅在 thresholdKind=3 或 thresholdKind=4 時使用。 |
| ThresholdKind | 閾值型別。可能的型別包括：• 1：p 值閾值，• 2：z 分數閾值，• 3：q 值閾值，• 4：q 值閾值。 |
| foldChange | MA 圖上的倍數變化閾值（對於方法：FC）。 |
| outputDir | 輸出目錄。 |
| normalMethod | 歸一化方法："none"、"loess"、"median"。建議："none"。 |
| DepthKind | 1- 以每個複製樣本中唯一比對到基因組的總讀取數作為深度，0: 以每個複製樣本中唯一比對到所有註釋基因的總讀取數作為深度。 我們建議使用 depthKind=1，尤其是在註釋檔案中的基因是所有基因的一部分時。 |
| replicate1 | 包含來自複製批次 1 的唯一比對讀取的檔案（僅在 method="MATR" 時使用）。 |
| replicate2 | 包含來自複製批次 2 的唯一比對讀取的檔案（僅在 method="MATR" 時使用）。 |
| ReplicateLabel1 | 圖中重複批次 1 的標籤（僅在 method="MATR" 時使用）。 |
| ReplicateLabel2 | 圖中重複批次 2 的標籤（僅在 method="MATR" 時使用）。 |
示例
> kidneyR1L1 <- system.file("data", "kidneyChr21.bed.txt", package = "DEGseq") > liverR1L2 <- system.file("data", "liverChr21.bed.txt", package = "DEGseq") > refFlat <- system.file("data", "refFlatChr21.txt", package = "DEGseq") > mapResultBatch1 <- c(kidneyR1L1) > mapResultBatch2 <- c(liverR1L2) > outputDir <- file.path(tempdir(), "DEGseqExample") > DEGseq(mapResultBatch1, mapResultBatch2, fileFormat = "bed",refFlat = refFlat, outputDir = outputDir, method = "LRT")
請稍候...
mapResultBatch1: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt mapResultBatch2: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/liverChr21.bed.txt 檔案格式: bed refFlat: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/refFlatChr21.txt 在計算比對到基因的讀取數時，忽略鏈資訊！ 正在計算比對到每個基因的讀取數 ... 這將需要幾分鐘，請耐心等待！
請稍候...
樣本檔案: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt 正在計算比對到每個基因的讀取數。這將需要幾分鐘。
請稍候 …
總共 259 個獨特基因已處理 0 個讀取 (kidneyChr21.bed.txt) 已處理 10000 個讀取 (kidneyChr21.bed.txt) 已處理 20000 個讀取 (kidneyChr21.bed.txt) 已處理 30000 個讀取 (kidneyChr21.bed.txt) 已處理 34304 個讀取 (kidneyChr21.bed.txt) 總共耗時 0.328000 秒！
請稍候...
樣本檔案: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/liverChr21.bed.txt 正在計算比對到每個基因的讀取數。這將需要幾分鐘。
請稍候 …
總共 259 個獨特基因已處理 0 個讀取 (liverChr21.bed.txt) 已處理 10000 個讀取 (liverChr21.bed.txt) 已處理 20000 個讀取 (liverChr21.bed.txt) 已處理 30000 個讀取 (liverChr21.bed.txt) 已處理 30729 個讀取 (liverChr21.bed.txt) 總共耗時 0.297000 秒！
請稍候...
geneExpFile1: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample/group1.exp geneExpFile1 中的基因 ID 列: 1 geneExpFile1 中的表達值列: 2 從樣本 1 獲取的唯一比對到基因組的總讀取數: 34304 geneExpFile2: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample/group2.exp geneExpFile2 中的基因 ID 列: 1 geneExpFile2 中的表達值列: 2 從樣本 2 獲取的唯一比對到基因組的總讀取數: 30729 用於識別差異表達基因的方法: LRT p 值閾值: 0.001 輸出目錄: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample
請稍候 …
正在識別差異表達基因 ... 請耐心等待 ... 輸出 ...
已完成 ... 結果可在以下目錄中檢視: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGs
getGeneExp
[edit | edit source]描述
此函式用於計算每個基因的讀取數並計算 RPKM。它以來自樣本的 RNA-seq 資料的唯一比對讀取以及基因註釋檔案作為輸入。因此，使用者應提前將讀取 (從樣本的測序文庫獲得) 比對到相應的基因組。
用法
getGeneExp( mapResultBatch, fileFormat="bed", readLength=32, strandInfo=FALSE, refFlat, output=paste(mapResultBatch[1],".exp",sep=""), min.overlapPercent= 1 )
引數
| 引數 | 描述 |
|---|---|
| mapResultBatch | 包含樣本的唯一比對結果檔案的向量。注意: 樣本可以具有多個技術重複。 |
| fileFormat | 檔案格式："bed" 或 "eland"。 "bed" 格式示例：chr12 7 38 readID 2 "eland" 格式示例：readID chr12.fa 7 U2 F 注意：欄位分隔符字元為 TAB。並且檔案必須遵循其中一個示例的格式。 |
| readLength | 讀取的長度 (僅在 fileFormat="eland" 時使用)。 strandInfo 鏈資訊在 cDNA 克隆過程中是否保留。 • "TRUE" : 保留， • "FALSE": 未保留。 |
| refFlat | UCSC refFlat 格式的基因註釋檔案。 |
| output | 輸出檔案。 |
| min.overlapPercent | 讀取與外顯子重疊長度與讀取自身長度的最小百分比，用於計算外顯子的讀取數。應 <=1。0: 讀取與外顯子之間至少重疊 1 個鹼基。 |
示例
> kidneyR1L1 <- system.file("data", "kidneyChr21.bed.txt", package = "DEGseq") > refFlat <- system.file("data", "refFlatChr21.txt", package = "DEGseq") > mapResultBatch <- c(kidneyR1L1) > output <- file.path(tempdir(), "kidneyChr21.bed.exp") > getGeneExp(mapResultBatch, refFlat = refFlat, output = output) 請稍候... 樣本檔案: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt 正在計算比對到每個基因的讀取數。這將需要幾分鐘。請稍候 ... 總共 259 個獨特基因已處理 0 個讀取 (kidneyChr21.bed.txt) 已處理 10000 個讀取 (kidneyChr21.bed.txt) 已處理 20000 個讀取 (kidneyChr21.bed.txt) 已處理 30000 個讀取 (kidneyChr21.bed.txt) 已處理 34304 個讀取 (kidneyChr21.bed.txt) 總共耗時 0.328000 秒！ > exp <- readGeneExp(file = output, geneCol = 1, valCol = c(2, + 3), label = c("raw count", "RPKM")) > exp[30:32, ] raw count RPKM C21orf131 0 0.000 C21orf15 0 0.000 C21orf2 51 665.789
readGeneExp
[edit | edit source]描述
此方法用於將基因表達值從檔案讀取到 R 工作空間中的矩陣。以便矩陣可以用作其他軟體包的輸入，例如 edgeR。該方法的輸入是一個包含基因表達值的檔案。
用法
readGeneExp(file, geneCol=1, valCol=2, label = NULL, header=TRUE, sep="")
引數
引數描述 file 包含基因表達值的檔案。 GeneCol 檔案中的基因 ID 列。 valCol 檔案中要讀取的表達值列。 label 列的標籤。 Header 一個邏輯值，指示檔案是否將其第一行作為變數名稱。 sep 欄位分隔符字元。如果 sep = "" (read.table 的預設值)，則分隔符是空白，即一個或多個空格、製表符、換行符或回車符。
示例
> geneExpFile <- system.file("data", "GeneExpExample1000.txt", + package = "DEGseq") > exp <- readGeneExp(file = geneExpFile, geneCol = 1, valCol = c(7, + 9, 12, 15, 18, 8, 10, 11, 13, 16)) > exp[30:32, ]
R1L1Kidney
R1L3Kidney
R1L7Kidney
R2L2Kidney
R2L6Kidney ENSG00000188976 73 77 68 70 82 ENSG00000187961 15 15 13 12 15 ENSG00000187583 1 1 3 0 3
R1L2Liver
R1L4Liver
R1L6Liver
R1L8Liver
R2L3Liver ENSG00000188976 34 56 45 55 42 ENSG00000187961 8 13 11 12 20 ENSG00000187583 0 1 0 0 2
samWrapper
[edit | edit source]描述
此函式是 samr 中函式的包裝器。它用於識別兩組具有多個重複的樣本或來自不同個體 (例如，疾病樣本與對照樣本) 的兩組樣本的差異表達基因。
用法
samWrapper( geneExpFile1, geneCol1=1, expCol1=2, measure1=rep(1, length(expCol1) ), geneExpFile2, geneCol2=1, expCol2=2, measure2=rep(2, length(expCol2) ), header=TRUE, sep="", paired=FALSE, s0=NULL, s0.perc=NULL, nperms=100, testStatistic= c("standard","wilcoxon"), max.qValue=1e-3, in.foldchange=logged2=FALSE, output )
引數
引數描述 geneExpFile1 包含組 1 的基因表達值的檔案。 geneCol1 geneExpFile1 中的基因 ID 列。 expCol1 geneExpFile1 中的表達值列。參見示例。 measure1
numeric vector of outcome measurements for group1. like c(1,1,1...) when paired=FALSE, or like c(-1,-2,-3,...) when paired=TRUE.
geneExpFile2 包含組 2 的基因表達值的檔案。 geneCol2 geneExpFile2 中的基因 ID 列。 ExpCol2 geneExpFile2 中的表達值列。參見示例。 Measure2 組 2 的結果測量數值向量。如果 paired=FALSE，則類似 c(2,2,2...)，如果 paired=TRUE，則類似 c(1,2,3,...)。 header 一個邏輯值，指示 geneExpFile1 和 geneExpFile2 是否將其第一行作為變數名稱。 sep 欄位分隔符字元。如果 sep = "" (read.table 的預設值)，則分隔符是空白，即一個或多個空格、製表符、換行符或回車符。 paired 一個邏輯值，指示樣本是否配對。 s0 測試統計量分母的可交換性因子；預設是自動選擇。 s0.perc 用於 s0 的標準偏差值的百分位數；預設是自動選擇；-1 表示 s0=0 (不同於 s0.perc=0，表示 s0=標準偏差值的第零百分位數= sd 值的最小值。 Nperms 用於估計錯誤發現率的置換次數。 TestStatistic 在兩類非配對情況下使用的測試統計量。 "standard" (t 統計量) 或 "wilcoxon" (雙樣本 Wilcoxon 或 Mann-Whitney 檢驗)。建議使用 "standard"。 max.qValue 所需的最大 q 值；應 <1。 min.foldchange 所需的最小倍數變化；應 >1。預設值為零，表示不應用倍數變化標準。 logged2 一個邏輯值，指示表達值是否為 log2。 output 輸出檔案。
示例
> geneExpFile <- system.file("data", "GeneExpExample1000.txt", + package = "DEGseq") > set.seed(100) > geneExpFile1 <- geneExpFile > geneExpFile2 <- geneExpFile > output <- file.path(tempdir(), "samWrapperOut.txt") > expCol1 = c(7, 9, 12, 15, 18) > expCol2 = c(8, 10, 11, 13, 16) > measure1 = c(-1, -2, -3, -4, -5) > measure2 = c(1, 2, 3, 4, 5) > samWrapper(geneExpFile1 = geneExpFile1, geneCol1 = 1, expCol1 = expCol1, measure1 = measure1, geneExpFile2 = geneExpFile2, geneCol2 = 1, expCol2 = expCol2, measure2 = measure2, nperms = 100, min.foldchange = 2, max.qValue = 1e-04, output = output, paired = TRUE)
參考文獻
[edit | edit source]Jiang,H. and Wong,W.H. (2009) Statistical inferences for isoform expression in RNASeq. Bioinformatics, 25, 1026–1032.
Likun Wang, Zhixing Feng, Xi Wang, Xiaowo Wang, and Xuegong Zhang. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data Bioinformatics Advance Access published on October 24, 2009, DOI 10.1093/bioinformatics/btp612.
王立昆; 馮志興; 王曦; 王曉沃; 張雪峰. DEGseq 手冊. 清華大學. 北京, 中國. 2009(B). 可從 <http://www.bioconductor.org/packages/devel/bioc/manuals/DEGseq/man/DEGseq.pdf> 訪問, 訪問日期: 2009 年 11 月 18 日
王立昆; 王曦. 如何使用 Degseq 包. 清華大學自動化系/ TNLIST 生物資訊學實驗室和生物資訊學部門. 北京, 中國. 2009(A). 可從 <http://bioinfo.au.tsinghua.edu.cn/software/degseq/DEGseq.pdf> 訪問, 訪問日期: 2009 年 11 月 18 日