#-----
# For the FDR adjustment in MaAslin2
#-----



## load packages
rm(list=ls())
require(microbiome)

require(tidyverse)
#require(curatedomcis)
require(curatedMetagenomicData)
#BiocManager::install("curatedMetagenomicData")
require(clusterProfiler)
require(readxl)
require(data.table)
require(Maaslin2)

## rfor the test of maaslin2 
input_data = system.file("extdata", "HMP2_taxonomy.tsv", package="Maaslin2") # The abundance table file
input_data

input_metadata = system.file("extdata", "HMP2_metadata.tsv", package="Maaslin2") # The metadata table file
input_metadata

#get the pathway (functional) data - place holder
download.file("https://raw.githubusercontent.com/biobakery/biobakery_workflows/master/examples/tutorial/stats_vis/input/pathabundance_relab.tsv", "./pathabundance_relab.tsv")

## then we test the file as data.frame

df_input_data = read.table(file = input_data, header = TRUE, sep = "\t", row.names = 1, stringsAsFactors = FALSE)
df_input_data[1:5, 1:5]

df_input_metadata = read.table(file = input_metadata, header = TRUE, sep = "\t", row.names = 1, stringsAsFactors = FALSE)
df_input_metadata[1:5, ]

df_input_path = data.frame(read.csv("./pathabundance_relab.tsv", sep = "\t", stringsAsFactors = F, row.names = 1))
df_input_path[1:5, 1:5]

##
fit_data2 = Maaslin2(
  input_data = df_input_data, 
  input_metadata = df_input_metadata, 
  output = "demo_output2", 
  fixed_effects = c("diagnosis", "dysbiosis"),
  reference = c("diagnosis,nonIBD"))

all_res<-read_tsv("C:\\Users\\kaluo\\Documents\\demo_output2\\all_results.tsv")

## here we focus on dysbiosis 

result_dysbiosis<-all_res%>%filter(metadata=="dysbiosis")

result_dysbiosis$fdr_p_only_dysbiosis<-p.adjust(result_dysbiosis$pval,method = "BH")

# for the all_res

all_res$qval_lk<-p.adjust(all_res$pval,method = "BH")
all_res%>%as.data.table()%>%print() 
# note that qval_lk and qval is equal, suggesting that the post-hoc fdr adjustment is for all variables (including the main interest one,like dysbiosis here)


result_dysbiosis%>%as.data.table()%>%print() 
# note that qval and fdr_p_only_dysbiosis is not equal, indicating that the poct-hoc fdr adjustment for one specific variable is not equal to the default qvalues in MaAslin2


####
# Question
####

#**************
# As shown above, is it somewhat overadjustment in MaAslin2 when it comes to FDR correction?
# Usually, we are only interested in one specific variable, so would it be reasonable to perform
# the FDR adjustment only for the selected variable?
#**************