ChiLiubio / microeco

An R package for data analysis in microbial community ecology
GNU General Public License v3.0
194 stars 56 forks source link

Error while differential analysis with aldex2 with humann stratified files #383

Open RoyDibakar opened 1 month ago

RoyDibakar commented 1 month ago

Hello, I am trying to analyze the humann output files to obtain the differentially abundant pathways in each group using aldex2_kw.

however I got this error Error in aldex.clr.function(reads, conds, mc.samples, denom, verbose, : not all reads are integers

Here, i used the example files that are given in the tutorial.

Thank you.

ChiLiubio commented 1 month ago

Hi. Yes. The aldex2 method only accepts the integers like the so-called counts from amplicon sequencing. You can use round function convert the data to intergers if you want to use this method. It is also feasible to use other methods, as humann output is RPK (reads per kilobase) normalized data.

your_microtable$otu_table %<>% round

Best, Chi

RoyDibakar commented 1 month ago

Thank you for your answer. However, now i got this error

Error in aldex.clr.function(reads, conds, mc.samples, denom, verbose, : one or more reads are negative

while computing the differential abundance of taxa column ("Species") but it work perfectly when i took the "pathway" column.

ChiLiubio commented 1 month ago

Could you please attach your scripts? I will reproduce the issue to check how it happened.

RoyDibakar commented 1 month ago

Hello, I am attaching the codes as well as, a part of my input data. I want to use aldex2_kw to find the differential microbial pathways.

Thank you. The codes are- sample_file_path <- "path/Sample.txt" match_file_path <- "path/Match.txt" abund_file_path <- "path/Pathway.txt"

humann2meco(abund_file_path, db = "MetaCyc") humann2meco(abund_file_path, db = "MetaCyc", sample_table = sample_file_path, match_table = match_file_path)

testy <- humann2meco(abund_file_path, db = "MetaCyc", sample_table = sample_file_path, match_table = match_file_path) testy$tidy_dataset()

taxa biomarker

test$cal_abund(select_cols = 4:10, rel = TRUE) test$taxa_abund$Phylum %<>% .[!grepl("unclass", rownames(.)), ]

ALDEx2_kw

test1 <- trans_diff$new(testy, method = "ALDEx2_kw", group = "Group", taxa_level = "Species", p_adjust_method = "none")

Error in aldex.clr.function(reads, conds, mc.samples, denom, verbose, : one or more reads are negative

Match.txt pathway.txt Sample.txt

ChiLiubio commented 1 month ago

Hi. The reason is ALDEx2_kw method invokes the data in feature abundance table, i.e. otu_table, instead of relative abundance in taxa_abund list. So we need to generate a prepared feature abundance table. Here is the modified code.

library(microeco)
library(file2meco)
library(magrittr)

sample_file_path <- "path/Sample.txt"
match_file_path <- "path/Match.txt"
abund_file_path <- "path/Pathway.txt"

test <- humann2meco(abund_file_path, db = "MetaCyc", sample_table = sample_file_path, match_table = match_file_path)
test$tidy_dataset()

# clone microtable object to avoid modify raw object
tmp_test <- clone(test)
# only select the taxonomy
tmp_test$tax_table %<>% .[, 4:10]
# merge otu_table to generate a new microtable object
tmp_species <- tmp_test$merge_taxa("Species")
# optimize the rownames; use speices names
rownames(tmp_species$otu_table) <- rownames(tmp_species$tax_table) <- tmp_species$tax_table$Species 
tmp_species$otu_table %<>% .[rownames(.) != "unclassified", ]
# use integers
tmp_species$otu_table %<>% round
tmp_species$tidy_dataset()

#ALDEx2_kw
test1 <- trans_diff$new(tmp_species, method = "ALDEx2_kw", group = "Group", taxa_level = "Species", p_adjust_method = "none")
RoyDibakar commented 1 month ago

Hello, Thank you very much for your reply. It works perfectly now. But, I have one more query.

What if I merge pathway and species column in the tax table, and the try to find (using ALDEx2_kw) the differential pathways and their contributing taxa, will it work? It works with Linda, lefse algorithm.

Thank you.

ChiLiubio commented 1 month ago

I have not tried to use ALDEx2 to analyze such data. Maybe it will work. Please have a try.

RoyDibakar commented 1 month ago

I tried actually. And it gave the same error.

Error in aldex.clr.function(reads, conds, mc.samples, denom, verbose, : one or more reads are negative.

Is there a way to solve this problem?

Thank you.

ChiLiubio commented 1 month ago

Please also attach the steps so that I can reproduce it?

RoyDibakar commented 1 month ago

Okay. Here are the steps.

library(microeco) library(file2meco) library(magrittr)

sample_file_path <- "path/Sample.txt" match_file_path <- "path/Match.txt" abund_file_path <- "path/Pathway.txt"

test <- humann2meco(abund_file_path, db = "MetaCyc", sample_table = sample_file_path, match_table = match_file_path) test$tidy_dataset()

make a dataframe with a new column

a = as.data.frame(test $taxtable) a$new = paste(a$pathway, a$Species, sep="") colnames(a)[11] <- "Pathway_Species"

add the tax table back to microtable object

meco_object <- microtable$new(sample_table = test$sample_table, otu_table =test$otu_table, tax_table = a) meco_object$tidy_dataset()

clone microtable object to avoid modify raw object

tmp_test <- clone(meco_object )

only select the taxonomy

tmp_test$tax_table %<>% .[, 4:11]

merge otu_table to generate a new microtable object

tmp_species <- tmp_test$merge_taxa("Pathway_Species")

optimize the rownames; use speices names

rownames(tmp_species$otu_table) <- rownames(tmp_species$tax_table) <- tmp_species$tax_table$Pathway_Species tmp_species$otu_table %<>% .[rownames(.) != "unclassified", ]

use integers

tmp_species$otu_table %<>% round tmp_species$tidy_dataset()

ALDEx2_kw

test1 <- trans_diff$new(tmp_species, method = "ALDEx2_kw", group = "Group", taxa_level = "Pathway_Species", p_adjust_method = "none")

Thank you.

ChiLiubio commented 1 month ago

Yes. It is a problem as you say. Please switch to another method.

RoyDibakar commented 1 month ago

Okay. In Ancombc method, can we plot only top 25 pathways? The use_number() function is not working for ancombc.

Thank you.

ChiLiubio commented 1 month ago

As you show, the method currently donot support this. I will try to update it. Thanks.

RoyDibakar commented 1 month ago

Yes please. It will be a great help.

Thank you