Question: model for LRT
0
gravatar for takashi.k1201
3 months ago by
takashi.k12010 wrote:

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 • 127 views
ADD COMMENTlink modified 3 months ago by Michael Love22k • written 3 months ago by takashi.k12010
Answer: model for LRT
1
gravatar for Michael Love
3 months ago by
Michael Love22k
United States
Michael Love22k wrote:

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 COMMENTlink written 3 months ago by Michael Love22k

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 REPLYlink modified 3 months ago • written 3 months ago by takashi.k12010

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 REPLYlink written 3 months ago by Michael Love22k

OK! Thanks!

ADD REPLYlink written 3 months ago by takashi.k12010
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 177 users visited in the last hour