Goal: To ensure "the sum of the counts for [my samples] is the same as the counts in the [samples] columns in ddsColl" after collapsing replicates using DESeq2.
Before performing differential expression (DE) analysis on my RNA-seq data stored in a DESeqDataSet object, I need to collapse the read counts observed for each of my technical replicates. This should be simple enough using the DESeq2 manual and vignettes available as reference.
This is how you collapse replicates:
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
# examine the colData of the collapsed data
colData(ddsColl)
# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
However, when I try this, I get the following error with no apparent resolution reported in the manual or on online forums:
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[, .... is not TRUE
How might I circumvent this error and continue on with the analysis?
Input files:
> head(cts)
KO1_P1_L001.Counts KO1_P1_L002.Counts KO1_P2_L001.Counts KO1_P2_L002.Counts KO1_P3_L001.Counts KO1_P3_L002.Counts KO2_P1_L001.Counts KO2_P1_L002.Counts KO2_P2_L001.Counts KO2_P2_L002.Counts KO2_P3_L001.Counts KO2_P3_L002.Counts WT_P1_L001.Counts WT_P1_L002.Counts WT_P2_L001.Counts WT_P2_L002.Counts WT_P3_L001.Counts WT_P3_L002.Counts
DDX11L1 1 0 1 1 2 1 0 1 2 1 3 1 7 2 1 1 0 0
WASH7P 124 144 128 151 102 118 39 41 103 87 105 106 125 120 13 26 104 136
MIR6859-1 8 6 4 5 4 6 2 2 7 2 7 6 6 6 2 2 2 9
MIR1302-2HG 1 4 0 0 0 2 2 2 2 2 0 4 1 2 0 3 2 2
MIR1302-2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
FAM138A 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> print(coldata)
cell_type condition sample run
KO1_P1_L001.Counts human longKO l1 long_P1_L001
KO1_P1_L002.Counts human longKO l1 long_P1_L002
KO1_P2_L001.Counts human longKO l2 long_P2_L001
KO1_P2_L002.Counts human longKO l2 long_P2_L002
KO1_P3_L001.Counts human longKO l3 long_P3_L001
KO1_P3_L002.Counts human longKO l3 long_P3_L002
KO2_P1_L001.Counts human shortKO s1 short_P1_L001
KO2_P1_L002.Counts human shortKO s1 short_P1_L002
KO2_P2_L001.Counts human shortKO s2 short_P2_L001
KO2_P2_L002.Counts human shortKO s2 short_P2_L002
KO2_P3_L001.Counts human shortKO s3 short_P3_L001
KO2_P3_L002.Counts human shortKO s3 short_P3_L002
WT_P1_L001.Counts human WT w1 wt_P1_L001
WT_P1_L002.Counts human WT w1 wt_P1_L002
WT_P2_L001.Counts human WT w2 wt_P2_L001
WT_P2_L002.Counts human WT w2 wt_P2_L002
WT_P3_L001.Counts human WT w3 wt_P3_L001
WT_P3_L002.Counts human WT w3 wt_P3_L002
My workflow:
> dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
> dds
class: DESeqDataSet
dim: 59368 18
metadata(1): version
assays(1): counts
rownames(59368): DDX11L1 WASH7P ... MT-TT MT-TP
rowData names(0):
colnames(18): KO1_P1_L001.Counts KO1_P1_L002.Counts ... WT_P3_L001.Counts WT_P3_L002.Counts
colData names(4): cell_type condition sample run
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> ddsColl
class: DESeqDataSet
dim: 59368 9
metadata(1): version
assays(1): counts
rownames(59368): DDX11L1 WASH7P ... MT-TT MT-TP
rowData names(0):
colnames(9): l1 l2 ... w2 w3
colData names(5): cell_type condition sample run runsCollapsed
> colData(ddsColl)
DataFrame with 9 rows and 5 columns
cell_type condition sample run runsCollapsed
<character> <factor> <character> <character> <character>
l1 human longKO l1 long_P1_L001 long_P1_L001,long_P1_L002
l2 human longKO l2 long_P2_L001 long_P2_L001,long_P2_L002
l3 human longKO l3 long_P3_L001 long_P3_L001,long_P3_L002
s1 human shortKO s1 short_P1_L001 short_P1_L001,short_P1_L002
s2 human shortKO s2 short_P2_L001 short_P2_L001,short_P2_L002
s3 human shortKO s3 short_P3_L001 short_P3_L001,short_P3_L002
w1 human WT w1 wt_P1_L001 wt_P1_L001,wt_P1_L002
w2 human WT w2 wt_P2_L001 wt_P2_L001,wt_P2_L002
w3 human WT w3 wt_P3_L001 wt_P3_L001,wt_P3_L002
> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
logical(0)
> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[, .... is not TRUE
The only difference I see between my DESeq2 workflow output is the result of matchFirstLevel (see below for expected result). Why might this be for my dataset, but not the example dataset?
Example workflow:
> dds <- makeExampleDESeqDataSet(m=12)
> dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
> dds$run <- paste0("run",1:12)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
condition sample run runsCollapsed
<factor> <factor> <character> <character>
sample1 B sample1 run10 run10,run11
sample2 A sample2 run5 run5
sample3 A sample3 run6 run6
sample4 A sample4 run1 run1,run4
sample5 B sample5 run7 run7
sample6 B sample6 run9 run9
sample7 A sample7 run3 run3,run12
sample8 B sample8 run8 run8
sample9 A sample9 run2 run2
> colnames(ddsColl)
[1] "sample1" "sample2" "sample3" "sample4" "sample5" "sample6" "sample7" "sample8" "sample9"
> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE
> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))