Question

model for LRT

0

Entering edit mode

takashi.k1201 • 0

@takashik1201-18730

Last seen 5.6 years ago

Hi,

I have 12 RNA-Seq data sequenced at a time and another 3 data sequenced later. All of them have sex and treatment/untreatment information. I designed phenotype data as follows.

"ids","sex","treatment","seq","group"
"ind01","female","untreated","seq1","untreatedFemale"
"ind02","female","untreated","seq1","untreatedFemale"
"ind03","female","untreated","seq1","untreatedFemale"
"ind04","male","untreated","seq1","untreatedMale"
"ind05","male","untreated","seq1","untreatedMale"
"ind06","male","untreated","seq1","untreatedMale"
"ind07","female","treated","seq1","treatedFemale"
"ind08","female","treated","seq1","treatedFemale"
"ind09","female","treated","seq1","treatedFemale"
"ind10","male","treated","seq1","treatedMale"
"ind11","male","treated","seq1","treatedMale"
"ind12","male","treated","seq1","treatedMale"
"ind13","female","treated","seq2","treatedFemale"
"ind14","female","treated","seq2","treatedFemale"
"ind15","female","treated","seq2","treatedFemale"

I could see DEGs in untreated male VS treated male, untreated female VS treated female, and treated male VS treated female comparison by Wald test. For further biological implication, I would like to see genes which are specifically differentially expressed in one treated male/female group but not in the other treated sex.

I tried two full and reduced models. One is

full=~ sex + treatment + seq + group, reduced =~ sex + treatment + seq

But this was not full rank. Next, I tried

full=~ sex + treatment + seq + sex:treatment, reduced =~ sex + treatment + seq

and it could be analyzed.

My questions are:

What is the differences between these two models?
Is second model appropriate for my purpose?

Thank you for your kind helps!

deseq2 LRT • 948 views

ADD COMMENT • link updated 5.6 years ago by Michael Love 42k • written 5.6 years ago by takashi.k1201 • 0

score 1 · Answer 1 · 2018-12-12

1

Entering edit mode

Michael Love 42k

@mikelove

Last seen 3 hours ago

United States

Group, being a combination of sex and treatment cannot appear in the model with those two. When people add a variable like group they are removing the other two variables, as group now takes the place of those two variables and their interaction.

It sounds like, from the comparisons you are interested in, the best design would be ~batch + group. If you want to find genes which are DE in one treatment group but not the other, you will have to build both gene sets and then perform set operations, there is not a contrast which will provide this gene list.

ADD COMMENT • link 5.6 years ago Michael Love 42k

0

Entering edit mode

Thank you for your valuable comments! I now understand what I want cannot be tested. So, I changed phenotype data from what I showed to

"ids","group","seq"
"ind01","untreated","seq1"
"ind02","untreated","seq1"
"ind03","untreated","seq1"
"ind04","untreated","seq1"
"ind05","untreated","seq1"
"ind06","untreated","seq1"
"ind07","treatedFemale","seq1"
"ind08","treatedFemale","seq1"
"ind09","treatedFemale","seq1"
"ind10","treatedMale","seq1"
"ind11","treatedMale","seq1"
"ind12","treatedMale","seq1"
"ind13","treatedFemale","seq2"
"ind14","treatedFemale","seq2"
"ind15","treatedFemale","seq2"

because I didn't see expression difference between untreated male and female both in clustering (PCA, Heatmap, fuzzy c-means) and two-group comparison (Wald test).

By using the revised phenotype data, I ran LRT with full=~ seq + group, reduced =~ seq and got some results.

My another question is:

I think new phenotype data is good for my purpose but I cannot be 100% sure. Could you give some comments for new test?
I think I should categorize genes in res object obtained by results(dds, filterFun=ihw, alpha=0.05) which genes up/down-regulated in which comparison. Is there good way to do that?

Thank you for your another helps!

ADD REPLY • link 5.6 years ago takashi.k1201 • 0

0

Entering edit mode

The support site is more for specific questions about Bioconductor software, not for guidance on your analysis per se. You can reach out to a bioinformatician or post to another forum such as Biostars.