Question: Phyloseq to metagenomeSeq
0
gravatar for nabiyogesh
5 months ago by
nabiyogesh0
nabiyogesh0 wrote:

Hi,

I am new to R and trying to understand metagenomeSeq, I will be very much thankful if you could help me to understand this:

I do have three treatment condition (T1, T2, and T3) and want to see how these treatments affect the microbial diversity in rice roots. I am interested in the pairwise comparison between them.

I used below code but I am not sure whether it is right or wrong, I do also have few queries regarding below code.

> root_data <- phyloseq_to_metagenomeSeq(root_M_filter)



DO I need to used normalise count for further processing? what is setting?



> normalizedMatrix <- MRcounts(root_data, norm = TRUE)





> settings = zigControl(maxit = 1, verbose = FALSE)



Is it correct way to make model matrix?



> mod = model.matrix(~root_data$Treatment)



> colnames(mod) = levels(root_data$Treatment)



> res = fitZig(obj = root_data, mod = mod, control = settings)

> res

An object of class "fitZigResults"

Slot "fit":

An object of class "MArrayLM"

$coefficients

                                         T1           T2          T3

71f84e7910006f22684121564206e8ca -2.8129112 -0.321366092  0.06113294

54e204fb99c80764e964456dadd6a0e5  0.5643061  0.019018800  1.01351768

55cd1fc570879d645bbf7a3642e9b0a8 -3.1159264  0.072348076  0.02205983

65f5c31e12c277aec319e2096463f9d2 -0.1355134 -0.004567206  0.53119153

7bed62f0fef250fd831dcf13bf43f4fc  1.9469670 -0.596315540 -0.35124338

                                 scalingFactor

71f84e7910006f22684121564206e8ca    1.62101725

54e204fb99c80764e964456dadd6a0e5   -0.24292482

55cd1fc570879d645bbf7a3642e9b0a8    1.50462644

65f5c31e12c277aec319e2096463f9d2    0.05833637

7bed62f0fef250fd831dcf13bf43f4fc   -0.56247683



What is design here?



> finalMod = slot(res, "fit")$design



It is correct option for contract.make, or should I also involved some more options?



> contrast.matrix = makeContrasts(T1 - T2, levels = finalMod)

> fit2 = contrasts.fit(zigFit, contrast.matrix)

> fit2 = eBayes(fit2)

> topTable(fit2)

                                     logFC   AveExpr         t      P.Value

3147790f0d5a78316fb9dd64f53b9473 11.262064 14.091887 16.331949 3.482815e-15

d2ec9f3b77975c0f457e4b7413b217ff  9.466087 12.212476 10.606122 6.159311e-11

7b908379d1fd2e75d84fd9aaecfd3f77  6.907182  9.784971  9.163780 1.267984e-09

d5fdc1df6dfd31acce9a34b9d8a76590  5.404459  8.373729  7.035926 1.802316e-07

19ddf8e59673c884734cf6093c672da0  5.381863  8.449413  7.011910 1.912363e-07

Kind Regards

limma metagenomeseq • 211 views
ADD COMMENTlink modified 5 months ago • written 5 months ago by nabiyogesh0

Thank you for writing, please excuse the delay in answering:

1) Not knowing how you got to the phyloseq object, I am unable to answer the normalization question. In general, we recommend that you normalize data using cumNorm or (newly added wrenchNorm) functions.

2) You should not set maxit=1 in zigControl as that may not allow the zero-inflated mixture model EM to converge for all features

3) design is the model matrix used for differential abundance. We use limma for the count side of the zero-inflated mixture so you should refer to their documentation for more detail.

Everything else looks ok, Thanks

ADD REPLYlink written 5 months ago by Hector Corrada Bravo60

Thank you very much for your response. I have not done normalization in phyloseq so I am doing it with below commands in metagenomeSeq, here I used maxit =10, could you please suggest to help me to understand the maxit value, how should I know whether I am using correct maxit value or not?

#calculating the normalization factors
#p2 = cumNormStat(obj, pFlag = TRUE, main = "Normalization factors")
obj = cumNorm(obj, p = cumNormStatFast(obj))
normFactor = normFactors(obj)
normFactor = log2(normFactor/median(normFactor) + 1)
settings = zigControl(maxit = 10, verbose = TRUE)
Accession = pData(obj)$Accession
mod = model.matrix(~Accession)
colnames(mod) = levels(Accession)
colnames(mod)
ADD REPLYlink modified 5 months ago • written 5 months ago by nabiyogesh0

You should leave as default. If there was a convergence issue, you will see a warning, in which case you can set to a larger value. In my experience, things converge very quickly for most datasets.

ADD REPLYlink written 5 months ago by Hector Corrada Bravo60

You should leave as default. If there was a convergence issue, you will see a warning, in which case you can set to a larger value. In my experience, things converge very quickly for most datasets.

ADD REPLYlink written 5 months ago by Hector Corrada Bravo60

Thanks Hector, do you mean I should use this command like below:

settings = zigControl()

thanks for your time and help.

nabiyogesh

ADD REPLYlink written 5 months ago by nabiyogesh0

That's right, that should work.

ADD REPLYlink written 5 months ago by Hector Corrada Bravo60

Dear Hector,

Thanks for your response, I also want to compare abundant ASV from metagenomeSeq with normalized ASV count, for that should I use output from below command, will it give me normalized ASV count? and then I want to convert it into log2 value and can cross check the metagenomeSeq result.

obj = cumNorm(obj, p = cumNormStatFast(obj))

ADD REPLYlink written 5 months ago by nabiyogesh0

Use the MRcounts function to get counts, it has arguments to specify normalized and/or logged counts.

ADD REPLYlink written 5 months ago by Hector Corrada Bravo60

thanks Hector for your all help.

obj <- MRcounts(obj, norm = TRUE, log = TRUE).

do you also suggest while running MetagenomeSeq should I also use MRcount for normalisation instead of these three line R code:

obj = cumNorm(obj, p = cumNormStatFast(obj))
normFactor = normFactors(obj)
normFactor = log2(normFactor/median(normFactor) + 1)
ADD REPLYlink written 5 months ago by nabiyogesh0

They are equivalent more or less. The latter gives you more control over how normalization is done. Also, please consider using wrenchNorm instead of cumNorm which we have seen works well.

ADD REPLYlink written 4 months ago by Hector Corrada Bravo60
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: 204 users visited in the last hour