[z_Deprecated] Data generation & collection

1.1.1 Simulation data generation

A comprehensive evaluation requires stringent simulation and corresponding analyses. Here, we employ MIDASim, an intuitive microbiome data simulator that recapitulates microbiome data structure without distributional assumptions while allowing users to manipulate variables of interest, in our case batch effect (caused by multiple studies) and condition (biological) effect. We further edited codes to enable the simulation to include varying odds ratio (between biological signals and batch in each taxa), existing and non-existing relationship between batch effect and library size, as well as the generation of both count and relative abundance data.

To reproduce one iteration of our data using MIDASim, one can do the following to generate minimally viable simulated dataset

• First, in command line:

cd benchmark # enter the benchmark directory
R # to initiate the R environment

• You've now entered R, in which you can now generate the data using the following commands:

source('./generate_data_MIDASim.R') # load generation file
# the odds ratio/condition (biological) effect/batch effect you wish to generate the data for
or_l = c(1.25)
cond_effect_val_l = c(0.5)
batch_effect_val_l = c(0.5)
output_root = './trial_new/simulate' # please make sure the output_root exists
# run functions
scaled_slurm_midas_data_generation(output_root, otu_original, n, or_l, cond_effect_val_l, batch_effect_val_l, iter=GLOBAL_ITER, batch_libsize_related = FALSE, libsize_l=sampled_libsize_l)

• This will help you generate a csv file containing the simulated count matrix, a csv file containing the simulated relative abundance matrix, a csv file containing the simulated meta data, and two txt files containing the simulated batch and condition assignments. Examples are available here. Note that this is just an example that we did not generate with a set random seed - if you wish to generate reproducible outputs. you can do so by adding random seed to the top of the generate_data_MIDASim.R script.

To explain the arguments above:

• otu_original is real-world microbiome data we generate the simulation upon (default is a cleaned 150-sample 301-taxa ibd cohort from the Human Microbiome Project, loaded by ibd_150.Rdata in the same directory).
• n is the number of samples in the simulated dataset (default=450).
• or_l, cond_effect_val_l, batch_effect_val_l: three parameter lists to generate data for.
• iter is simulation iteration for data generation (default=1 without user selection).
• batch_libsize_related, a boolean decides whether the simulated dataset has an established relationship between batch effect and library, if TRUE then we turn the odds ratio into binary correlation between the occurence likehood of batch effect on a taxa and library size, otherwise the relationship is random.
• libsize_l, a list of library_sizes for samples defined, only relevant if batch_libsize_related == TRUE.

This code chunk will generate a count and its corresponding relative abundance datasets for each paratemer combination in the three lists defined above, along with the metadata for samples (which ones in which batch) as well as metadata for taxa (which ones are ground truth perturbed biomarkers) for later experiments.

Note that you may always replace the input parameter list with a complete list so that you generate data for all parameter combinations, which might take a bit longer:

or_l = c(1, 1.25, 1.5)
cond_effect_val_l = c(0, 0.25, 0.5, 0.75, 1)
batch_effect_val_l = c(0, 0.25, 0.5, 0.75, 1)

To run the simulation script in scale, in the folder benchmark/slurm_bash_scripts there is a bash script called step0_run_simulate_sim.sh, which one can revise to generate their own slurm bash scripts based on and move back into the benchmark folder for running the generate_data_MIDASim.R script in scale with slurm. Note that for running MIDASim data generation in slurm, you might need to uncomment the commented lines at the bottom of the generate_data_MIDASim.R script.

1.2 Real world microbiome data collection and cleaning

A comprehensive evaluation also requires well-collected and well-preprocessed real-world microbiome data to ensure that the simulation experiment results work in practice so that in the real-world datasets, biomarker and trend discovery in integrated datasets can be achieved while removing the artifact from combining multiple studies.

To this end, we collect and clean both count and relative abundance real world datasets from two databases, MicrobiomeHD and CuratedMetagenomicsData (CMD). The real world data preparation can be done by running the following code in the command line in the data directory:

• pre-clean the two microbiomeHD datasets, moving them from different batches (studies) to the same files

ipython # enter interactive python mode

• and you can do the following to load the microbiomeHD datasets

from step0_microbiomeHD_precleaning import * 
overall_path = '/athena/linglab/scratch/chf4012'
# autism 2 microbiomeHD
output_dir_path = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/intermediate_autism_2_microbiomeHD/autism_2_microbiomeHD'
address_directory = overall_path+'/mic_bc_benchmark/data/autism_2_microbiomeHD'
data_mat, meta_data = load_data_microbiomeHD(address_directory, output_dir_path)
# cdi 3 microbiomeHD
output_dir_path = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/intermediate_cdi_3_microbiomeHD/cdi_3_microbiomeHD'
address_directory = overall_path+'/mic_bc_benchmark/data/cdi_3_microbiomeHD'
data_mat, meta_data = load_data_microbiomeHD(address_directory, output_dir_path

• then, loads the four datasets (two microbiomeHD datasets from above and two CuratedMetagenomicsData datasets automatically loadable with step1_clean_and_prune.R) into phyloseq objects, clean and prune samples and taxa below 0.05% abundance, and save to a folder

R # enter intractive R mode

then,

source('step1_clean_and_prune.R')
overall_path = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/'
autism_phyloseq_obj <- load_phyloseq_from_merged_microbiomeHD(paste0(overall_path, 'intermediate_autism_2_microbiomeHD/autism_2_microbiomeHD_count_data.csv'), paste0(overall_path, 'intermediate_autism_2_microbiomeHD/autism_2_microbiomeHD_meta_data.csv'))
clean_prune_save_phyloseq(autism_phyloseq_obj, 'autism_2_microbiomeHD', 0.05, 0.05, save = TRUE, save_to = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/pruned_autism_2_microbiomeHD')

cdi_phyloseq_obj <- load_phyloseq_from_merged_microbiomeHD(paste0(overall_path, 'intermediate_cdi_3_microbiomeHD/cdi_3_microbiomeHD_count_data.csv'), paste0(overall_path, 'intermediate_cdi_3_microbiomeHD/cdi_3_microbiomeHD_meta_data.csv'))
clean_prune_save_phyloseq(cdi_phyloseq_obj, 'cdi_2_microbiomeHD', 0.05, 0.05, save = TRUE, save_to = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/pruned_cdi_3_microbiomeHD')

ibd_phyloseq_obj <- load_phyloseq_from_merged_CMD(c("HMP_2019_ibdmdb", "LiJ_2014", "NielsenHB_2014"), c("IBD", 'healthy'))
clean_prune_save_phyloseq(ibd_phyloseq_obj[[1]], ibd_phyloseq_obj[[2]], 0.05, 0.05, save = TRUE, save_to = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/pruned_ibd_3_CMD')

crc_phyloseq_obj <- load_phyloseq_from_merged_CMD(c("FengQ_2015", "HanniganGD_2017", "ThomasAM_2018a", "YachidaS_2019", "ZellerG_2014"), c("adenoma", "healthy"))
clean_prune_save_phyloseq(crc_phyloseq_obj[[1]], crc_phyloseq_obj[[2]], 0.05, 0.05, save = TRUE, save_to = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data/pruned_crc_8_CMD')

• load the phyloseq formatted data from step1 to save to a directory in a format easily digestible by later steps of the pipeline, along with generating a complete confounding check table

ipython

then,

from step2_preprocessing_summarystats import *
overall_path = '/athena/linglab/scratch/chf4012/mic_bc_benchmark/data'
# autism_2_microbiomeHD
data_mat, meta_data = preprocess_data_phyloseq(f'{overall_path}/pruned_autism_2_microbiomeHD', f'{overall_path}/cleaned_data/autism_2_microbiomeHD/autism_2_microbiomeHD', id = 'Sam_id', covar_l = [], relab = False)
data_mat, meta_data = load_results_from_benchmarked_methods(f'{overall_path}/cleaned_data/autism_2_microbiomeHD/autism_2_microbiomeHD_count_data.csv', f'{overall_path}/cleaned_data/autism_2_microbiomeHD/autism_2_microbiomeHD_meta_data.csv')
check_complete_confounding(meta_data, 'Dataset', 'DiseaseState', f'{overall_path}/cleaned_data/autism_2_microbiomeHD/autism_2_microbiomeHD')

# cdi_3_microbiomeHD
data_mat, meta_data = preprocess_data_phyloseq(f'{overall_path}/pruned_cdi_3_microbiomeHD', f'{overall_path}/cleaned_data/cdi_3_microbiomeHD/cdi_3_microbiomeHD', id = 'Sam_id', covar_l = [], relab = False)
data_mat, meta_data = load_results_from_benchmarked_methods(f'{overall_path}/cleaned_data/cdi_3_microbiomeHD/cdi_3_microbiomeHD_count_data.csv', f'{overall_path}/cleaned_data/cdi_3_microbiomeHD/cdi_3_microbiomeHD_meta_data.csv')
check_complete_confounding(meta_data, 'Dataset', 'DiseaseState', f'{overall_path}/cleaned_data/cdi_3_microbiomeHD/cdi_3_microbiomeHD')

# ibd_3_CMD
data_mat, meta_data = preprocess_data_phyloseq(f'{overall_path}/pruned_ibd_3_CMD', f'{overall_path}/cleaned_data/ibd_3_CMD/ibd_3_CMD', id = 'Sam_id', covar_l = ['disease', 'gender', 'age_category'])
data_mat, meta_data = load_results_from_benchmarked_methods(f'{overall_path}/cleaned_data/ibd_3_CMD/ibd_3_CMD_count_data.csv', f'{overall_path}/cleaned_data/ibd_3_CMD/ibd_3_CMD_meta_data.csv')
check_complete_confounding(meta_data, "study_name", "disease", f'{overall_path}/cleaned_data/ibd_3_CMD/ibd_3_CMD')

# crc_8_CMD
data_mat, meta_data = preprocess_data_phyloseq(f'{overall_path}/pruned_crc_8_CMD', f'{overall_path}/cleaned_data/crc_8_CMD/crc_8_CMD', id = 'Sam_id', covar_l = [])
data_mat, meta_data = load_results_from_benchmarked_methods(f'{overall_path}/cleaned_data/crc_8_CMD/crc_8_CMD_count_data.csv', f'{overall_path}/cleaned_data/crc_8_CMD/crc_8_CMD_meta_data.csv')
check_complete_confounding(meta_data, "study_name", "disease", f'{overall_path}/cleaned_data/crc_8_CMD/crc_8_CMD')

These code chunks result in the four subfolders in data/cleaned_data (example here), each containing three files: the count (relative abundance) OTU data, the corresponding metadata, and the complete confounding checking table (number of samples in each biological condition in each cleaned dataset).