added python-based beta diversity analysis

Author: KunDHuang

docs/beta_diversity_analysis.md

@@ -128,4 +128,69 @@ Below, we are showcasing how to inspect the beta diversity of microbiomes from t
 
 ## Python-based method
 
+#### R packages required
+* [scikit-bio >= 0.5.6](https://scikit.bio/)
+* [pandas >= 1.3.5](https://pandas.pydata.org/)
+* [numpy >= 1.23.5](https://numpy.org/)
+* [matplotlib >= 3.5.0](https://matplotlib.org/)
+* [seaborn >= 0.11.2](https://seaborn.pydata.org/)
+
+#### Beta diversity analysis with PCoA plotting integrating maximum three variables
+Here, we introduce a python script `multi_variable_pcoa_plot.py` in the `path_to_the_package/KunDH-2023-CRM-MSM_metagenomics/scripts` to perform PCoA analysis:
+
+```{python}
+usage: multi_variable_pcoa_plot.py [-h] [--abundance_table [ABUNDANCE_TABLE]] [--metadata [METADATA]] [--transformation [TRANSFORMATION]] [--metric [METRIC]] [--amplifier [AMPLIFIER]] [--sample_column [SAMPLE_COLUMN]] [--variable1 [VARIABLE1]] [--variable2 [VARIABLE2]] [--variable3 [VARIABLE3]]
+                                   [--marker_palette [MARKER_PALETTE]] [--marker_shapes [MARKER_SHAPES]] [--marker_sizes [MARKER_SIZES]] [--output_figure [OUTPUT_FIGURE]] [--test [TEST]] [--df_opt [DF_OPT]] [--font_style [FONT_STYLE]] [--font_size [FONT_SIZE]]
+
+This program is to do PCoA analysis on microbial taxonomic or functional abundance data integrating maximum three variables together.
+
+optional arguments:
+  -h, --help            show this help message and exit
+  --abundance_table [ABUNDANCE_TABLE]
+                        Input the merged abundance table generated by MetaPhlAn.
+  --metadata [METADATA]
+                        Input a tab-delimited metadata file.
+  --transformation [TRANSFORMATION]
+                        Specify the tranformation function applied on data points in the original table. For abundance table, you can choose <sqrt>/<log>. Default setting is <None>.
+  --metric [METRIC]     Specify the metric you want to use for calculating beta diversity in the case of as input using abundance table.<braycurtis>/<unweighted_unifrac>/<jaccard>/<weighted_unifrac>. Default setting is <braycurtis>
+  --amplifier [AMPLIFIER]
+                        Specify how much you want to amplify your original data point. For example, <--amplifier 100> indicates that all original data point times 100. Default is 1.
+  --sample_column [SAMPLE_COLUMN]
+                        Specify the header of column containing metagenome sample names in the metadata file.
+  --variable1 [VARIABLE1]
+                        Specify the header of the variable in the metadata table you want to assess. This variable will be represented by colors.
+  --variable2 [VARIABLE2]
+                        Specify the header of second variable in the metadata table you want to assess. This variable will be represented by marker shapes.
+  --variable3 [VARIABLE3]
+                        Specify the header of the third variable in the metadata table you want to assess. This variable will be represented by marker sizes.
+  --marker_palette [MARKER_PALETTE]
+                        Input a tab-delimited mapping file where 1st column contains group names and 2nd column contains color codes. default: [None] (automatic handling)
+  --marker_shapes [MARKER_SHAPES]
+                        Input a tab-delimited mapping file where 1st column contains group names and 2nd column contains marker shapes. default: [None] (automatic handling)
+  --marker_sizes [MARKER_SIZES]
+                        Input a tab-delimited mapping file where values are group names and keys are marker size. default: [None] (automatic handling)
+  --output_figure [OUTPUT_FIGURE]
+                        Specify the name for the output figure. For example, output_figure.svg
+  --test [TEST]         Specify an output file for saving permanova test results. For example, project_name
+  --df_opt [DF_OPT]     Specify the output name for saving coordinates (PC1 and PC2) for each sample. For example, project_name_coordinates.tsv
+  --font_style [FONT_STYLE]
+                        Specify the font style which is composed by font family and font type, delimited with a comma. default: [sans-serif,Arial]
+  --font_size [FONT_SIZE]
+                        Specify the font size. default: [11]
+
+examples:
+pcoa_painter.py --abundance_table <merged_metaphlan_table> --metadata <metadata> --sample_column <sample_header> --variable1 <variable1_name> --variable2 <variable2_name> --variable3 <variable3_name> --output_figure <output.png>
+```
+
+To demostrate the usage of `multi_variable_pcoa_plot.py`, we will drawa PCoA plot based on [microbiome compistion of samples](../example_data/mvpp_mpa_species_relab.tsv.bz2) from 11 populations grouped as *W (Westernization)*, *NW (Non-Westernization)*, *NWU (Non-Westernization(Urban))* and *MSM (Men-having-sex-with-men)*. Different populations will be assigned with custom colors using a [color map file](../example_data/mvpp_color_map.tsv) and *MSM* population will be highlighted with larger marker size using a [marker size map file](../example_data/mvpp_marker_size_map.tsv). The metadata of each sample is provided by a [metadata file](../example_data/mvpp_metadata.tsv).
+
+Example command:
+~~~
+$multi_variable_pcoa_plot.py --abundance_table mvpp_mpa_species_relab.tsv --metadata mvpp_metadata.tsv --sample_column sample --variable1 country --variable2 westernization --variable3 country --output_figure mvpp_pcoa.png --test mvpp_permanova.tsv --df_opt mvpp_coordinates_df.tsv --marker_palette mvpp_color_map.tsv --marker_sizes mvpp_marker_size_map.tsv
+~~~
+
+![Multiple variable PCoA plot](../images/mvpp_pcoa.png)
+
+As optional ouputs, `multi_variable_pcoa_plot.py` also generates non-adjustment PERMANOVA test (e.g. [mvpp_permanova.tsv](../example_data/mvpp_permanova.tsv)) and coordinates of PC1 and PC2 (e.g. [mvpp_coordinates.tsv](../example_data/mvpp_coordinates.tsv)) which can be used in visualization in other ways we will discuss shortly below.
+
 ## A method mixing R and Python

example_data/mvpp_color_map.tsv

@@ -0,0 +1,11 @@
+SWE	#544179
+ITA	#6166B3
+USA	#32C1CD
+JPN	#17D7A0
+FJI	#E02401
+MDG	#F78812
+TZA	#AB6D23
+GHA	#FFCCD2
+ETH	#FF5C58
+MSM	#888888
+CHN	#000000