I can't get the 2 by 2 comparisons that I need because the matrix is not full of rank
1
0
Entering edit mode
@user-24822
Last seen 3.1 years ago

Hello,

I have an RNA-seq experiment whose colData looks like that (with more samples but this is a representation of all possible groups we have):

Sample State AAV Sex group

C001CH0 KI GFP M KI_GFP

C001CH2 WT GFP M WT_GFP

C001CH4 KI JAK F KI_JAK2

I would like to get the comparison of the three groups at a time with a one-way-anova approach from the column 'group'. Also, I am interested on studying these comparisons with the results function: list(c("AAV_JAK2_vs_GFP","StateKI.AAVJAK2"))) to get KI-GFP vs KI-JAK2, contrast=c("State","KI","WT")) to get GFP-WT vs GFP-KI and name="StateKI.AAVJAK2" to get WT-GFP vs KI-JAK2. The Sex column will be included in the formula to modulate the batch effect. Then, the design that I am considering for this is the following:

#1
dds_anova <- DESeqDataSetFromMatrix(countData = all_count,
                      colData = design_all,
                      design= ~ Sex + group)
dds = DESeq(dds_anova, test = "LRT", reduced = ~ 1) #to test the 3 groups all together
#This one apparently works

#2
dds_2by2 <- DESeqDataSetFromMatrix(countData = all_count,
                      colData = design_all,
                      design= ~ Sex + State + AAV + State:AAV) #to test the 2 by 2 after in the results

Error in checkFullRank(modelMatrix) : 
  the model matrix is not full rank, so the model cannot be fit as specified.
  Levels or combinations of levels without any samples have resulted in
  column(s) of zeros in the model matrix.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')

sessionInfo( )

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

DESeq2_1.28.1

I understand this is because deseq2 would expect also the WT-JAK2 existing when it creates automatically the possible contrasts. But that's not the case in my design and then, the matrix is not full of rank. I could split the matrix and the colData to have only the individual 2 by 2 comparisons individually. However, I don't know if the normalisation of the data in the full matrix containing all the information for the anova in comparison with the reduced version of the matrix could affect the statistics. Hence, I am not really being consistent with the normalisation to study the different situations.

  • So: what do you think it could be the best approach to solve or hatch the problem? How could I avoid splitting the data to use all the information available in the matrix in order to normalise before obtaining the results for the different contrasts?

Thank you in advance,

Miriam

DESeq2 • 817 views
ADD COMMENT
0
Entering edit mode
ATpoint ★ 4.0k
@atpoint-13662
Last seen 1 hour ago
Germany

If it is not full rank then sex is at some point nested with group, meaning that for some group you have either only M or F in sex. Can you show the full colData (please highlight with markdown using the CODE button)?

ADD COMMENT
0
Entering edit mode

this is the full colData:

Sample State AAV Sex group
C001CH0 KI  GFP M   KI_GFP
C001CH2 WT  GFP M   WT_GFP
C001CH3 KI  GFP M   KI_GFP
C001CH4 KI  JAK M   KI_JAK2
C001CH5 WT  GFP M   WT_GFP
C001CH6 KI  GFP M   KI_GFP
C001CH7 KI  JAK M   KI_JAK2
C001CH8 KI  GFP F   KI_GFP
C001CH9 KI  JAK F   KI_JAK2
C001CHA KI  GFP M   KI_GFP
C001CHB KI  JAK M   KI_JAK2
C001CHC WT  GFP F   WT_GFP
C001CHE KI  JAK F   KI_JAK2
C001CHF WT  GFP F   WT_GFP
C001CHG KI  GFP F   KI_GFP
C001CHH KI  JAK F   KI_JAK2
C001CHI KI  GFP F   KI_GFP
C001CHJ KI  JAK F   KI_JAK2
C001CHK WT  GFP M   WT_GFP
C001CHL KI  GFP F   KI_GFP
C001CHM KI  JAK F   KI_JAK2
C001CHN KI  GFP M   KI_GFP
C001CHO KI  JAK M   KI_JAK2
C001CHP KI  GFP M   KI_GFP
C001CHQ KI  JAK M   KI_JAK2
C001CHS KI  JAK M   KI_GFP
C001CHT KI  GFP M   KI_JAK2
C001CHU KI  JAK M   KI_JAK2
ADD REPLY
0
Entering edit mode

So no WT_JAK2? That's a problem

ADD REPLY
0
Entering edit mode

I see. So, in your opinion, there is no way to skip this problem ?

ADD REPLY
1
Entering edit mode

You have three groups, not a 2x2. Just do ~sex + group, and use results() with contrast and test="Wald" to make comparisons of pairs of the three groups.

ADD REPLY
0
Entering edit mode

Thank you very much. I will do so. I was making things even more complicated that they were indeed.

ADD REPLY

Login before adding your answer.

Traffic: 1128 users visited in the last hour
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