Analyzing RNA-seq data with DESeq2翻譯（1）

01-30

為了深入學習DESeq包以及差異基因表達分析，我開始對Bioconductor網站上的DESeq包的說明進行翻譯。以下是部分選譯內容。

Quick start

DESeq2分析的基本步驟如下：

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~batch + condition)ndds <- DESeq(dds)nres <- results(dds, contrast=c("condition","treated","control"))n

cts是count matrix

coldata是table of sample information

design是指明如何對樣本進行比較的信息。上述代碼中是為了在控制批次效應的同時研究condition間的效應。

數據輸入

Salmon, Sailfish, kallito等軟體輸出的transcript quantification轉錄本定量文件，可以用DESeqDataSetFromTximport導入。
htseq-count軟體輸出的結果可以用DESeqDataSetFromHTSeq導入。
RangedSummarizedExperiment用DESeqDataSet導入。

輸入數據

為何必須輸入非標準化（非均一化）的counts值？

EDSeq2要求輸入數據是高通量測序實驗得到的，整數值組成的矩陣。該矩陣中第i行和第j列的交點代表在樣品j中有多少條reads落在基因i上（就是說矩陣中樣品為列，基因為行）。

矩陣中的數據應當是沒有標準化後的，這樣才能保證DESeq2所用的統計模型的分析結果的準確性。DESeq2的統計模型會自動根據矩陣大小進行校正，因此根據矩陣大小轉化後的或者標準化（均一化）的矩陣不能作為DESeq2的輸入數據。

DESeqDataSet

DESeqDataSet是DESeq2包中儲存read counts以及統計分析過程中的數據的一個「對象」，在代碼中常表示為「dds」。

一個DESeqDataSet對象必須關聯相應的design公式。design公式指明了要對哪些變數進行統計分析。該公式（上文中的design = ~batch + condition）以短波浪字元開頭，中間用加號連接變數。design公式可以修改，但是相應的差異表達分析就需要重新做，因為design公式是用來估計統計模型的分散度以及log2 fold change的。

注意：要將你感興趣的變數放在design公式的後面，對照變數放在前面。

以下介紹4種不同pipiline下構建DESeqDataSet方法：

從轉錄本丰度文件以及tximport構建
從count matrix 構建
從htseq-count文件構建
從SummarizedExperiment對象構建

導入轉錄本丰度文件以及tximport文件生成的矩陣

該方法是最新的及推薦的方法。具體是將上游的快速轉錄本丰度定量軟體獲得的基因水平的count matrix，然後用tximport包導入，這種方法支持的軟體有：

Salmon
Sailfish
kailisto
RSEM

用以上軟體進行轉錄本丰度估計的優點有：

校正了潛在的不同樣本間基因長度導致的錯誤。
運行速度快，與alignment-based的方法相比需要的內存和硬碟空間少。
避免將那些比對到帶有同源序列的多個基因上的片段丟棄，從而增加敏感度。

以下我們演示如何將Salmon軟體導出的定量文件「quant.sf」導入到DESeq2中。

library("DESeq2")nddsTxi <- DESeqDataSetFormTximport(txi, colData = samples, design =~condition)n

導入count matrix矩陣

用DESeqDataSetFromMatrix命令導入count matrix。除count matrix外，使用者還需要將count matrix的列名單獨建立一個DataFrame，以及一個design formula。以下我們演示DESeqDataSetFromMatrix的用法，主要過程是從pasilla包中載入count data，讀取的count matrix命名為cts，樣品信息表命名為coldata。下邊演示的方法是從featureCounts的輸出結果中提取以上信息。

library("pasilla")nnpasCts <- system.file("extdata", "pasilla_gene_counts.tsv",package="pasilla",mustWork=TRUE)nnpasAnno <- system.file("extdata", "pasilla_sample_annotation.csv", packages="pasilla", mustWork=TRUE)nncts <- as.matrix(read.csv(pasCts, sep="t", row.names="gene_id"))nncoldata <- read.csv(pasAnno,row.names=1)nncoldata <- coldata[,c("condition", "type")]n

看一下cts和coldata的內容

head(cts,2)n

untreated1untreated2untreated3untreated4treated1treated2FBgn0000003000000FBgn000000892161767014088

coldatan

conditiontypetreated1fbtreatedsingle-readtreated2fbtreatedpaired-endtreated3fbtreatedpaired-enduntreated1fbuntreatedsingle-readuntreated2fbuntreatedsingle-readuntreated3fbuntreatedpaired-enduntreated4fbuntreatedpaired-end

注意count matirx 的列與coldata的行的順序一致，DESeq2默認這兩者間的順序是一致的。

上面的結果中可以看出cts和coldata中列和行的順序並不一致，另外coldata中的行名也多了「fb」兩個字母，需要調整順序並且去掉多餘的字母。

rownames(coldata) <- sub("fb", , rownames(coldata)) #去掉fbnall(rownames(coldata)%in%colnames(cts)) n[1] TRUE #此行為運行結果nall(rowname(coldata)==colname(cts))n[1]FALSE #此行為運行結果ncts <- cts[, rownames(coldata)] #直接把列名提換了嗎？該列的數據提換了嗎？nall(rowname(coldata)==colname(cts))n

有了count matrix（cts）和sample information（coldata），我們就可以構建DESeqDataSet矩陣了

library("DESeq2")ndds <- DESeqDataSetFromMatrix(countData=cts, colData = coldata, design = ~condition)nddsn

導入htseq-count數據

通過HTSeq軟體進行read count的結果數據，可以用DESeqDataSetFromHTSeqCount導入到DESeq2中。首先要將htseq-count軟體的結果文件的位置信息保存在一個變數中。下文中我們用pasilla包中的htseq-count結果文件作為例子。

directory <- systme.file("extdata", package = "pasilla", mustWork =TRUE)nnsampleFiles <- grep("treated",list.files(directory), value=TRUE) #用list.file指明which file to read，再用grep選擇包含"treated"的文件nnsampleCondition <- sub("(.*treated).*". "1", sampleFiles)#用sub函數從sample filename得到condition statusnnsampleTable <- data.frame(sampleName =sampleFiles, fileName =sampleFiles, condition =sampleCondition)nn#build the DESeqDataSetnnlibrary("DESeq2")nnddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory= directory, design=~condition)nnddsHTseq #顯示該矩陣n

class: DESeqDataSet
dim: 70463 7 metadata(1): version assays(1): counts rownames(70463): FBgn0000003:001 FBgn0000008:001 ... FBgn0261575:001 FBgn0261575:002 rowData names(0): colnames(7): treated1fb.txt treated2fb.txt ... untreated3fb.txt untreated4fb.txt colData names(1): condition

預過濾

雖然提前去掉low count基因不是運行DESeq進行差異表達分析必須的，但是提前過濾掉low count基因有兩個好處：減少dds矩陣的大小，提高DESeq運行的速度。這裡我們演示如何去掉read count在10以下的基因。注意：更嚴格的過濾會在DESeq的results功能中自動實現

keep <- rowSums(counts(dds)) >= 10ndds <- dds[keep,]n

關於因子的說明

默認情況下，R會安裝字母順序選擇一個level作為參考level。就是說，如果沒有給DEseq2函數指定要與哪個level進行比較（即指定哪個level是對照組），默認會安裝字母順序選擇排在最前面的level作為對照組。

如果要指定哪個level是對照組，可以直接在函數#results#的參數contrast中指明對照level（下文中會具體介紹），或者直接在factor level中指明。

注意，更改design公式中變數的factor level只能在運行DESeq2分析前，分析中或分析後都不能再進行更改。更改dds矩陣中的factor level可以用以下兩個方法：

#用factorndds$condition <- factor(dds$condition, levels = c("untreated", "treated"))n#或者用relevel，直接指定參考levelndds$condition <- relevel(dds$condition, ref = "untreated")n

如果需要從dds中提取一個子集，比如說從dds中去掉部分樣品，這時有可能會同時去掉所有樣品的一個或者多個level。這種情況下，可以用droplevles函數去掉沒有對應sample的level：

dds$condition <- droplevels(dds$condition)n

摺疊技術重複

DESeq2包中的collapseReplicates函數可以用於將技術重複的定量數據摺疊成一個樣品的。注意不能把這個函數用到生物學重複上，更多信息請查看collapseReplicates的幫助頁面。