Question

Phyloseq to metagenomeSeq

0

Entering edit mode

Jitendra ▴ 10

@nabiyogesh-11718

Last seen 4 months ago

United Kingdom

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 • 2.2k views

ADD COMMENT • link 5.5 years ago Jitendra ▴ 10

0

Entering edit mode

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 REPLY • link 5.5 years ago Hector Corrada Bravo ▴ 60

0

Entering edit mode

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 REPLY • link 5.5 years ago Jitendra ▴ 10

0

Entering edit mode

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 REPLY • link 5.5 years ago Hector Corrada Bravo ▴ 60

0

Entering edit mode

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 REPLY • link 5.5 years ago Hector Corrada Bravo ▴ 60

0

Entering edit mode

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

settings = zigControl()

thanks for your time and help.

nabiyogesh

ADD REPLY • link 5.5 years ago Jitendra ▴ 10

0

Entering edit mode

That's right, that should work.

ADD REPLY • link 5.5 years ago Hector Corrada Bravo ▴ 60

0

Entering edit mode

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 REPLY • link 5.5 years ago Jitendra ▴ 10

0

Entering edit mode

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

ADD REPLY • link 5.5 years ago Hector Corrada Bravo ▴ 60

0

Entering edit mode

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 REPLY • link 5.5 years ago Jitendra ▴ 10

0

Entering edit mode

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 REPLY • link 5.4 years ago Hector Corrada Bravo ▴ 60