<jupyter_start><jupyter_text>Cross-Study analysis, using MGnifyRThe [MGnify API](https://www.ebi.ac.uk/metagenomics/api/v1) returns data and relationships as JSON. [MGnifyR](https://github.com/beadyallen/MGnifyR) is a package to help you read MGnify data into your R analyses.**This example shows you how to perform an analysis across multiple Studies. It is an example of comparing taxonomic diversity from two places in different Studies.**You can also discover more about the "API" using the [Browsable API interface in your web browser](https://www.ebi.ac.uk/metagenomics/api/v1).This is an interactive code notebook (a Jupyter Notebook).To run this code, click into each cell and press the ▶ button in the top toolbar, or press `shift+enter`.---<jupyter_code>library(IRdisplay)
display_markdown(file = '../_resources/mgnifyr_help.md')<jupyter_output><empty_output><jupyter_text>Load libraries:<jupyter_code>library(vegan)
library(ggplot2)
library(phyloseq)
library(MGnifyR)

mg <- mgnify_client(usecache = T, cache_dir = '/tmp/mgnify_cache')<jupyter_output><empty_output><jupyter_text>Example: compare taxonomic abundances of two soil studies*This example is based on a [MGnify workshop exercise](https://beadyallen.github.io/MGnifyR/Exercises.html) created by [<NAME>](https://github.com/beadyallen) (the author of MGnifyR).*  Fetch the MGnify Analyses accessions for each of two Studies(one with samples from Malaysia, one with samples from Panama)<jupyter_code>panama <- mgnify_analyses_from_studies(mg, 'MGYS00003920')
malaysia <- mgnify_analyses_from_studies(mg, 'MGYS00003918')<jupyter_output><empty_output><jupyter_text>Join the Analyses accession lists<jupyter_code>accessions <- c(panama, malaysia)
sprintf('There are %d accessions between the studies', length(accessions))<jupyter_output><empty_output><jupyter_text>Fetch metadata for all of the Analyses from the MGnify API<jupyter_code>metadata <- mgnify_get_analyses_metadata(mg, accessions)
head(metadata)<jupyter_output><empty_output><jupyter_text>Taxonomic analysisFirst, build a [phyloseq](https://joey711.github.io/phyloseq/) object of the Analyses metadata<jupyter_code>ps <- mgnify_get_analyses_phyloseq(mg, metadata$analysis_accession)<jupyter_output><empty_output><jupyter_text>Filter out low-abundance samplesFilter out samples with low abundances, leaving just those in the "normal" distribution. This is particularly important given the normalisation approach we're about to use ([`phyloseq`'s slightly controversial `rarefy_even_depth`](https://www.rdocumentation.org/packages/phyloseq/versions/1.16.2/topics/rarefy_even_depth).)We make a histogram and note that samples with abundances $\lt 10^3$ i.e. $\lt 1000$ seem to be outliers.<jupyter_code>hist(log10(sample_sums(ps)), breaks = 50)
ps_good <- subset_samples(ps, sample_sums(ps) > 1000)
hist(log10(sample_sums(ps_good)), breaks = 50)<jupyter_output><empty_output><jupyter_text>Estimate richnessUse `estimate_richness` to calculate various diversity measures for each analysis.To read the documentation on this method, enter `?estimate_richness` in a Code Cell and run it.<jupyter_code>#  ?estimate_richness
diversity = estimate_richness(ps_good)
head(diversity)

diversity[diversity$Observed == max(diversity$Observed),]<jupyter_output><empty_output><jupyter_text>Normalise the data by rarefication<jupyter_code>ps_rare <- rarefy_even_depth(ps_good)
div_rare <- estimate_richness(ps_rare)
merged_df <- merge(div_rare, metadata, by = 0, all.y = F)
head(merged_df)
options(repr.plot.width=12, repr.plot.height=8)
ggplot(merged_df, aes(x=`sample_geo-loc-name`, y=Observed)) + geom_boxplot() + theme(text = element_text(size = 20))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1060893098782138, 'fraction_numerical': 0.013802435723951286, 'mean_word_length': 7.676056338028169, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 39, 'https://': 7, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '414897', 'license': 'Apache-2.0', 'n_tokens_mistral': 1311, 'n_tokens_neox': 1203, 'n_words': 302, 'path': 'notebooks-src/notebooks/R Examples/Cross-study analysis.ipynb', 'repo_name': 'EBI-Metagenomics/notebooks', 'size': '3698'}