DESeq2 paired design
1
0
Entering edit mode
Raj • 0
@raj-9784
Last seen 3.1 years ago
USA

Hi Mike,

I have a question about the paired design. I have gone through vignette about linear combinations but could not achieve Pair adjustment with another factor. I am getting model being not full rank. Could you please help? Below is my design and what I try to achieve,


prep = c("L","D","D","D","L","L","D","L","D")
lb1 = rep(paste('Lib1','Tre',sep="_"),9)
lb2 = rep(paste('Lib2','Tre',sep="_"),9)

df.cd = data.frame(SampleID     = LETTERS[1:18],
                   Prep         = c(prep,prep),
                   cond         = c(lb1,lb2),
                   Pair         = factor(c(seq(1:9),seq(1:9))) )

rownames(df.cd) = df.cd$SampleID

all.equal(colnames(mat),rownames(df.cd))#[1] TRUE

dds = DESeqDataSetFromMatrix(countData = mat,
                            colData = df.cd,
                            design = ~Pair+Prep+cond)
****ERROR** Please help***


dds2  = DESeq(dds)

res   = results(dds2,contrast=c("cond","Lib1_Tre",'Lib2_Tre'))

sessionInfo( )
R version 4.1.0 (2021-05-18)
DESeq2_1.32.0
DESeq2 • 1.5k views
ADD COMMENT •
0
Entering edit mode

Please let me know if anything is not clear.

ADD REPLY •
1
Entering edit mode
@mikelove
Last seen 21 hours ago
United States

For questions about statistical design and analysis, I recommend consulting with a local statistician or someone familiar with linear models in R. I unfortunately have to limit my time to software related questions.
ADD COMMENT •
0
Entering edit mode

Hi Mike, I reread your vignette again and did the trick of refactoring like the following but still the issue.

prep = c("L","D","D","D","L","L","D","L","D") lb1 = rep(paste('Lib1','Tre',sep="_"),9) lb2 = rep(paste('Lib2','Tre',sep="_"),9)

df.cd = data.frame(SampleID = LETTERS[1:18], Prep = c(prep,prep), cond = c(lb1,lb2), Pair = factor(c(seq(1:9),seq(1:9))) )

rownames(df.cd) = df.cd$SampleID

df.cd = df.cd[order(df.cd$Prep,df.cd$Pair),] # reordering to distinguish Prep groups

df.cd$ind.n = factor(c(rep(rep(1:5),each=2),rep(rep(1:4),each=2)))

all.equal(colnames(mat),rownames(df.cd))#[1] TRUE

analysis of group-specific condition effects, while controlling for differences in individual.

dds = DESeqDataSetFromMatrix(countData = mat, colData = df.cd, design = ~Prep + Prep:ind.n + Prep:cond)

ADD REPLY •
0
Entering edit mode

Hi Mike,

I went through the vignette section "Levels without samples" and worked out my own model matrix. Now, the question is about the contrasting levels of cond. Please suggest?

prep = c("L","D","D","D","L","L","D","L","D") lb1 = rep("Lib1",9) lb2 = rep("Lib2",9)

df.cd = data.frame(SampleID = LETTERS[1:18], Prep = factor(c(prep,prep)), cond = factor(c(lb1,lb2)), Pair = factor(c(seq(1:9),seq(1:9))) )

rownames(df.cd) = df.cd$SampleID

reordering for the purposes of better creation of column which can distinguish individuals nested within a group

df.cd = df.cd[order(df.cd$Prep,df.cd$Pair),]

creating a column

df.cd$ind.n = factor(c(rep(rep(1:5),each=2),rep(rep(1:4),each=2)))

making sure matrix and condition data frame are ordered same

all.equal(colnames(mat),rownames(df.cd))

[1] TRUE

dds = DESeqDataSetFromMatrix(countData = mat, colData = df.cd, design = ~Prep + Prep:ind.n + Prep:cond) m1 = model.matrix(~Prep + Prep:ind.n + Prep:cond, df.cd)

all.zero = apply(m1, 2, function(x) all(x==0)) all.zero idx = which(all.zero) m1 = m1[,-idx]

dds = DESeqDataSetFromMatrix(countData = mat, colData = df.cd, design = m1)

dds2 = DESeq(dds)

---wanted to contrast two libs....any suggestions?

res = results(dds2, contrast=list("Lib1",'Lib2'))

--PS: Please ignore the below answer as it is an attempt at better formatting

ADD REPLY •
1
Entering edit mode

Sorry, I really don't have time to help people work out their designs and contrasts to match their experimental goals, I have to limit myself to issues with the software, not statistical support. You can try working with a statistician perhaps or someone familiar with linear models in R.

ADD REPLY •
1
Entering edit mode

@Raj, You can post it over at biostars.org which has a much larger audience than this forum here and is intended for all kinds of bioinfo-related questions rather than Bioc support forum which is mainly about technical help with the packages. Just add a link to this thread here for reference.

ADD REPLY •
0
Entering edit mode

Hi ATpoint, Thank you for the suggestion. I will post it in biostars as well.

Hi Mike,

 I have a specific question about the vignette section. I read from vignette.

"Above, the terms grpX.cndB and grpY.cndB give the group-specific condition effects, in other words, the condition B vs A effect for group X samples, and likewise for group Y samples. These terms control for all of the six individual effects. "

The terms does not have cndA but vignette says condition effects is achieved for the condition B vs A effect for group X samples, and likewise for group Y samples. These terms control for all of the six individual effects. Could you please explain about this vignette snippet?

Thank you very much for your time

ADD REPLY •
1
Entering edit mode

The way that R's model.matrix names these linear model terms, they are just named with cndB. But the meaning is B vs A. For main effects, we rename the term so it has the "vs" and the reference level for clarity, but for higher order terms it is directly from R's model.matrix.

Again, for understanding the meaning of terms in a linear model, that is something you need to consult with someone at your institute, if you aren't familiar, and if you want to get the analysis right. It's not the responsibility of the open source software developer to explain these details of linear model terms.

ADD REPLY •
0
Entering edit mode

Hi Mike, I reread your vignette again and did the trick of refactoring like the following but still the issue. ```

prep = c("L","D","D","D","L","L","D","L","D")
lb1 = rep(paste('Lib1','Tre',sep="_"),9)
lb2 = rep(paste('Lib2','Tre',sep="_"),9)

df.cd = data.frame(SampleID     = LETTERS[1:18],
                   Prep         = c(prep,prep),
                   cond         = c(lb1,lb2),
                   Pair         = factor(c(seq(1:9),seq(1:9))) )

rownames(df.cd) = df.cd$SampleID

df.cd = df.cd[order(df.cd$Prep,df.cd$Pair),] # reordering to distinguish Prep groups

df.cd$ind.n = factor(c(rep(rep(1:5),each=2),rep(rep(1:4),each=2)))
as.data.frame(df.cd)

all.equal(colnames(mat),rownames(df.cd))#[1] TRUE

#analysis of group-specific condition effects, while controlling for differences in individual.
dds = DESeqDataSetFromMatrix(countData = mat,
                            colData = df.cd,
                            design = ~Prep + Prep:ind.n + Prep:cond)
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 619 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6