Alternative SubsetFun

[JackDettoLoco]

Thank you for your work on SqueezeMeta and for providing the documentation and tools.

I ran SqueezeMeta following the wiki instructions to process metagenomic and metatranscriptomic data together. At this point I have selected the metatranscriptomic samples in R for DESeq2 analysis (I’m not sure if there are better alternatives).

Using subsetFun I filtered on KEGG map 00361 and continued the analysis. This worked as expected. However, I also wanted to explore other pathways present in the same organisms that carry my function of interest. Using subsetFun alone is not sufficient for this, because it removes all functions that are not the one I’m targeting.

To address this, I devised the following procedure. I’d like to confirm whether this approach makes sense:

filtered <- Hadza[["misc"]][["samples"]][grep("^R.*IT", Hadza[["misc"]][["samples"]])] #filter
r_IT<-subsetSamples(SQM = Hadza,samples = filtered) #filtering original sqm
r_IT1 <- subsetFun(SQM = r_IT, fun = "00361")  # filter on KEGG function
r_IT2 <- subsetBins(SQM = r_IT, bins = r_IT1[["bins"]][["table"]] %>% rownames())  # filter by bins carrying that function
r_IT2[["functions"]][["KEGG"]][["abund"]] -> dfdeseq2  # extract table for DESeq2

Thanks in advance for any feedback or suggestions.

[fpusan]

So what is your function of interest exactly? Every function that is present in KEGG map 00361?

[JackDettoLoco]

My goal here is essentially to reduce the noise in the DESeq2 analysis coming from bacteria that are not involved in HCH degradation.

More specifically, the idea is to focus on a sub-community of bacteria that potentially contribute to HCH degradation. To do that, I select bacteria that contain at least one KEGG Ortholog (KO) belonging to KEGG pathway 00361 (the KEGG pathway related to HCH degradation). From that subset of bacteria, I then compare all their functional KOs, including those that are not part of the KEGG 00361 pathway.

Conceptually, this means that the analysis focuses on the functional dynamics of the bacterial community that actually carries the genetic potential to degrade HCH, rather than on the entire microbial community. This allows me, for example, to investigate whether other functions in that sub-community become enriched or persistent, such as transposases or other xenobiotic degradation pathways, even if they are outside the specific KEGG pathway of interest.

An alternative approach would be to run DESeq2 on the entire community and only later filter KOs belonging to the sub-community of interest. However, in that case the abundance of each KO would still be influenced by reads coming from bacteria that do not belong to the HCH-related sub-community, which is exactly the source of noise I am trying to minimize.

Thanks again for the feedback and for taking the time to look into this.

[fpusan]

Your code is doing what you want but is probably not selecting the right data,

What works

1. Your first call to subsetFun will get you all the ORFs containing the substring "00361" in any of their functional annotations.
2. The bin table of that result will include all the bins that contained at least one such ORFs.
3. Using subsetBins on the original object with that list of bins will give you all the ORFs associated to those bins.

As far as I understand it this would do what you described in your previous post.

What is probably not working

subsetFun filters the ORFs based on the content of SQM$orfs$tables, more precisely it looks for regex matches in the following columns by default: 'Gene name', 'KEGG ID', 'KEGGFUN', 'KEGGPATH', 'COG ID', 'COGFUN', 'COGPATH', 'PFAM' (the columns parameter can be used to restrict the search to certain columns).

The right column in your case would be KEGGPATH, but in there we have full pathway names ("Chlorocyclohexane and chlorobenzene degradation" in your case) rather than KEGG pathway ids ("00361").

So I suspect subsetFun(SQM = r_IT, fun = "00361") is instead retrieving ORFs annotated with KEGG ID "K00361" (nasB; nitrite reductase [NAD(P)H] [EC:1.7.1.4]).

I would instead use:

subsetFun(SQM = r_IT, fun = "Chlorocyclohexane and chlorobenzene degradation", columns="KEGGPATH")

Another tricky thing

Selecting bins that contain at least one KEGG Ortholog (KO) belonging to KEGG pathway of interest may give you many false positives, as some reactions are very promiscuous and can participate in a lot of other pathways. As far as I can tell this will not be a big problem for pathway 00361, since the reactions in there seem to be very specific to aromatic compounds, but it's something to keep in mind.