Paired Difference Analysis

[ninaxhua]

Hello,

I wanted to ask if using the following model design was appropriate for the question I'm trying to answer. I have samples that are paired (BG and AG status). The samples come from 8 individuals and each individual has a BG and AG sample. Additionally, there are 4 treatment groups (A, C, T, E) and a control (W).

I wanted to compare the difference between the two statuses across the treatments and control. I was able to use this model in DESeq2, but my local statistician suggested that I compare the results to a package that uses a different distribution.

DESeq2 model: ~treatment + treatment:status results(dds, contrast=list("TreatmentC.StatusBG","TreatmentA.StatusBG")) Based on this site, I created a similar model, but received this output:

`# Define the metadata categories TB.treatment <- pData(n_inclusive_MRE)$Treatment TB.ID <- pData(n_inclusive_MRE)$id TB.bs <- pData(n_inclusive_MRE)$BrushStatus

Define the normalisation factor (earlier)

TB.norm.factor = normFactors(n_inclusive_MRE) TB.norm.factor = log2(TB.norm.factor/median(TB.norm.factor)+1)

Create the model

TB.mod <- model.matrix(~ TB.treatment + TB.bs + TB.ID + TB.norm.factor) settings <- zigControl(maxit = 10, verbose = TRUE) TB.fit <- fitZig(obj = n_inclusive_MRE, mod = TB.mod, useCSSoffset = FALSE, control = settings)

output

Coefficients not estimable: TB.IDW21AG TB.IDW21BG TB.norm.factor TB.IDC21BG TB.IDE21BG TB.IDT21BG it= 0, nll=NaN, log10(eps+1)=Inf, stillActive=1141 Coefficients not estimable: TB.IDW21AG TB.IDW21BG TB.norm.factor TB.IDC21BG TB.IDE21BG TB.IDT21BG Error in lm.fit(mmZero, qlogis(pi)) : NA/NaN/Inf in 'y' In addition: Warning messages: 1: Partial NA coefficients for 1141 probe(s) 2: Partial NA coefficients for 1123 probe(s)`

Thanks for the help!

[hcorrada]

One way to debug this is to try using limma with the same model matrix. If that fails, then there is an issue with the design. If it doesn’t fail, it may be an issue with the zero inflation part of the model. Let us know.

On Feb 26, 2018, 11:01 AM -0500, ninaxhua notifications@github.com, wrote:

Hello, I wanted to ask if using the following model design was appropriate for the question I'm trying to answer. I have samples that are paired (BG and AG status). The samples come from 8 individuals and each individual has a BG and AG sample. Additionally, there are 4 treatment groups (A, C, T, E) and a control (W).

I wanted to compare the difference between the two statuses across the treatments and control. I was able to use this model in DESeq2, but my local statistician suggested that I compare the results to a package that uses a different distribution. DESeq2 model: ~treatment + treatment:status results(dds, contrast=list("TreatmentC.StatusBG","TreatmentA.StatusBG")) Based on this site, I created a similar model, but received this output: # Define the metadata categories TB.treatment <- pData(n_paired_merged_MRE)$Treatment TB.ID <- pData(n_paired_merged_MRE)$id TB.bs <- pData(n_paired_merged_MRE)$BrushStatus Define the normalisation factor normfactor_MRE = normFactors(n_paired_merged_MRE) normfactor_MRE = log2(normfactor_MRE/median(normfactor_MRE)+1) Create the model TB.mod <- model.matrix(~ TB.treatment + TB.bs + TB.ID + normfactor_MRE) settings <- zigControl(maxit = 10, verbose = TRUE) TB.fit <- fitZig(obj = n_paired_merged_MRE, mod = TB.mod, useCSSoffset = FALSE, control = settings) output Coefficients not estimable: TB.IDW21AG TB.IDW21BG normfactor_MRE TB.IDC21BG TB.IDE21BG TB.IDT21BG it= 0, nll=NaN, log10(eps+1)=Inf, stillActive=1141 Coefficients not estimable: TB.IDW21AG TB.IDW21BG normfactor_MRE TB.IDC21BG TB.IDE21BG TB.IDT21BG Error in lm.fit(mmZero, qlogis(pi)) : NA/NaN/Inf in 'y' In addition: Warning messages: 1: Partial NA coefficients for 1141 probe(s) 2: Partial NA coefficients for 1137 probe(s) Thanks for the help! — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.