Part 1: Integrating Parkinson’s disease cohorts#
In this notebook we demo the utility of mgnipy in curating cross-study datasets from MGnify for secondary analysis.
Specifically, we integrate the gut microbiome profiles of multiple Parkinson’s disease vs. healthy control cohorts from various MGnify studies, relying on the sample metadata available on MGnify , ENA or BioSamples for the disease status label.
The mgnipy curated abundance dataset with metadata is then further preprocessed in Part 2.
# uncomment if colab
# !pip install mgnipy
Searching for studies using MGnifier#
To start we configure our MGnipy client and access the MGnify API Studies resource.
We will filter our query to studies of the gut microbiome that mention “parkinson”s disease.
We can preview the resulting query urls via .explain()
from mgnipy import MGnipy
# Initialize MGnipy with a cache directory
MG = MGnipy(cache_dir="downloads")
# Search for studies related to Parkinson's disease in the human gut microbiome
pd_studies = MG.studies(
search="parkinson",
biome_lineage="root:Host-associated:Human:Digestive system:Large intestine:Fecal",
)
# Show all of the request urls for the search i.e., the query set
pd_studies.explain()
https://www.ebi.ac.uk/metagenomics/api/v2/studies?biome_lineage=root%3AHost-associated%3AHuman%3ADigestive+system%3ALarge+intestine%3AFecal&search=parkinson&page=1
looks good. we can proceed with actually executing the list query/queries via .get(). To enrich our list of studies with metadata details we can do this in bulk using .enrich_details() or asynchronously via .aenrich_details()
# as http client manager
async with MG:
# populate study list
await pd_studies.aget()
# enrich study list with metadta
await pd_studies.aenrich_details()
# can view as pandas or even save to file if you prefer
study_meta = pd_studies.metadata
# taking a look
study_meta.to_pandas(expand_nested_dicts=True)
Now that we found some studies that match our sesarch criteria, we can take a look at their datasets.
Using MGazine to access the study datasets#
we can access the mgazine of datasets (kinda like a list of available datasets) via .datasets attribute. The study details we retrieved above will also be passed on to the mgazine
# access mgazine
MZ = pd_studies.datasets
# take a look
print(MZ)
Notice in “Nonempty metadata sets” we can see that the study details we collected above are preserved in the mgazine.
The mgnify_studies attribute is a MGnifyMetadata object so contains all the same methods for viewing e.g.:
MZ_SSU.mgnify_studies.to_pandas(expand_nested_dicts=True)... .to_list()... .to_polars()... .records()etc.
Later on in Using MGnetizer to colllect more metadata we will assign additional sets of metadata to .mgnify_runs and .biosamples_metadata which will also convert the lists of records into a MGnifyMetadata object
Filtering the dataset list#
We can filter by the pipeline version and short descriptions of the datasets.
For the ABaCo demo we will use the taxonomic analyses and we will use v4 onwards due to differences in pipeline versions and specifically SILVA databases that were used for the taxonomic analysis
# we can filter by passing as index
V3 = MZ["Taxonomic assignments"]
V4_5 = MZ["Taxonomic assignments SSU"]
V6 = MZ["Summary of SILVA-SSU taxonomies"]
print(V3, V4_5, V6)
Combining dataset lists#
# can add magazines
MZ_SSU = V3 + V4_5 + V6
# keeping latest pipeline version if multiple output files
dedupe_downloads: list[dict] = (
MZ_SSU.downloads_df()
.sort_values(by='pipeline_version', ascending=False)
.drop_duplicates(subset='accession', keep='first')
).to_dict(orient='records')
MZ_SSU.downloads = dedupe_downloads
# print still works
print(MZ_SSU)
Now that we have filtered our list of datasets a bit, let’s actually get and merge the taxonomic datasets.
(Lazy)loading into one taxonomic dataset#
Currently availble in mgnipy are special MGazines for handling taxonomic assignment datasets in the classic taxa x sample/run format TaxaMGazine or in Darwin core-ready format DWCTaxaMGazine
we can also access these from an existing MGazine instance via .taxonomic or .taxonomic_dwc_ready
taxo = MZ_SSU.taxonomic
now load where the datasets will be merged as polars.LazyFrame
# lazyload the mgnify taxanomic assignments datasets
taxo.load()
# calling to_pandas or to_polars will collect the data and return a dataframe
taxo.to_pandas().head()
additionally there is a method .taxonomic_metadata() that parses “taxonomy” into the taxonomic ranks, returning as pandas or polars dataframe which is configured via arg df_engine=. The default is pandas.
also any run and sample metadata relevant to the observations (i.e., by run accessions) can be merged and viewed using .obs_metadata() again as polars or pandas dataframes. As we know metadata() for the observations in the taxonomic mgazine is not available:
# see first 5 sample's metadata
display(taxo.obs_metadata().head())
# recall the "Nonempty metadata set:"
print(taxo)
if we recall from earlier and again in the print statement, we only had .mgnify_studies enriched.
We still need to enrich with run and/or sample metadata, which we will do next :)
Using MGnetizer to collect more metadata#
MGnetizer is designed to retrieve the rich metadata from MGnify for a given list of MGnify accessions/ids.
First we get the ids to pass on
# separate run and assembly ids
runs_ids: list[str] = [x for x in taxo.runs_accessions if not x.startswith('ERZ')]
assembly_ids: list[str] = [x for x in taxo.runs_accessions if x.startswith('ERZ')]
print(f"Runs: {len(runs_ids)}")
print(f"Assemblies: {len(assembly_ids)}")
Runs: 1024
Assemblies: 952
Now we will instantiate the MGnetizers and pass the accessions/ids.
# initialize mgnetizer for runs and assemblies
mnet_run = MG.mgnetizer(resource="run", all_ids=runs_ids)
mnet_acc = MG.mgnetizer(resource="assembly", all_ids=assembly_ids)
# now making the API calls with context manager
async with MG:
await mnet_run.aenrich(limit=None)
await mnet_acc.aenrich(limit=None)
.metadata attribute of MGnetizer returns a MGnifyMetadata object just like with the MGnifiers. the MGnifyMetadata’s can be combined:
# we can combine the metadata from both runs and assemblies into one MGnifyMetadata
run_metadata = mnet_run.metadata + mnet_acc.metadata
if wanting to do some cleaning, do so (which below I do as a pandas df) and then convert back to list of dicts for the MGazine .mgnify_run property.
Now can pass tidied metadata back to our MGazine of taxonomic data:
# now back to TaxaMGazine instance as list of records
taxo.mgnify_runs = df_run.to_dict(orient="records")
# and now
print(taxo)
# also the metadata is updated
taxo.obs_metadata().info()
However, the available run/sample metadata is not consistent for all e.g. sample__sample_title which has lots missing. We can try to get additional sample metadata from the BioSamples database.
Collecting even more sample metadata using BioSampler#
BioSampler is designed to retrieve the rich sample metadata from BioSamples for a list of Run or Sample ENA accessions.
we again start from the mgnipy instance to automatically pass on the configuration
# getting list of sample ids to go to BioSamples
sample_ids = taxo.mgnify_runs.to_pandas()["sample_accession"].to_list()
# init biosampler
bios = MG.biosampler(sample_ids=sample_ids)
# now making the API calls with context manager
async with bios:
await bios.aenrich(limit=None, incl_ena=False)
again we can pass this metadata as list of records to the MGazine instance to .biosamples_metadata for merging:
# assign the biosamples metadata to the taxo object
taxo.biosamples_metadata = bios.metadata.to_list(drop_duplicates=True)
# now can see the additional metadata set
print(taxo)
TaxaMGazine containing:
- MGnify pipeline versions: ['v3', 'v5', 'v6']
- Number of downloads: 8
- Short descriptions: ['Summary of SILVA-SSU taxonomies',
'Taxonomic assignments',
'Taxonomic assignments SSU']
- Nonempty metadata sets: .mgnify_runs, .mgnify_studies, .biosamples_metadata
and now if we were to look at the observations metadata .obs_metadata() there is even more info.
Finding disease label in sample (obs) metadata#
inspecting the metadata we found disease status distrbuted between the following columns, which we normalise into a new label column has_parkinsons_disease with Y and N
samples for which disease status could not be identified from the metadata are then excluded.
Filtered down to 1606 samples with disease status metadata.
PD samples: 975, Non-PD samples: 631
Studies: ['MGYS00005129' 'MGYS00005601' 'MGYS00001650' 'MGYS00006121'
'MGYS00006759' 'MGYS00005755']
We will only include the studies with samples with disease status metadata and will move on with this demonstration.
Also using .obs we can assign our cleaned metadata set which will be given priority over the other sets:
# filter datasets in mgazine to the filtered studies
taxo.downloads = [x for x in taxo.downloads if x['accession'] in filt_studies]
# also add on the observation metadata that we had prepared just before
taxo.obs = df_obs_filt.reset_index().to_dict(orient="records")
print(taxo)
TaxaMGazine containing:
- MGnify pipeline versions: ['v3', 'v5', 'v6']
- Number of downloads: 6
- Short descriptions: ['Summary of SILVA-SSU taxonomies',
'Taxonomic assignments',
'Taxonomic assignments SSU']
- Nonempty metadata sets: .mgnify_runs, .mgnify_studies, .biosamples_metadata, .obs
then the TaxaMGazine can handle the rest and we can export to_anndata() if we want
ad_tax = taxo.to_anndata()
ad_tax
AnnData object with n_obs × n_vars = 1606 × 3247
obs: 'experiment_type', 'instrument_model', 'instrument_platform', 'sample_accession', 'study_accession', 'updated_at', 'run_accession', 'status', 'pipeline_version', 'sample__accession', 'sample__ena_accessions', 'sample__sample_title', 'sample__biome', 'sample__updated_at', 'study__accession', 'study__ena_accessions', 'study__title', 'study__updated_at', 'study__biome.biome_name', 'study__biome.lineage', 'study__metadata.study_name', 'study__metadata.center_name', 'study__metadata.study_title', 'study__metadata.study_accession', 'study__metadata.study_description', 'study__metadata.secondary_study_accession', 'GivenID', 'RunID', 'SRA accession', 'name', 'taxid', 'ENA first public', 'ENA-CHECKLIST', 'External Id', 'INSDC center name', 'INSDC last update', 'INSDC status', 'Submitter Id', 'broad-scale environmental context', 'collection date', 'description', 'environmental medium', 'geographic location (country and/or sea)', 'geographic location (latitude)', 'geographic location (longitude)', 'host age', 'host diet', 'host disease status', 'host family relationship', 'host sex', 'host subject id', 'local environmental context', 'organism', 'project name', 'scientific_name', 'sequencing method', 'title', 'ENA-FIRST-PUBLIC', 'ENA-LAST-UPDATE', 'INSDC first public', 'age', 'body product', 'disease status', 'environment (biome)', 'environment (feature)', 'environment (material)', 'geographic location (countryand/orsea,region)', 'human gut environmental package', 'investigation type', 'medical history performed', 'miscellaneous parameter', 'multiplex identifiers', 'pcr primers', 'sex', 'target gene', 'target subfragment', 'host body product', 'parkinson', 'timepoint', 'gastrointestinal tract disorder', 'sample collection device or method', 'sample storage duration', 'sample storage temperature', 'INSDC secondary accession', 'NCBI submission model', 'NCBI submission package', 'env_broad_scale', 'env_local_scale', 'env_medium', 'geo loc name', 'host', 'host_phenotype', 'isolate', 'status__biosamples_metadata', 'isolation source', 'lat lon', 'collection_date', 'descrip', 'geo_loc_name', 'isolation_source', 'lat_lon', 'Age_at_collection', 'Anti_inflammatory_drugs', 'Antibiotics_current', 'Antibiotics_past_3_months', 'Antihistamines', 'Asthma_or_COPD_med', 'BMI', 'BioSampleModel', 'Birth_control_or_estrogen', 'Blood_pressure_med', 'Blood_thinners', 'Bristol_stool_chart', 'Case_status', 'Celiac_disease', 'Cholesterol_med', 'Co_Q_10', 'Colitis', 'Constipation', 'Crohns_disease', 'Day_of_stool_collection_abdominal_pain', 'Day_of_stool_collection_bloating', 'Day_of_stool_collection_constipation', 'Day_of_stool_collection_diarrhea', 'Day_of_stool_collection_digestion_issue', 'Day_of_stool_collection_excess_gas', 'Depression_anxiety_mood_med', 'Diabetes_med', 'Diarrhea', 'Do_you_drink_alcohol', 'Do_you_drink_caffeinated_beverages', 'Do_you_smoke', 'GI_cancer_past_3_months', 'Gained_10lbs_in_last_year', 'Hispanic_or_Latino', 'How_often_do_you_eat_FRUITS_or_VEGETABLES', 'How_often_do_you_eat_GRAINS', 'How_often_do_you_eat_NUTS', 'How_often_do_you_eat_POULTRY_BEEF_PORK_SEAFOOD_EGGS', 'How_often_do_you_eat_YOGURT', 'IBD', 'IBS', 'INSDC center alias', 'Indigestion_drugs', 'Intestinal_disease', 'Jewish_ancestry', 'Laxatives', 'Loss_10lbs_in_last_year', 'Pain_med', 'Probiotic', 'Race', 'Radiation_Chemo', 'SIBO', 'Sex', 'Sleep_aid', 'Thyroid_med', 'Ulcer_past_3_months', 'broker name', 'collection_method', 'has_parkinsons_disease'
var: 'Superkingdom', 'Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species'
Using MGni.py we could find Parkinson’s disease cohorts and integrate their abundance tables along with sample metadata from MGnify and BioSamples.
In Part 2 we further pre-process the data including cleaning and normalising the taxonomic matrix, followed by ABaCo for batch correction between the studies.
# exporting to h5ad file
ad_tax.obs = ad_tax.obs.astype(str) #workaround for h5ad export issue with mixed types in obs
ad_tax.write_h5ad('pd.h5ad')