comparison decoupler_aucell_score.py @ 6:ed2a77422e00 draft

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

5:893ff9213a34

    # Convert the list of dictionaries to a DataFrame
    return pd.DataFrame(gene_sets)


def score_genes_aucell(
    adata: anndata.AnnData, gene_list: list, score_name: str, use_raw=False
):
    """Score genes using Aucell.

    Parameters
    ----------

            "gene_id": gene_list,
            "geneset": score_name,
        }
    )
    # run decoupler's run_aucell
    dc.run_aucell(
        adata, net=geneset_df, source="geneset", target="gene_id", use_raw=use_raw
    )
    # copy .obsm['aucell_estimate'] matrix columns to adata.obs using the column names
    adata.obs[score_name] = adata.obsm["aucell_estimate"][score_name]


def run_for_genelists(
    adata, gene_lists, score_names, use_raw=False, gene_symbols_field="gene_symbols"
):
    if len(gene_lists) == len(score_names):
        for gene_list, score_names in zip(gene_lists, score_names):
            genes = gene_list.split(",")
            ens_gene_ids = adata.var[adata.var[gene_symbols_field].isin(genes)].index
            score_genes_aucell(
                adata,
                ens_gene_ids,
                f"AUCell_{score_names}",
                use_raw,
            )
    else:
        raise ValueError(
            "The number of gene lists (separated by :) and score names (separated by :) must be the same"
        )

        "--gene_symbols_field",
        type=str,
        help="Name of the gene symbols field in the AnnData object",
        required=True,
    )
    parser.add_argument("--use_raw", action="store_true", help="Use raw data")
    parser.add_argument(
        "--write_anndata", action="store_true", help="Write the modified AnnData object"
    )

    if args.gmt_file is not None:
        # Load MSigDB file in GMT format
        msigdb = read_gmt(args.gmt_file)

        gene_sets_to_score = args.gene_sets_to_score.split(",") if args.gene_sets_to_score else []
        # Score genes by their ensembl ids using the score_genes_aucell function
        for _, row in msigdb.iterrows():
            gene_set_name = row["gene_set_name"]
            if not gene_sets_to_score or gene_set_name in gene_sets_to_score:
                genes = row["genes"].split(",")
                # Convert gene symbols to ensembl ids by using the columns gene_symbols and index in adata.var specific to the gene set
                ens_gene_ids = adata.var[
                    adata.var[args.gene_symbols_field].isin(genes)
                ].index
                score_genes_aucell(
                    adata, ens_gene_ids, f"AUCell_{gene_set_name}", args.use_raw
                )
    elif args.gene_lists_to_score is not None and args.score_names is not None:
        gene_lists = args.gene_lists_to_score.split(":")
        score_names = args.score_names.split(",")
        run_for_genelists(
            adata, gene_lists, score_names, args.use_raw, args.gene_symbols_field
        )

    # Save the modified AnnData object or generate a file with cells as rows and the new score_names columns
    if args.write_anndata:
        adata.write_h5ad(args.output_file)

6:ed2a77422e00

    # Convert the list of dictionaries to a DataFrame
    return pd.DataFrame(gene_sets)


def score_genes_aucell(
    adata: anndata.AnnData, gene_list: list, score_name: str, use_raw=False, min_n_genes=5
):
    """Score genes using Aucell.

    Parameters
    ----------

            "gene_id": gene_list,
            "geneset": score_name,
        }
    )
    # run decoupler's run_aucell
    # catch the value error
    try:
        dc.run_aucell(
            adata, net=geneset_df, source="geneset", target="gene_id", use_raw=use_raw
        )
        # copy .obsm['aucell_estimate'] matrix columns to adata.obs using the column names
        adata.obs[score_name] = adata.obsm["aucell_estimate"][score_name]
    except ValueError as ve:
        print(f"Gene list {score_name} failed, skipping: {str(ve)}")


def run_for_genelists(
    adata, gene_lists, score_names, use_raw=False, gene_symbols_field="gene_symbols", min_n_genes=5
):
    if len(gene_lists) == len(score_names):
        for gene_list, score_names in zip(gene_lists, score_names):
            genes = gene_list.split(",")
            ens_gene_ids = adata.var[adata.var[gene_symbols_field].isin(genes)].index
            score_genes_aucell(
                adata,
                ens_gene_ids,
                f"AUCell_{score_names}",
                use_raw,
                min_n_genes
            )
    else:
        raise ValueError(
            "The number of gene lists (separated by :) and score names (separated by :) must be the same"
        )

        "--gene_symbols_field",
        type=str,
        help="Name of the gene symbols field in the AnnData object",
        required=True,
    )
    # argument for min_n Minimum of targets per source. If less, sources are removed.
    parser.add_argument(
        "--min_n",
        type=int,
        required=False,
        default=5,
        help="Minimum of targets per source. If less, sources are removed.",
    )
    parser.add_argument("--use_raw", action="store_true", help="Use raw data")
    parser.add_argument(
        "--write_anndata", action="store_true", help="Write the modified AnnData object"
    )

    if args.gmt_file is not None:
        # Load MSigDB file in GMT format
        msigdb = read_gmt(args.gmt_file)

        gene_sets_to_score = (
            args.gene_sets_to_score.split(",") if args.gene_sets_to_score else []
        )
        # Score genes by their ensembl ids using the score_genes_aucell function
        for _, row in msigdb.iterrows():
            gene_set_name = row["gene_set_name"]
            if not gene_sets_to_score or gene_set_name in gene_sets_to_score:
                genes = row["genes"].split(",")
                # Convert gene symbols to ensembl ids by using the columns gene_symbols and index in adata.var specific to the gene set
                ens_gene_ids = adata.var[
                    adata.var[args.gene_symbols_field].isin(genes)
                ].index
                score_genes_aucell(
                    adata, ens_gene_ids, f"AUCell_{gene_set_name}", args.use_raw, args.min_n
                )
    elif args.gene_lists_to_score is not None and args.score_names is not None:
        gene_lists = args.gene_lists_to_score.split(":")
        score_names = args.score_names.split(",")
        run_for_genelists(
            adata, gene_lists, score_names, args.use_raw, args.gene_symbols_field, args.min_n
        )

    # Save the modified AnnData object or generate a file with cells as rows and the new score_names columns
    if args.write_anndata:
        adata.write_h5ad(args.output_file)