changeset 0:59a7f3f83aec draft

planemo upload for repository https://github.com/ebi-gene-expression-group/container-galaxy-sc-tertiary/ commit 20f4a739092bd05106d5de170523ad61d66e41fc

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/decoupler_pseudobulk.py	Sun Sep 24 08:44:24 2023 +0000
@@ -0,0 +1,365 @@
+import argparse
+
+import anndata
+import decoupler
+import pandas as pd
+
+
+def get_pseudobulk(
+    adata,
+    sample_col,
+    groups_col,
+    layer=None,
+    mode="sum",
+    min_cells=10,
+    min_counts=1000,
+    use_raw=False,
+):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> adata.X = abs(adata.X).astype(int)
+    >>> pseudobulk = get_pseudobulk(adata, "bulk_labels", "louvain")
+    """
+
+    return decoupler.get_pseudobulk(
+        adata,
+        sample_col=sample_col,
+        groups_col=groups_col,
+        layer=layer,
+        mode=mode,
+        use_raw=use_raw,
+        min_cells=min_cells,
+        min_counts=min_counts,
+    )
+
+
+def prepend_c_to_index(index_value):
+    if index_value and index_value[0].isdigit():
+        return "C" + index_value
+    return index_value
+
+
+# write results for loading into DESeq2
+def write_DESeq2_inputs(pdata, layer=None, output_dir="", factor_fields=None):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> adata.X = abs(adata.X).astype(int)
+    >>> pseudobulk = get_pseudobulk(adata, "bulk_labels", "louvain")
+    >>> write_DESeq2_inputs(pseudobulk)
+    """
+    # add / to output_dir if is not empty or if it doesn't end with /
+    if output_dir != "" and not output_dir.endswith("/"):
+        output_dir = output_dir + "/"
+    obs_for_deseq = pdata.obs.copy()
+    # replace any index starting with digits to start with C instead.
+    obs_for_deseq.rename(index=prepend_c_to_index, inplace=True)
+    # avoid dash that is read as point on R colnames.
+    obs_for_deseq.index = obs_for_deseq.index.str.replace("-", "_")
+    obs_for_deseq.index = obs_for_deseq.index.str.replace(" ", "_")
+    col_metadata_file = f"{output_dir}col_metadata.csv"
+    # write obs to a col_metadata file
+    if factor_fields:
+        # only output the index plus the columns in factor_fields in that order
+        obs_for_deseq[factor_fields].to_csv(col_metadata_file, sep=",", index=True)
+    else:
+        obs_for_deseq.to_csv(col_metadata_file, sep=",", index=True)
+    # write var to a gene_metadata file
+    pdata.var.to_csv(f"{output_dir}gene_metadata.csv", sep=",", index=True)
+    # write the counts matrix of a specified layer to file
+    if layer is None:
+        # write the X numpy matrix transposed to file
+        df = pd.DataFrame(pdata.X.T, index=pdata.var.index, columns=obs_for_deseq.index)
+    else:
+        df = pd.DataFrame(
+            pdata.layers[layer].T, index=pdata.var.index, columns=obs_for_deseq.index
+        )
+    df.to_csv(f"{output_dir}counts_matrix.csv", sep=",", index_label="")
+
+
+def plot_pseudobulk_samples(
+    pseudobulk_data,
+    groupby,
+    figsize=(10, 10),
+    save_path=None,
+):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> adata.X = abs(adata.X).astype(int)
+    >>> pseudobulk = get_pseudobulk(adata, "bulk_labels", "louvain")
+    >>> plot_pseudobulk_samples(pseudobulk, groupby=["bulk_labels", "louvain"], figsize=(10, 10))
+    """
+    fig = decoupler.plot_psbulk_samples(
+        pseudobulk_data, groupby=groupby, figsize=figsize, return_fig=True
+    )
+    if save_path:
+        fig.savefig(f"{save_path}/pseudobulk_samples.png")
+    else:
+        fig.show()
+
+
+def plot_filter_by_expr(
+    pseudobulk_data, group, min_count=None, min_total_count=None, save_path=None
+):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> adata.X = abs(adata.X).astype(int)
+    >>> pseudobulk = get_pseudobulk(adata, "bulk_labels", "louvain")
+    >>> plot_filter_by_expr(pseudobulk, group="bulk_labels", min_count=10, min_total_count=200)
+    """
+    fig = decoupler.plot_filter_by_expr(
+        pseudobulk_data,
+        group=group,
+        min_count=min_count,
+        min_total_count=min_total_count,
+        return_fig=True,
+    )
+    if save_path:
+        fig.savefig(f"{save_path}/filter_by_expr.png")
+    else:
+        fig.show()
+
+
+def filter_by_expr(pdata, min_count=None, min_total_count=None):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> adata.X = abs(adata.X).astype(int)
+    >>> pseudobulk = get_pseudobulk(adata, "bulk_labels", "louvain")
+    >>> pdata_filt = filter_by_expr(pseudobulk, min_count=10, min_total_count=200)
+    """
+    genes = decoupler.filter_by_expr(
+        pdata, min_count=min_count, min_total_count=min_total_count
+    )
+    return pdata[:, genes].copy()
+
+
+def check_fields(fields, adata, obs=True, context=None):
+    """
+    >>> import scanpy as sc
+    >>> adata = sc.datasets.pbmc68k_reduced()
+    >>> check_fields(["bulk_labels", "louvain"], adata, obs=True)
+    """
+
+    legend = ""
+    if context:
+        legend = f", passed in {context},"
+    if obs:
+        if not set(fields).issubset(set(adata.obs.columns)):
+            raise ValueError(
+                f"Some of the following fields {legend} are not present in adata.obs: {fields}. Possible fields are: {list(set(adata.obs.columns))}"
+            )
+    else:
+        if not set(fields).issubset(set(adata.var.columns)):
+            raise ValueError(
+                f"Some of the following fields {legend} are not present in adata.var: {fields}. Possible fields are: {list(set(adata.var.columns))}"
+            )
+
+
+def main(args):
+    # Load AnnData object from file
+    adata = anndata.read_h5ad(args.adata_file)
+
+    # Merge adata.obs fields specified in args.adata_obs_fields_to_merge
+    if args.adata_obs_fields_to_merge:
+        fields = args.adata_obs_fields_to_merge.split(",")
+        check_fields(fields, adata)
+        adata = merge_adata_obs_fields(fields, adata)
+
+    check_fields([args.groupby, args.sample_key], adata)
+
+    factor_fields = None
+    if args.factor_fields:
+        factor_fields = args.factor_fields.split(",")
+        check_fields(factor_fields, adata)
+
+    print(f"Using mode: {args.mode}")
+    # Perform pseudobulk analysis
+    pseudobulk_data = get_pseudobulk(
+        adata,
+        sample_col=args.sample_key,
+        groups_col=args.groupby,
+        layer=args.layer,
+        mode=args.mode,
+        use_raw=args.use_raw,
+        min_cells=args.min_cells,
+        min_counts=args.min_counts,
+    )
+
+    # Plot pseudobulk samples
+    plot_pseudobulk_samples(
+        pseudobulk_data,
+        args.groupby,
+        save_path=args.save_path,
+        figsize=args.plot_samples_figsize,
+    )
+
+    plot_filter_by_expr(
+        pseudobulk_data,
+        group=args.groupby,
+        min_count=args.min_counts,
+        min_total_count=args.min_total_counts,
+        save_path=args.save_path,
+    )
+
+    # Filter by expression if enabled
+    if args.filter_expr:
+        filtered_adata = filter_by_expr(
+            pseudobulk_data,
+            min_count=args.min_counts,
+            min_total_count=args.min_total_counts,
+        )
+
+        pseudobulk_data = filtered_adata
+
+    # Save the pseudobulk data
+    if args.anndata_output_path:
+        pseudobulk_data.write_h5ad(args.anndata_output_path, compression="gzip")
+
+    write_DESeq2_inputs(
+        pseudobulk_data, output_dir=args.deseq2_output_path, factor_fields=factor_fields
+    )
+
+
+def merge_adata_obs_fields(obs_fields_to_merge, adata):
+    """
+    Merge adata.obs fields specified in args.adata_obs_fields_to_merge
+
+    Parameters
+    ----------
+    obs_fields_to_merge : str
+        Fields in adata.obs to merge, comma separated
+    adata : anndata.AnnData
+        The AnnData object
+
+    Returns
+    -------
+    anndata.AnnData
+        The merged AnnData object
+
+    docstring tests:
+    >>> import scanpy as sc
+    >>> ad = sc.datasets.pbmc68k_reduced()
+    >>> ad = merge_adata_obs_fields(["bulk_labels","louvain"], ad)
+    >>> ad.obs.columns
+    Index(['bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score',
+           'G2M_score', 'phase', 'louvain', 'bulk_labels_louvain'],
+          dtype='object')
+    """
+    field_name = "_".join(obs_fields_to_merge)
+    for field in obs_fields_to_merge:
+        if field not in adata.obs.columns:
+            raise ValueError(f"The '{field}' column is not present in adata.obs.")
+        if field_name not in adata.obs.columns:
+            adata.obs[field_name] = adata.obs[field].astype(str)
+        else:
+            adata.obs[field_name] = (
+                adata.obs[field_name] + "_" + adata.obs[field].astype(str)
+            )
+    return adata
+
+
+if __name__ == "__main__":
+    # Create argument parser
+    parser = argparse.ArgumentParser(
+        description="Perform pseudobulk analysis on an AnnData object"
+    )
+
+    # Add arguments
+    parser.add_argument("adata_file", type=str, help="Path to the AnnData file")
+    parser.add_argument(
+        "-m",
+        "--adata_obs_fields_to_merge",
+        type=str,
+        help="Fields in adata.obs to merge, comma separated",
+    )
+    parser.add_argument(
+        "--groupby",
+        type=str,
+        required=True,
+        help="The column in adata.obs that defines the groups",
+    )
+    parser.add_argument(
+        "--sample_key",
+        required=True,
+        type=str,
+        help="The column in adata.obs that defines the samples",
+    )
+    # add argument for layer
+    parser.add_argument(
+        "--layer",
+        type=str,
+        default=None,
+        help="The name of the layer of the AnnData object to use",
+    )
+    # add argument for mode
+    parser.add_argument(
+        "--mode",
+        type=str,
+        default="sum",
+        help="The mode for Decoupler pseudobulk analysis",
+        choices=["sum", "mean", "median"],
+    )
+    # add boolean argument for use_raw
+    parser.add_argument(
+        "--use_raw",
+        action="store_true",
+        default=False,
+        help="Whether to use the raw part of the AnnData object",
+    )
+    # add argument for min_cells
+    parser.add_argument(
+        "--min_cells",
+        type=int,
+        default=10,
+        help="Minimum number of cells for pseudobulk analysis",
+    )
+    parser.add_argument(
+        "--save_path", type=str, help="Path to save the plot (optional)"
+    )
+    parser.add_argument(
+        "--min_counts",
+        type=int,
+        help="Minimum count threshold for filtering by expression",
+    )
+    parser.add_argument(
+        "--min_total_counts",
+        type=int,
+        help="Minimum total count threshold for filtering by expression",
+    )
+    parser.add_argument(
+        "--anndata_output_path",
+        type=str,
+        help="Path to save the filtered AnnData object or pseudobulk data",
+    )
+    parser.add_argument(
+        "--filter_expr", action="store_true", help="Enable filtering by expression"
+    )
+    parser.add_argument(
+        "--factor_fields",
+        type=str,
+        help="Comma separated list of fields for the factors",
+    )
+    parser.add_argument(
+        "--deseq2_output_path",
+        type=str,
+        help="Path to save the DESeq2 inputs",
+        required=True,
+    )
+    parser.add_argument(
+        "--plot_samples_figsize",
+        type=int,
+        default=[10, 10],
+        nargs=2,
+        help="Size of the samples plot as a tuple (two arguments)",
+    )
+    parser.add_argument("--plot_filtering_figsize", type=int, default=[10, 10], nargs=2)
+
+    # Parse the command line arguments
+    args = parser.parse_args()
+
+    # Call the main function
+    main(args)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/decoupler_pseudobulk.xml	Sun Sep 24 08:44:24 2023 +0000
@@ -0,0 +1,161 @@
+<tool id="decoupler_pseudobulk" name="Decoupler pseudo-bulk" version="1.4.0+galaxy0" profile="20.05">
+    <description>aggregates single cell RNA-seq data for running bulk RNA-seq methods</description>
+    <requirements>
+        <requirement type="package" version="1.4.0">decoupler</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+mkdir deseq_output_dir &&
+mkdir plots_output_dir &&
+python '$__tool_directory__/decoupler_pseudobulk.py' $input_file
+    #if $adata_obs_fields_to_merge:
+    --adata_obs_fields_to_merge $adata_obs_fields_to_merge
+    #end if
+    --groupby $groupby
+    --sample_key $sample_key
+    #if $layer:
+    --layer $layer
+    #end if
+    --mode $mode
+    #if $use_raw:
+    --use_raw
+    #end if
+    #if $min_cells:
+    --min_cells $min_cells
+    #end if
+    #if $produce_plots:
+    --save_path plots_output_dir
+    #end if
+    #if $min_counts:
+    --min_counts $min_counts
+    #end if
+    #if $min_total_counts:
+    --min_total_counts $min_total_counts
+    #end if
+    #if $produce_anndata:
+    --anndata_output_path $pbulk_anndata
+    #end if
+    #if $filter_expr:
+    --filter_expr
+    #end if
+    #if $factor_fields:
+    --factor_fields '$factor_fields'
+    #end if
+    --deseq2_output_path deseq_output_dir
+    --plot_samples_figsize $plot_samples_figsize
+    --plot_filtering_figsize $plot_filtering_figsize
+]]></command>
+    <environment_variables>
+        <environment_variable name="NUMBA_CACHE_DIR">\$_GALAXY_JOB_TMP_DIR</environment_variable>
+        <environment_variable name="MPLCONFIGDIR">\$_GALAXY_JOB_TMP_DIR</environment_variable>
+    </environment_variables>
+    <inputs>
+        <param type="data" name="input_file" format="data" label="Input AnnData file"/>
+        <param type="text" name="adata_obs_fields_to_merge" label="Obs Fields to Merge" optional="true" help="Fields in adata.obs to merge, comma separated (optional). They will be available as field1_field2_field3 in the AnnData Obs dataframe."/>
+        <param type="text" name="groupby" label="Groupby column" help="The column in adata.obs that defines the groups. Merged columns in the above field are available here."/>
+        <param type="text" name="sample_key" label="Sample Key column" help="The column in adata.obs that defines the samples. Merged columns in the above field are available here."/>
+        <param type="text" name="layer" label="Layer" optional="true" help="The name of the layer of the AnnData object to use. It needs to be present in the AnnData object."/>
+        <param type="select" name="mode" label="Decoupler pseudobulk Mode" optional="true" help="Determines how counts are aggregated across cells with the specificied groups: sum, mean or median.">
+            <option value="sum" selected="true">sum</option>
+            <option value="mean">mean</option>
+            <option value="median">median</option>
+        </param>
+        <param type="text" name="factor_fields" label="Factor Fields" optional="true" help="Fields in adata.obs to use as factors, comma separated (optional). For EdgeR make sure that the first field is the main contrast field desired and the rest of the fields are the covariates desired. Decoupler produces two fields in the intermediate AnnData (which can be added here if desired for covariates): psbulk_n_cells and psbulk_counts."/>
+        <param type="boolean" name="use_raw" label="Use Raw" optional="true"/>
+        <param type="integer" name="min_cells" label="Minimum Cells" optional="true"/>
+        <param type="boolean" name="produce_plots" label="Produce plots"/>
+        <param type="boolean" name="produce_anndata" label="Produce AnnData with Pseudo-bulk"/>
+        <param type="integer" name="min_counts" label="Minimum Counts" optional="true"/>
+        <param type="integer" name="min_total_counts" label="Minimum Total Counts" optional="true"/>
+        <param type="boolean" name="filter_expr" label="Enable Filtering by Expression"/>
+        <param type="text" name="plot_samples_figsize" label="Plot Samples Figsize" value="10 10" help="X and Y sizes in points separated by a space"/>
+        <param type="text" name="plot_filtering_figsize" label="Plot Filtering Figsize" value="10 10" help="X and Y sizes in points separated by a space"/>
+    </inputs>
+    <outputs>
+        <data name="pbulk_anndata" format="h5ad" label="${tool.name} on ${on_string}: Pseudo-bulk AnnData">
+            <filter>produce_anndata</filter>
+        </data>
+        <data name="count_matrix" format="csv" label="${tool.name} on ${on_string}: Count Matrix" from_work_dir="deseq_output_dir/counts_matrix.csv"/>
+        <data name="samples_metadata" format="csv" label="${tool.name} on ${on_string}: Samples Metadata (factors file)" from_work_dir="deseq_output_dir/col_metadata.csv"/>
+        <data name="genes_metadata" format="csv" label="${tool.name} on ${on_string}: Genes Metadata" from_work_dir="deseq_output_dir/gene_metadata.csv"/>
+        <data name="plot_output" format="png" label="${tool.name} on ${on_string}: Pseudobulk plot" from_work_dir="plots_output_dir/pseudobulk_samples.png">
+            <filter>produce_plots</filter>
+        </data>
+        <data name="filter_by_expr_plot" format="png" label="${tool.name} on ${on_string}: Filter by Expression plot" from_work_dir="plots_output_dir/filter_by_expr.png">
+            <filter>produce_plots</filter>
+        </data>
+    </outputs>
+    <tests>
+        <test expect_num_outputs="6">
+            <param name="input_file" value="mito_counted_anndata.h5ad"/>
+            <param name="adata_obs_fields_to_merge" value="batch,sex"/>
+            <param name="groupby" value="batch_sex"/>
+            <param name="sample_key" value="genotype"/>
+            <param name="factor_fields" value="genotype,batch_sex"/>
+            <param name="mode" value="sum"/>
+            <param name="min_cells" value="10"/>
+            <param name="produce_plots" value="true"/>
+            <param name="produce_anndata" value="true"/>
+            <param name="min_counts" value="10"/>
+            <param name="min_total_counts" value="1000"/>
+            <param name="filter_expr" value="true"/>
+            <param name="plot_samples_figsize" value="10 10"/>
+            <param name="plot_filtering_figsize" value="10 10"/>
+            <output name="pbulk_anndata" ftype="h5ad">
+                <assert_contents>
+                    <has_h5_keys keys="obs/psbulk_n_cells"/>
+                </assert_contents>
+            </output>
+            <output name="count_matrix" ftype="csv">
+                <assert_contents>
+                    <has_n_lines n="3620"/>
+                </assert_contents>
+            </output>
+            <output name="samples_metadata" ftype="csv">
+                <assert_contents>
+                    <has_n_lines n="8"/>
+                </assert_contents>
+            </output>
+            <output name="genes_metadata" ftype="csv">
+                <assert_contents>
+                    <has_n_lines n="3620"/>
+                </assert_contents>
+            </output>
+            <output name="plot_output" ftype="png">
+                <assert_contents>
+                    <has_size value="31853" delta="3000"/>
+                </assert_contents>
+            </output>
+            <output name="filter_by_expr_plot" ftype="png">
+                <assert_contents>
+                    <has_size value="21656" delta="2000"/>
+                </assert_contents>
+            </output>
+        </test>
+    </tests>
+    <help>
+        <![CDATA[
+        This tool performs pseudobulk analysis and filtering using Decoupler-py. Provide the input AnnData file and specify the necessary parameters.
+
+        - Input AnnData file: The input AnnData file to be processed.
+        - Obs Fields to Merge: Fields in adata.obs to merge, comma separated (optional).
+        - Groupby column: The column in adata.obs that defines the groups.
+        - Sample Key column: The column in adata.obs that defines the samples.
+        - Layer (optional): The name of the layer of the AnnData object to use.
+        - Mode: The mode for Decoupler pseudobulk analysis (sum, mean, median). Sum by default.
+        - Factor Fields (optional): Fields in adata.obs to use as factors, comma separated (optional). For EdgeR make sure that the first field is the main contrast field desired and the rest of the fields are the covariates desired.
+        - Use Raw: Whether to use the raw part of the AnnData object.
+        - Minimum Cells: Minimum number of cells for pseudobulk analysis (optional).
+        - Minimum Counts: Minimum count threshold for filtering by expression (optional).
+        - Minimum Total Counts: Minimum total count threshold for filtering by expression (optional).
+        - Enable Filtering by Expression: Check this box to enable filtering by expression.
+        - Plot Samples Figsize: Size of the samples plot as a tuple (two arguments).
+        - Plot Filtering Figsize: Size of the filtering plot as a tuple (two arguments).
+
+        The tool will output the filtered AnnData, count matrix, samples metadata, genes metadata (in DESeq2 format), and the pseudobulk plot and filter by expression plot (if enabled).
+
+        ]]>
+    </help>
+    <citations>
+        <citation type="doi">doi.org/10.1093/bioadv/vbac016</citation>
+    </citations>
+</tool>

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/get_test_data.sh	Sun Sep 24 08:44:24 2023 +0000
@@ -0,0 +1,21 @@
+#!/usr/bin/env bash
+
+BASENAME_FILE='mito_counted_anndata.h5ad'
+
+MTX_LINK='https://zenodo.org/record/7053673/files/Mito-counted_AnnData'
+
+# convenience for getting data
+function get_data {
+  local link=$1
+  local fname=$2
+
+  if [ ! -f $fname ]; then
+    echo "$fname not available locally, downloading.."
+    wget -O $fname --retry-connrefused --waitretry=1 --read-timeout=20 --timeout=15 -t 3 $link
+  fi
+}
+
+# get matrix data
+mkdir -p test-data
+pushd test-data
+get_data $MTX_LINK $BASENAME_FILE