|
DESeq2 接受raw count的定量表格,，然后根據(jù)樣本分組進行差異分析,，具體步驟如下 1. 讀取數(shù)據(jù)讀取基因的表達量表格和樣本的分組信息兩個文件，其中表達量的文件示例如下 gene_id ctrl-1 ctrl-2 ctrl-3 case-1 case-2 case-3
geneA 14 0 11 4 0 12
geneB 125 401 442 175 59 200每一行為一個基因,，每一列代表一個樣本,。 分組信息的文件示例如下 sample condition
ctrl-1 control
ctrl-2 control
ctrl-3 control
case-1 case
case-2 case
case-3 case第一列為樣本名，第二列為樣本的分組信息,。 讀取文件的代碼如下 # 讀取表達量的表格
count <- read.table(
"gene.counts.tsv",
header=T,
sep="\t",
row.names=1,
comment.char="",
check.names=F)
# 預處理,，過濾低豐度的數(shù)據(jù)
countData <- count[apply(count, 1, sum) > 0 , ]
# 讀取樣本分組信息
colData <- read.table(
"sample.group.tsv",
header=T,
sep="\t",
row.names=1,
comment.char="",
check.names=F)
# 構(gòu)建DESeq2中的對象
dds <- DESeqDataSetFromMatrix(
countData = countData,
colData = colData,
design = ~ condition)
# 指定哪一組作為control
dds$condition <- relevel(dds$condition, ref = "control")在讀取數(shù)據(jù)的過程中，有兩點需要注意,，第一個就是根據(jù)表達量對基因進行過濾,，通常是過濾低表達量的基因，這一步是可選的,，閾值可以自己定義,；另外一個就是指定哪一組作為control組，在計算log2FD時 ,，需要明確control組,，默認會字符串順序?qū)Ψ纸M的名字進行排序，排在前面的作為control組,，這種默認行為選出的control可能與我們的實驗設計不同,，所以必須明確指定control組。 2. 歸一化計算每個樣本的歸一化系數(shù),，代碼如下 dds <- estimateSizeFactors(dds)
3. 估計基因的離散程度DESeq2假定基因的表達量符合負二項分布,，有兩個關(guān)鍵參數(shù),，總體均值和離散程度α值， 如下圖所示 這個 通過下列代碼估算每個基因的α值 dds <- estimateDispersions(dds)4. 差異分析代碼如下 dds <- nbinomWaldTest(dds)
res <- results(dds)為了簡化調(diào)用,，將第二部到第四部封裝到了 dds <- DESeq(dds)
res <- results(dds)
write.table(
res,
"DESeq2.diff.tsv",
sep="\t",
quote=F,
col.names = NA)在DESeq2差異分析的過程中,，已經(jīng)考慮到了樣本之間已有的差異,，所以可以發(fā)現(xiàn)，最終結(jié)果里的log2FD值和我們拿歸一化之后的表達量計算出來的不同, 示意如下 > head(results(dds)[, 1:2])
log2 fold change (MLE): condition B vs A
DataFrame with 6 rows and 2 columns
baseMean log2FoldChange
<numeric> <numeric>
gene1 7.471250 -0.8961954
gene2 18.145279 0.4222174
gene3 2.329461 -2.3216915
gene4 165.634256 -0.1974001
gene5 38.300621 1.3573162
gene6 7.904819 1.8384322提取歸一化之后的表達量，自己計算baseMean和logFoldChange, 示例數(shù)據(jù)包含了12個樣本,，其中前6個樣本為control, 后6個樣本為case , 代碼如下 > nor_count <- counts(dds, nor = T)
> res <- data.frame(
baseMean = apply(nor_count, 1, mean),
log2FoldChange = apply(nor_count, 1, function(t){ mean(t[7:12]) / mean(t[1:6])})
)
> head(res)
baseMean log2FoldChange
gene1 7.471250 0.5380191
gene2 18.145279 1.3404422
gene3 2.329461 0.1991979
gene4 165.634256 0.8719078
gene5 38.300621 2.5621035
gene6 7.904819 3.5365201對比DESeq2和自己計算的結(jié)果,，可以發(fā)現(xiàn)，baseMeans是一致的,，而log2Foldchange 差異很大,，甚至連數(shù)值的正負都發(fā)生了變化。 log2FD 反映的是不同分組間表達量的差異,，這個差異由兩部分構(gòu)成,，一種是樣本間本身的差異，比如生物學重復樣本間基因的表達量就有一定程度的差異,，另外一部分就是我們真正感興趣的,，由于分組不同或者實驗條件不同造成的差異。 用歸一化之后的數(shù)值直接計算出的log2FD包含了以上兩種差異,，而我們真正感興趣的只有分組不同造成的差異,，DESeq2在差異分析的過程中已經(jīng)考慮到了樣本本身的差異，其最終提供的log2FD只包含了分組間的差異,，所以會與手動計算的不同,。 ·end· |
|
|
