Documentation for Preprocessor

scdataloader.preprocess.Preprocessor

Prepare data into training, valid and test split. Normalize raw expression values, binning or using other transform into the preset model input format.

Initializes the preprocessor and configures the workflow steps.

Your dataset should contain at least the following obs: - organism_ontology_term_id with the ontology id of the organism of your anndata - gene names in the var.index field of your anndata that map to the ensembl_gene nomenclature or the hugo gene symbols nomenclature (if the later, set is_symbol to True)

Parameters:
  • filter_gene_by_counts (int or bool, default: False ) –

    Determines whether to filter genes by counts. If int, filters genes with counts. Defaults to False.

  • filter_cell_by_counts (int or bool, default: False ) –

    Determines whether to filter cells by counts. If int, filters cells with counts. Defaults to False.

  • normalize_sum (float or bool, default: 10000.0 ) –

    Determines whether to normalize the total counts of each cell to a specific value. Defaults to 1e4.

  • log1p (bool) –

    Determines whether to apply log1p transform to the normalized data. Defaults to True.

  • n_hvg_for_postp (int or bool, default: 0 ) –

    Determines whether to subset to highly variable genes for the PCA. Defaults to False.

  • hvg_flavor (str, default: 'seurat_v3' ) –

    Specifies the flavor of highly variable genes selection. See :func:scanpy.pp.highly_variable_genes for more details. Defaults to "seurat_v3".

  • binning (int, default: None ) –

    Determines whether to bin the data into discrete values of number of bins provided.

  • result_binned_key (str, default: 'X_binned' ) –

    Specifies the key of :class:~anndata.AnnData to store the binned data. Defaults to "X_binned".

  • length_normalize (bool, default: False ) –

    Determines whether to length normalize the data. Defaults to False.

  • force_preprocess (bool, default: False ) –

    Determines whether to bypass the check of raw counts. Defaults to False.

  • min_dataset_size (int, default: 100 ) –

    The minimum size required for a dataset to be kept. Defaults to 100.

  • min_valid_genes_id (int, default: 10000 ) –

    The minimum number of valid genes to keep a dataset. Defaults to 10_000.

  • min_nnz_genes (int, default: 200 ) –

    The minimum number of non-zero genes to keep a cell. Defaults to 200.

  • maxdropamount (int, default: 50 ) –

    The maximum amount of dropped cells per dataset. (2 for 50% drop, 3 for 33% drop, etc.) Defaults to 2.

  • madoutlier (int, default: 5 ) –

    The maximum absolute deviation of the outlier samples. Defaults to 5.

  • pct_mt_outlier (int, default: 8 ) –

    The maximum percentage of mitochondrial genes outlier. Defaults to 8.

  • batch_key (str) –

    The key of :class:~anndata.AnnData.obs to use for batch information. This arg is used in the highly variable gene selection step.

  • skip_validate (bool, default: False ) –

    Determines whether to skip the validation step. Defaults to False.

  • keepdata (bool, default: False ) –

    Determines whether to keep the data in the AnnData object. Defaults to False.

Source code in scdataloader/preprocess.py
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def __init__(
    self,
    filter_gene_by_counts: Union[int, bool] = False,
    filter_cell_by_counts: Union[int, bool] = False,
    normalize_sum: float = 1e4,
    n_hvg_for_postp: int = 0,
    use_layer: Optional[str] = None,
    is_symbol: bool = False,
    hvg_flavor: str = "seurat_v3",
    binning: Optional[int] = None,
    result_binned_key: str = "X_binned",
    length_normalize: bool = False,
    force_preprocess: bool = False,
    min_dataset_size: int = 100,
    min_valid_genes_id: int = 10_000,
    min_nnz_genes: int = 200,
    maxdropamount: int = 50,
    madoutlier: int = 5,
    pct_mt_outlier: int = 8,
    batch_keys: list[str] = [
        "assay_ontology_term_id",
        "self_reported_ethnicity_ontology_term_id",
        "sex_ontology_term_id",
        "donor_id",
        "suspension_type",
    ],
    skip_validate: bool = False,
    additional_preprocess: Optional[Callable[[AnnData], AnnData]] = None,
    additional_postprocess: Optional[Callable[[AnnData], AnnData]] = None,
    do_postp: bool = True,
    organisms: list[str] = ["NCBITaxon:9606", "NCBITaxon:10090"],
    use_raw: bool = True,
    keepdata: bool = False,
) -> None:
    """
    Initializes the preprocessor and configures the workflow steps.

    Your dataset should contain at least the following obs:
    - `organism_ontology_term_id` with the ontology id of the organism of your anndata
    - gene names in the `var.index` field of your anndata that map to the ensembl_gene nomenclature
    or the hugo gene symbols nomenclature (if the later, set `is_symbol` to True)

    Args:
        filter_gene_by_counts (int or bool, optional): Determines whether to filter genes by counts.
            If int, filters genes with counts. Defaults to False.
        filter_cell_by_counts (int or bool, optional): Determines whether to filter cells by counts.
            If int, filters cells with counts. Defaults to False.
        normalize_sum (float or bool, optional): Determines whether to normalize the total counts of each cell to a specific value.
            Defaults to 1e4.
        log1p (bool, optional): Determines whether to apply log1p transform to the normalized data.
            Defaults to True.
        n_hvg_for_postp (int or bool, optional): Determines whether to subset to highly variable genes for the PCA.
            Defaults to False.
        hvg_flavor (str, optional): Specifies the flavor of highly variable genes selection.
            See :func:`scanpy.pp.highly_variable_genes` for more details. Defaults to "seurat_v3".
        binning (int, optional): Determines whether to bin the data into discrete values of number of bins provided.
        result_binned_key (str, optional): Specifies the key of :class:`~anndata.AnnData` to store the binned data.
            Defaults to "X_binned".
        length_normalize (bool, optional): Determines whether to length normalize the data.
            Defaults to False.
        force_preprocess (bool, optional): Determines whether to bypass the check of raw counts.
            Defaults to False.
        min_dataset_size (int, optional): The minimum size required for a dataset to be kept.
            Defaults to 100.
        min_valid_genes_id (int, optional): The minimum number of valid genes to keep a dataset.
            Defaults to 10_000.
        min_nnz_genes (int, optional): The minimum number of non-zero genes to keep a cell.
            Defaults to 200.
        maxdropamount (int, optional): The maximum amount of dropped cells per dataset. (2 for 50% drop, 3 for 33% drop, etc.)
            Defaults to 2.
        madoutlier (int, optional): The maximum absolute deviation of the outlier samples.
            Defaults to 5.
        pct_mt_outlier (int, optional): The maximum percentage of mitochondrial genes outlier.
            Defaults to 8.
        batch_key (str, optional): The key of :class:`~anndata.AnnData.obs` to use for batch information.
            This arg is used in the highly variable gene selection step.
        skip_validate (bool, optional): Determines whether to skip the validation step.
            Defaults to False.
        keepdata (bool, optional): Determines whether to keep the data in the AnnData object.
            Defaults to False.
    """
    self.filter_gene_by_counts = filter_gene_by_counts
    self.filter_cell_by_counts = filter_cell_by_counts
    self.normalize_sum = normalize_sum
    self.hvg_flavor = hvg_flavor
    self.binning = binning
    self.organisms = organisms
    self.result_binned_key = result_binned_key
    self.additional_preprocess = additional_preprocess
    self.additional_postprocess = additional_postprocess
    self.force_preprocess = force_preprocess
    self.min_dataset_size = min_dataset_size
    self.min_valid_genes_id = min_valid_genes_id
    self.min_nnz_genes = min_nnz_genes
    self.maxdropamount = maxdropamount
    self.madoutlier = madoutlier
    self.n_hvg_for_postp = n_hvg_for_postp
    self.pct_mt_outlier = pct_mt_outlier
    self.batch_keys = batch_keys
    self.length_normalize = length_normalize
    self.skip_validate = skip_validate
    self.use_layer = use_layer
    self.is_symbol = is_symbol
    self.do_postp = do_postp
    self.use_raw = use_raw
    self.keepdata = keepdata

scdataloader.preprocess.LaminPreprocessor

Bases: Preprocessor

Methods:

Name Description
__call__

format controls the different input value wrapping, including categorical

Source code in scdataloader/preprocess.py
430
431
432
433
434
435
436
437
438
439
440
441
def __init__(
    self,
    *args,
    cache: bool = True,
    keep_files: bool = True,
    force_preloaded: bool = False,
    **kwargs,
):
    super().__init__(*args, **kwargs)
    self.cache = cache
    self.keep_files = keep_files
    self.force_preloaded = force_preloaded

__call__

format controls the different input value wrapping, including categorical binned style, fixed-sum normalized counts, log1p fixed-sum normalized counts, etc.

Parameters:
  • adata (AnnData) –

    The AnnData object to preprocess.

  • batch_key (str) –

    The key of AnnData.obs to use for batch information. This arg is used in the highly variable gene selection step.

Source code in scdataloader/preprocess.py
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
def __call__(
    self,
    data: Union[ln.Collection, AnnData] = None,
    name: str = "preprocessed dataset",
    description: str = "preprocessed dataset using scprint",
    start_at: int = 0,
    version: str = "2",
):
    """
    format controls the different input value wrapping, including categorical
    binned style, fixed-sum normalized counts, log1p fixed-sum normalized counts, etc.

    Args:
        adata (AnnData): The AnnData object to preprocess.
        batch_key (str, optional): The key of AnnData.obs to use for batch information. This arg
            is used in the highly variable gene selection step.
    """
    files = []
    all_ready_processed_keys = set()
    if self.cache:
        for i in ln.Artifact.filter(description=description):
            all_ready_processed_keys.add(i.stem_uid)
    if isinstance(data, AnnData):
        return super().__call__(data)
    elif isinstance(data, ln.Collection):
        for i, file in enumerate(data.artifacts.all()[start_at:]):
            # use the counts matrix
            i = i + start_at
            print(i)
            if file.stem_uid in all_ready_processed_keys:
                print(f"{file.stem_uid} is already processed... not preprocessing")
                continue
            print(file)

            path = cache_path(file) if self.force_preloaded else file.cache()
            backed = read_h5ad(path, backed="r")
            if backed.obs.is_primary_data.sum() == 0:
                print(f"{file.key} only contains non primary cells.. dropping")
                # Save the stem_uid to a file to avoid loading it again
                with open("nonprimary.txt", "a") as f:
                    f.write(f"{file.stem_uid}\n")
                continue
            if backed.shape[1] < 1000:
                print(
                    f"{file.key} only contains less than 1000 genes and is likely not scRNAseq... dropping"
                )
                continue
            if file.size <= MAXFILESIZE:
                adata = backed.to_memory()
                print(adata)
            else:
                badata = backed
                print(badata)
            try:
                if file.size > MAXFILESIZE:
                    print(
                        f"dividing the dataset as it is too large: {file.size // 1_000_000_000}Gb"
                    )
                    num_blocks = int(np.ceil(file.size / (MAXFILESIZE / 2)))
                    block_size = int(
                        (np.ceil(badata.shape[0] / 30_000) * 30_000) // num_blocks
                    )
                    print("num blocks ", num_blocks)
                    for j in range(num_blocks):
                        if j == 0 and i == 390:
                            continue
                        start_index = j * block_size
                        end_index = min((j + 1) * block_size, badata.shape[0])
                        block = badata[start_index:end_index].to_memory()
                        print(block)
                        block = super().__call__(
                            block, dataset_id=file.stem_uid + "_p" + str(j)
                        )
                        myfile = ln.Artifact.from_anndata(
                            block,
                            description=description
                            + " n"
                            + str(i)
                            + " p"
                            + str(j)
                            + " ( revises file "
                            + str(file.key)
                            + " )",
                            version=version,
                        )
                        myfile.save()
                        if self.keep_files:
                            files.append(myfile)
                        else:
                            del myfile
                            del block

                else:
                    adata = super().__call__(adata, dataset_id=file.stem_uid)
                    myfile = ln.Artifact.from_anndata(
                        adata,
                        revises=file,
                        description=description + " p" + str(i),
                        version=version,
                    )
                    myfile.save()
                    if self.keep_files:
                        files.append(myfile)
                    else:
                        del myfile
                        del adata

            except ValueError as v:
                if v.args[0].startswith("we cannot work with this organism"):
                    print(v)
                    continue
                else:
                    raise v
            except Exception as e:
                if e.args[0].startswith("Dataset dropped due to"):
                    print(e)
                    continue
                else:
                    raise e

            # issues with KLlggfw6I6lvmbqiZm46
        if self.keep_files:
            dataset = ln.Collection(files, name=name, description=description)
            dataset.save()
            return dataset
        else:
            return
    else:
        raise ValueError("Please provide either anndata or ln.Collection")

scdataloader.preprocess.additional_preprocess

Source code in scdataloader/preprocess.py
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
def additional_preprocess(adata):
    adata.obs = adata.obs.replace(
        {
            "self_reported_ethnicity_ontology_term_id": {
                "multiethnic": "unknown",
                "American": "unknown",
                "Jewish Israeli": "unknown",
                "na": "unknown",
            }
        }
    )  # multi ethnic will have to get renamed
    adata.obs["cell_culture"] = False
    # if cell_type contains the word "(cell culture)" then it is a cell culture and we mark it as so and remove this from the cell type
    loc = adata.obs["cell_type_ontology_term_id"].str.contains(
        "(cell culture)", regex=False
    )
    if loc.sum() > 0:
        adata.obs["cell_type_ontology_term_id"] = adata.obs[
            "cell_type_ontology_term_id"
        ].astype(str)
        adata.obs.loc[loc, "cell_culture"] = True
        adata.obs.loc[loc, "cell_type_ontology_term_id"] = adata.obs.loc[
            loc, "cell_type_ontology_term_id"
        ].str.replace(" (cell culture)", "")

    loc = adata.obs["tissue_ontology_term_id"].str.contains(
        "(cell culture)", regex=False
    )
    if loc.sum() > 0:
        adata.obs.loc[loc, "cell_culture"] = True
        adata.obs["tissue_ontology_term_id"] = adata.obs[
            "tissue_ontology_term_id"
        ].astype(str)
        adata.obs.loc[loc, "tissue_ontology_term_id"] = adata.obs.loc[
            loc, "tissue_ontology_term_id"
        ].str.replace(" (cell culture)", "")

    loc = adata.obs["tissue_ontology_term_id"].str.contains("(organoid)", regex=False)
    if loc.sum() > 0:
        adata.obs.loc[loc, "cell_culture"] = True
        adata.obs["tissue_ontology_term_id"] = adata.obs[
            "tissue_ontology_term_id"
        ].astype(str)
        adata.obs.loc[loc, "tissue_ontology_term_id"] = adata.obs.loc[
            loc, "tissue_ontology_term_id"
        ].str.replace(" (organoid)", "")

    loc = adata.obs["tissue_ontology_term_id"].str.contains("CL:", regex=False)
    if loc.sum() > 0:
        adata.obs["tissue_ontology_term_id"] = adata.obs[
            "tissue_ontology_term_id"
        ].astype(str)
        adata.obs.loc[loc, "tissue_ontology_term_id"] = "unknown"
    return adata

scdataloader.preprocess.additional_postprocess

Source code in scdataloader/preprocess.py
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
def additional_postprocess(adata):
    # import palantir

    # define the "up to" 10 neighbors for each cells and add to obs
    # compute neighbors
    # need to be connectivities and same labels [cell type, assay, dataset, disease]
    # define the "neighbor" up to 10(N) cells and add to obs
    # define the "next time point" up to 5(M) cells and add to obs  # step 1: filter genes
    # if len(adata.obs["batches"].unique()) > 1:
    #    sc.external.pp.harmony_integrate(adata, key="batches")
    #    sc.pp.neighbors(adata, use_rep="X_pca_harmony")
    # else:
    sc.pp.neighbors(adata, use_rep="X_pca")
    sc.tl.leiden(adata, key_added="leiden_2", resolution=2.0)
    sc.tl.leiden(adata, key_added="leiden_1", resolution=1.0)
    sc.tl.leiden(adata, key_added="leiden_0.5", resolution=0.5)
    sc.tl.umap(adata)
    mid = adata.uns["dataset_id"] if "dataset_id" in adata.uns else "unknown_id"
    sc.pl.umap(
        adata,
        ncols=1,
        color=["cell_type", "batches"],
        save="_" + mid + ".png",
    )
    COL = "cell_type_ontology_term_id"
    NEWOBS = "clust_cell_type"
    MINCELLS = 10
    MAXSIM = 0.94
    from collections import Counter

    from .config import MAIN_HUMAN_MOUSE_DEV_STAGE_MAP

    adata.obs[NEWOBS] = (
        adata.obs[COL].astype(str) + "_" + adata.obs["leiden_1"].astype(str)
    )
    coun = Counter(adata.obs[NEWOBS])
    relab = {}
    for i in adata.obs[COL].unique():
        num = 0
        for n, c in sorted(coun.items(), key=lambda x: x[1], reverse=True):
            if i in n:
                if c < MINCELLS or num == 0:
                    relab[n] = i
                else:
                    relab[n] = i + "_" + str(num)
                num += 1

    adata.obs[NEWOBS] = adata.obs[NEWOBS].map(relab)

    cluster_means = pd.DataFrame(
        np.array(
            [
                adata.X[adata.obs[NEWOBS] == i].mean(axis=0)
                for i in adata.obs[NEWOBS].unique()
            ]
        )[:, 0, :],
        index=adata.obs[NEWOBS].unique(),
    )

    # Calculate correlation matrix between clusters
    cluster_similarity = cluster_means.T.corr()
    cluster_similarity.values[np.tril_indices(len(cluster_similarity), -1)] = 0

    # Get pairs with similarity > 0.95
    high_sim_pairs = []
    for i in range(len(cluster_similarity)):
        for j in range(i + 1, len(cluster_similarity)):
            if (
                cluster_similarity.iloc[i, j] > MAXSIM
                and cluster_similarity.columns[i].split("_")[0]
                == cluster_similarity.columns[j].split("_")[0]
            ):
                high_sim_pairs.append(
                    (
                        cluster_similarity.index[i],
                        cluster_similarity.columns[j],
                    )
                )
    # Create mapping for merging similar clusters
    merge_mapping = {}
    for pair in high_sim_pairs:
        if pair[0] not in merge_mapping:
            merge_mapping[pair[1]] = pair[0]
        else:
            merge_mapping[pair[1]] = merge_mapping[pair[0]]

    # Apply merging
    adata.obs[NEWOBS] = adata.obs[NEWOBS].map(merge_mapping).fillna(adata.obs[NEWOBS])
    adata.obs[NEWOBS] = adata.obs[NEWOBS].astype(str)
    coun = Counter(adata.obs[NEWOBS]).most_common()
    merge_mapping = {}
    for i in adata.obs[COL].unique():
        num = 0
        for j, c in coun:
            if i in j:
                merge_mapping[j] = i + "_" + str(num) if num > 0 else i
                num += 1
    adata.obs[NEWOBS] = adata.obs[NEWOBS].map(merge_mapping).fillna(adata.obs[NEWOBS])

    import bionty as bt

    stages = adata.obs["development_stage_ontology_term_id"].unique()
    if adata.obs.organism_ontology_term_id.unique() == ["NCBITaxon:9606"]:
        relabel = {i: i for i in stages}
        for stage in stages:
            stage_obj = bt.DevelopmentalStage.filter(ontology_id=stage).first()
            parents = set([i.ontology_id for i in stage_obj.parents.filter()])
            parents = parents - set(
                [
                    "HsapDv:0010000",
                    "HsapDv:0000204",
                    "HsapDv:0000227",
                ]
            )
            if len(parents) > 0:
                for p in parents:
                    if p in MAIN_HUMAN_MOUSE_DEV_STAGE_MAP:
                        relabel[stage] = p
        adata.obs["simplified_dev_stage"] = adata.obs[
            "development_stage_ontology_term_id"
        ].map(relabel)
    elif adata.obs.organism_ontology_term_id.unique() == ["NCBITaxon:10090"]:
        rename_mapping = {
            k: v for v, j in MAIN_HUMAN_MOUSE_DEV_STAGE_MAP.items() for k in j
        }
        relabel = {i: "unknown" for i in stages}
        for stage in stages:
            if stage in rename_mapping:
                relabel[stage] = rename_mapping[stage]
        adata.obs["simplified_dev_stage"] = adata.obs[
            "development_stage_ontology_term_id"
        ].map(relabel)
    else:
        raise ValueError("organism not supported")
    # palantir.utils.run_diffusion_maps(adata, n_components=20)
    # palantir.utils.determine_multiscale_space(adata)
    # terminal_states = palantir.utils.find_terminal_states(
    #    adata,
    #    celltypes=adata.obs.cell_type_ontology_term_id.unique(),
    #    celltype_column="cell_type_ontology_term_id",
    # )
    # sc.tl.diffmap(adata)
    # adata.obs["heat_diff"] = 1
    # for terminal_state in terminal_states.index.tolist():
    #    adata.uns["iroot"] = np.where(adata.obs.index == terminal_state)[0][0]
    #    sc.tl.dpt(adata)
    #    adata.obs["heat_diff"] = np.minimum(
    #        adata.obs["heat_diff"], adata.obs["dpt_pseudotime"]
    #    )
    return adata