Question

For comparimg 4 or more tissue types using DRIMSeq

0

Entering edit mode

rohitsatyam102 ▴ 20

@rohitsatyam102-24390

Last seen 6 weeks ago

India

Hi Everyone!!

I am analyzing 4 tissue types obtained from the 4 different layers of a plant, sampled at 4 different time points ( days after flowering (daf). Each tissue has 2-3 biological replicate. The samples are grouped as follows

sample_id	group
cellArep116daf A16daf
cellArep216daf	A16daf
cellArep316daf	A16daf
cellBrep116daf	B16daf
cellBrep216daf	B16daf
cellBrep14daf	B4daf
cellBrep24daf	B4daf
cellBrep18daf	B8daf
cellNrep116daf	N16daf
cellNrep216daf	N16daf
cellNrep316daf	N16daf
cellNrep14daf	N4daf
cellNrep24daf	N4daf
cellNrep34daf	N4daf
cellNrep18daf	N8daf
cellNrep28daf	N8daf
cellNrep38daf	N8daf
cellOrep116daf	O16daf
cellOrep14daf	O4daf
cellOrep24daf	O4daf
cellOrep34daf	O4daf
cellOrep18daf	O8daf
cellOrep28daf	O8daf

I wish to perform differential transcript usage analysis to understand how they differ between different time points for same cell type (say differences between B4daf vs B8daf and so on) and then for same time point but different cell types (B16daf vs A16daf). I am using DRIMSeq for that purpose.

I came across this post wherein the section Create the design matrix without intercept due to no clear 'ground state is somewhat unclear to me (the phrase itself is unclear to me). For a scenario like mine where there is no control case condition, how removing the intercept from the design matrix will affect my analysis, and is it desirable to do so because here it is mentioned that we should never remove the intercept term.

I am confused as to how should I go about it. Also setting one_way = FALSE would be appropriate in my case and how should I set contrast if I am comparing 4 tissues at a time harvested on the same day (say 16daf i.e A16daf vs B16daf vs N16daf vs O16daf)?

#DRIMSeq #DTU #contrast • 1.1k views

ADD COMMENT • link updated 4.1 years ago by Gosia Nowicka ▴ 90 • written 4.1 years ago by rohitsatyam102 ▴ 20

0

Entering edit mode

So I myself tried comparing the results obtained with/without intercept on my data and I found more genes with differential transcript usage while using intercept. Also, I observed more genes having significant p-values when using intercept. I used the following code:

    set.seed(1)
    d_with_intercept <- dmPrecision(d, design=design_full)
    ###design_null <- model.matrix(~ 1, data=DRIMSeq::samples(d))
    d_with_intercept <- dmFit(d_with_intercept, design = design_full, verbose = 1)

    my.contrasts_with_intercept <- makeContrasts(
        A.vs.B.16daf = Intercept-B16daf,
        A.vs.N.16daf = Intercept-E16daf,
        levels=design_full)
    d_with_intercept <- dmTest(d_with_intercept, contrast = my.contrasts_with_intercept[,"A.vs.N.16daf"], verbose = 1)
    plotPValues(d_with_intercept, level = "feature")
    plotPValues(d_with_intercept, level = "gene")

And without intercept

design_without_intercept <- model.matrix(~0+group, data=DRIMSeq::samples(d))
design_without_intercept
colnames(design_without_intercept) <- gsub('group', '', colnames(design_without_intercept))
colnames(design_without_intercept)
set.seed(1)
d_without_intercept <- dmPrecision(d, design=design_without_intercept)
###design_null <- model.matrix(~ 1, data=DRIMSeq::samples(d))
d_without_intercept <- dmFit(d_without_intercept, design = design_without_intercept, verbose = 1)

my.contrasts_without_intercept <- makeContrasts(
    A.vs.B.16daf = A16daf-B16daf,
    A.vs.N.16daf = A16daf-N16daf,
    levels=design_without_intercept)

d_without_intercept <- dmTest(d_without_intercept, contrast = my.contrasts_without_intercept[,"A.vs.N.16daf"], verbose = 1)

Will it be worthwhile to go with the intercept approach simply because it gives us more genes with differential transcript usage? In case yes, how do I explain this discrepancy?

enter image description here

ADD REPLY • link 4.1 years ago rohitsatyam102 ▴ 20

score 4 · Accepted Answer · 2020-10-20

Hello,

the set up of the design matrix (with/without intercept) and contrast depends on the question in hand. Based on my experience, I follow this strategy. Whenever I have a study with a reference group and then one or multiple other groups that I want to compare against that reference group, I use design with intercept and I test the remaining coefficients whether they are equal to zero. In case I want to do different pair-wise comparisons I find easier to interpret the design without intercept and setting up contrasts for each pair-wise comparison.

To illustrate the difference between the two designs with and without intercept I have prepared some plots for a case with 3 groups :D They show you the interpretation of the coefficients that are estimated with those two designs.

Just as a reminder the reference level is set up to the first level in your group variable which should be a factor. Let's assume it is A16daf.

With the design with intercept design_full <- model.matrix(~ 1 + group, data=DRIMSeq::samples(d)) if you want to test for the difference between B16daf and A16daf, you need to test only the B16daf coefficient.

my.contrasts_with_intercept <- makeContrasts(
    A.vs.B.16daf = B16daf,
    A.vs.N.16daf = E16daf,
    levels=design_full)

The comparison where you have a contrast Intercept-B16daf is not quite correct for what you want to ask... You compare the size of the blue and orange rectangle, which is asking if the transcript usage in group B is twice the usage in group A.

With the design without intercept design_full <- model.matrix(~ -1 + group, data=DRIMSeq::samples(d)) if you want to test for the difference between B16daf and A16daf, you need to set a contrast:

my.contrasts_without_intercept <- makeContrasts(
    A.vs.B.16daf = A16daf - B16daf,
    A.vs.N.16daf = A16daf - E16daf,
    levels=design_full)

You can define as many contrasts as you want for any pair-wise comparison you are interested in.

Let me know if you have any further questions.

sample_id	group	libray_layout
cellArep116daf A	16_daf
cellArep216daf	A	16daf
cellArep316daf	A	16_daf
cellBrep116daf	B	16_daf
cellBrep216daf	B	16_daf
cellBrep14daf	B	4_daf
cellBrep24daf	B	4_daf
cellBrep18daf	B	8_daf
cellNrep116daf	N	16_daf
cellNrep216daf	N	16_daf
cellNrep316daf	N	16_daf
cellNrep14daf	N	4_daf
cellNrep24daf	N	4_daf
cellNrep34daf	N	4_daf
cellNrep18daf	N	8_daf
cellNrep28daf	N	8_daf
cellNrep38daf	N	8_daf
cellOrep116daf	O	16_daf
cellOrep14daf	O	4_daf
cellOrep24daf	O	4_daf
cellOrep34daf	O	4_daf
cellOrep18daf	O	8_daf
cellOrep28daf	O	8_daf