Phyloseq to metagenomeSeq
0
0
Entering edit mode
nabiyogesh • 0
@nabiyogesh-11718
Last seen 5 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 • 681 views
0
Entering edit mode

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

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)

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.

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.

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

0
Entering edit mode

That's right, that should work.

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))

0
Entering edit mode

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

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)

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.