model for LRT
1
0
Entering edit mode
@takashik1201-18730
Last seen 2.8 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:

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

Thank you for your kind helps!

deseq2 LRT • 351 views
ADD COMMENT
1
Entering edit mode
@mikelove
Last seen 26 minutes 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
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:

  1. I think new phenotype data is good for my purpose but I cannot be 100% sure. Could you give some comments for new test?
  2. 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
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.

ADD REPLY
0
Entering edit mode

OK! Thanks!

ADD REPLY

Login before adding your answer.

Traffic: 310 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