**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][1] and [vignettes][2] 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:**
> cts <- read.table("counts.txt", header=T, sep="\t", row.names="Geneid")
> coldata <- read.csv(file="annotationFile.csv", sep=",", row.names=1)
> 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])))
[1]: https://bioconductor.org/packages/devel/bioc/manuals/DESeq2/man/DESeq2.pdf
[2]: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html