snputils

  21class SNPObject:
  22    """
  23    A class for Single Nucleotide Polymorphism (SNP) data, with optional support for 
  24    SNP-level Local Ancestry Information (LAI).
  25    """
  26    def __init__(
  27        self,
  28        calldata_gt: Optional[np.ndarray] = None,
  29        samples: Optional[np.ndarray] = None,
  30        variants_ref: Optional[np.ndarray] = None,
  31        variants_alt: Optional[np.ndarray] = None,
  32        variants_chrom: Optional[np.ndarray] = None,
  33        variants_filter_pass: Optional[np.ndarray] = None,
  34        variants_id: Optional[np.ndarray] = None,
  35        variants_pos: Optional[np.ndarray] = None,
  36        variants_qual: Optional[np.ndarray] = None,
  37        calldata_lai: Optional[np.ndarray] = None,
  38        ancestry_map: Optional[Dict[str, str]] = None
  39    ) -> None:
  40        """
  41        Args:
  42            calldata_gt (array, optional): 
  43                An array containing genotype data for each sample. This array can be either 2D with shape 
  44                `(n_snps, n_samples)` if the paternal and maternal strands are summed, or 3D with shape 
  45                `(n_snps, n_samples, 2)` if the strands are kept separate.
  46            samples (array of shape (n_samples,), optional): 
  47                An array containing unique sample identifiers.
  48            variants_ref (array of shape (n_snps,), optional): 
  49                An array containing the reference allele for each SNP.
  50            variants_alt (array of shape (n_snps,), optional): 
  51                An array containing the alternate allele for each SNP.
  52            variants_chrom (array of shape (n_snps,), optional): 
  53                An array containing the chromosome for each SNP.
  54            variants_filter_pass (array of shape (n_snps,), optional): 
  55                An array indicating whether each SNP passed control checks.
  56            variants_id (array of shape (n_snps,), optional): 
  57                An array containing unique identifiers (IDs) for each SNP.
  58            variants_pos (array of shape (n_snps,), optional): 
  59                An array containing the chromosomal positions for each SNP.
  60            variants_qual (array of shape (n_snps,), optional): 
  61                An array containing the Phred-scaled quality score for each SNP.
  62            calldata_lai (array, optional): 
  63                An array containing the ancestry for each SNP. This array can be either 2D with shape
  64                `(n_snps, n_samples*2)`, or 3D with shape (n_snps, n_samples, 2).
  65            ancestry_map (dict of str to str, optional): 
  66                A dictionary mapping ancestry codes to region names.
  67        """
  68        self.__calldata_gt = calldata_gt
  69        self.__samples = samples
  70        self.__variants_ref = variants_ref
  71        self.__variants_alt = variants_alt
  72        self.__variants_chrom = variants_chrom
  73        self.__variants_filter_pass = variants_filter_pass
  74        self.__variants_id = variants_id
  75        self.__variants_pos = variants_pos
  76        self.__variants_qual = variants_qual
  77        self.__calldata_lai = calldata_lai
  78        self.__ancestry_map = ancestry_map
  79
  80        self._sanity_check()
  81
  82    def __getitem__(self, key: str) -> Any:
  83        """
  84        To access an attribute of the class using the square bracket notation,
  85        similar to a dictionary.
  86        """
  87        try:
  88            return getattr(self, key)
  89        except:
  90            raise KeyError(f'Invalid key: {key}.')
  91
  92    def __setitem__(self, key: str, value: Any):
  93        """
  94        To set an attribute of the class using the square bracket notation,
  95        similar to a dictionary.
  96        """
  97        try:
  98            setattr(self, key, value)
  99        except:
 100            raise KeyError(f'Invalid key: {key}.')
 101
 102    @property
 103    def calldata_gt(self) -> np.ndarray:
 104        """
 105        Retrieve `calldata_gt`.
 106
 107        Returns:
 108            **array:** 
 109                An array containing genotype data for each sample. This array can be either 2D with shape 
 110                `(n_snps, n_samples)` if the paternal and maternal strands are summed, or 3D with shape 
 111                `(n_snps, n_samples, 2)` if the strands are kept separate.
 112        """
 113        return self.__calldata_gt
 114
 115    @calldata_gt.setter
 116    def calldata_gt(self, x: np.ndarray):
 117        """
 118        Update `calldata_gt`.
 119        """
 120        self.__calldata_gt = x
 121
 122    @property
 123    def samples(self) -> Optional[np.ndarray]:
 124        """
 125        Retrieve `samples`.
 126
 127        Returns:
 128            **array of shape (n_samples,):** 
 129                An array containing unique sample identifiers.
 130        """
 131        return self.__samples
 132
 133    @samples.setter
 134    def samples(self, x: Union[List, np.ndarray]):
 135        """
 136        Update `samples`.
 137        """
 138        self.__samples = np.asarray(x)
 139
 140    @property
 141    def variants_ref(self) -> Optional[np.ndarray]:
 142        """
 143        Retrieve `variants_ref`.
 144
 145        Returns:
 146            **array of shape (n_snps,):** An array containing the reference allele for each SNP.
 147        """
 148        return self.__variants_ref
 149
 150    @variants_ref.setter
 151    def variants_ref(self, x: np.ndarray):
 152        """
 153        Update `variants_ref`.
 154        """
 155        self.__variants_ref = x
 156
 157    @property
 158    def variants_alt(self) -> Optional[np.ndarray]:
 159        """
 160        Retrieve `variants_alt`.
 161
 162        Returns:
 163            **array of shape (n_snps,):** An array containing the alternate allele for each SNP.
 164        """
 165        return self.__variants_alt
 166
 167    @variants_alt.setter
 168    def variants_alt(self, x: np.ndarray):
 169        """
 170        Update `variants_alt`.
 171        """
 172        self.__variants_alt = x
 173
 174    @property
 175    def variants_chrom(self) -> Optional[np.ndarray]:
 176        """
 177        Retrieve `variants_chrom`.
 178
 179        Returns:
 180            **array of shape (n_snps,):** An array containing the chromosome for each SNP.
 181        """
 182        return self.__variants_chrom
 183
 184    @variants_chrom.setter
 185    def variants_chrom(self, x: np.ndarray):
 186        """
 187        Update `variants_chrom`.
 188        """
 189        self.__variants_chrom = x
 190
 191    @property
 192    def variants_filter_pass(self) -> Optional[np.ndarray]:
 193        """
 194        Retrieve `variants_filter_pass`.
 195
 196        Returns:
 197            **array of shape (n_snps,):** An array indicating whether each SNP passed control checks.
 198        """
 199        return self.__variants_filter_pass
 200
 201    @variants_filter_pass.setter
 202    def variants_filter_pass(self, x: np.ndarray):
 203        """
 204        Update `variants_filter_pass`.
 205        """
 206        self.__variants_filter_pass = x
 207
 208    @property
 209    def variants_id(self) -> Optional[np.ndarray]:
 210        """
 211        Retrieve `variants_id`.
 212
 213        Returns:
 214            **array of shape (n_snps,):** An array containing unique identifiers (IDs) for each SNP.
 215        """
 216        return self.__variants_id
 217
 218    @variants_id.setter
 219    def variants_id(self, x: np.ndarray):
 220        """
 221        Update `variants_id`.
 222        """
 223        self.__variants_id = x
 224
 225    @property
 226    def variants_pos(self) -> Optional[np.ndarray]:
 227        """
 228        Retrieve `variants_pos`.
 229
 230        Returns:
 231            **array of shape (n_snps,):** An array containing the chromosomal positions for each SNP.
 232        """
 233        return self.__variants_pos
 234
 235    @variants_pos.setter
 236    def variants_pos(self, x: np.ndarray):
 237        """
 238        Update `variants_pos`.
 239        """
 240        self.__variants_pos = x
 241
 242    @property
 243    def variants_qual(self) -> Optional[np.ndarray]:
 244        """
 245        Retrieve `variants_qual`.
 246
 247        Returns:
 248            **array of shape (n_snps,):** An array containing the Phred-scaled quality score for each SNP.
 249        """
 250        return self.__variants_qual
 251
 252    @variants_qual.setter
 253    def variants_qual(self, x: np.ndarray):
 254        """
 255        Update `variants_qual`.
 256        """
 257        self.__variants_qual = x
 258
 259    @property
 260    def calldata_lai(self) -> Optional[np.ndarray]:
 261        """
 262        Retrieve `calldata_lai`.
 263
 264        Returns:
 265            **array:** 
 266                An array containing the ancestry for each SNP. This array can be either 2D with shape
 267                `(n_snps, n_samples*2)`, or 3D with shape (n_snps, n_samples, 2).
 268        """
 269        return self.__calldata_lai
 270
 271    @calldata_lai.setter
 272    def calldata_lai(self, x: np.ndarray):
 273        """
 274        Update `calldata_lai`.
 275        """
 276        self.__calldata_lai = x
 277
 278    @property
 279    def ancestry_map(self) -> Optional[Dict[str, str]]:
 280        """
 281        Retrieve `ancestry_map`.
 282
 283        Returns:
 284            **dict of str to str:** A dictionary mapping ancestry codes to region names.
 285        """
 286        return self.__ancestry_map
 287
 288    @ancestry_map.setter
 289    def ancestry_map(self, x):
 290        """
 291        Update `ancestry_map`.
 292        """
 293        self.__ancestry_map = x
 294
 295    @property
 296    def n_samples(self) -> int:
 297        """
 298        Retrieve `n_samples`.
 299
 300        Returns:
 301            **int:** The total number of samples.
 302        """
 303        if self.__samples is not None:
 304            return len(self.__samples)
 305        elif self.__calldata_gt is not None:
 306            return self.__calldata_gt.shape[1]
 307        elif self.__calldata_lai is not None:
 308            if self.__calldata_lai.ndim == 2:
 309                return self.__calldata_lai.shape[1] // 2
 310            elif self.__calldata_lai.ndim == 3:
 311                return self.__calldata_lai.shape[1]
 312        else:
 313            raise ValueError("Unable to determine the total number of samples: no relevant data is available.")
 314
 315    @property
 316    def n_snps(self) -> int:
 317        """
 318        Retrieve `n_snps`.
 319
 320        Returns:
 321            **int:** The total number of SNPs.
 322        """
 323        # List of attributes that can indicate the number of SNPs
 324        potential_attributes = [
 325            self.__calldata_gt,
 326            self.__variants_ref,
 327            self.__variants_alt,
 328            self.__variants_chrom,
 329            self.__variants_filter_pass,
 330            self.__variants_id,
 331            self.__variants_pos,
 332            self.__variants_qual,
 333            self.__calldata_lai
 334        ]
 335
 336        # Check each attribute for its first dimension, which corresponds to `n_snps`
 337        for attr in potential_attributes:
 338            if attr is not None:
 339                return attr.shape[0]
 340
 341        raise ValueError("Unable to determine the total number of SNPs: no relevant data is available.")
 342
 343    @property
 344    def n_chrom(self) -> Optional[int]:
 345        """
 346        Retrieve `n_chrom`.
 347
 348        Returns:
 349            **int:** The total number of unique chromosomes in `variants_chrom`.
 350        """
 351        if self.variants_chrom is None:
 352            warnings.warn("Chromosome data `variants_chrom` is None.")
 353            return None
 354
 355        return len(self.unique_chrom)
 356
 357    @property
 358    def n_ancestries(self) -> int:
 359        """
 360        Retrieve `n_ancestries`.
 361
 362        Returns:
 363            **int:** The total number of unique ancestries.
 364        """
 365        if self.__calldata_lai is not None:
 366            return len(np.unique(self.__calldata_lai))
 367        else:
 368            raise ValueError("Unable to determine the total number of ancestries: no relevant data is available.")
 369
 370    @property
 371    def unique_chrom(self) -> Optional[np.ndarray]:
 372        """
 373        Retrieve `unique_chrom`.
 374
 375        Returns:
 376            **array:** The unique chromosome names in `variants_chrom`, preserving their order of appearance.
 377        """
 378        if self.variants_chrom is None:
 379            warnings.warn("Chromosome data `variants_chrom` is None.")
 380            return None
 381
 382        # Identify unique chromosome names and their first indexes of occurrence
 383        _, idx = np.unique(self.variants_chrom, return_index=True)
 384        # Return chromosome names sorted by their first occurrence to maintain original order
 385        return self.variants_chrom[np.sort(idx)]
 386
 387    @property
 388    def are_strands_summed(self) -> bool:
 389        """
 390        Retrieve `are_strands_summed`.
 391        
 392        Returns:
 393            **bool:** 
 394                True if the maternal and paternal strands have been summed together, which is indicated by 
 395                `calldata_gt` having shape `(n_samples, n_snps)`. False if the strands are stored separately, 
 396                indicated by `calldata_gt` having shape `(n_samples, n_snps, 2)`.
 397        """
 398        if self.calldata_gt is None:
 399            warnings.warn("Genotype data `calldata_gt` is None.")
 400            return None
 401        
 402        return self.calldata_gt.ndim == 2
 403
 404    def copy(self) -> SNPObject:
 405        """
 406        Create and return a copy of `self`.
 407
 408        Returns:
 409            **SNPObject:** 
 410                A new instance of the current object.
 411        """
 412        return copy.deepcopy(self)
 413
 414    def keys(self) -> List[str]:
 415        """
 416        Retrieve a list of public attribute names for `self`.
 417
 418        Returns:
 419            **list of str:** 
 420                A list of attribute names, with internal name-mangling removed, 
 421                for easier reference to public attributes in the instance.
 422        """
 423        return [attr.replace('_SNPObject__', '') for attr in vars(self)]
 424
 425    def allele_freq(
 426        self,
 427        sample_labels: Optional[Sequence[Any]] = None,
 428        ancestry: Optional[Union[str, int]] = None,
 429        laiobj: Optional["LocalAncestryObject"] = None,
 430        return_counts: bool = False,
 431        as_dataframe: bool = False,
 432    ) -> Any:
 433        """
 434        Compute per-SNP alternate allele frequencies from `calldata_gt`.
 435
 436        Args:
 437            sample_labels (sequence, optional):
 438                Population label per sample. If None, computes cohort-level frequencies.
 439            ancestry (str or int, optional):
 440                If provided, compute ancestry-masked frequencies using SNP-level LAI.
 441            laiobj (LocalAncestryObject, optional):
 442                Optional LAI object used when `self.calldata_lai` is not set.
 443            return_counts (bool, default=False):
 444                If True, also return called-allele counts with the same shape as frequencies.
 445            as_dataframe (bool, default=False):
 446                If True, return pandas DataFrame output.
 447
 448        Returns:
 449            Frequencies as a NumPy array (or DataFrame if `as_dataframe=True`).
 450            If `return_counts=True`, returns `(freq, counts)`.
 451        """
 452        if self.calldata_gt is None:
 453            raise ValueError("Genotype data `calldata_gt` is None.")
 454
 455        gt = np.asarray(self.calldata_gt)
 456        if gt.ndim not in (2, 3):
 457            raise ValueError("'calldata_gt' must be 2D or 3D array")
 458
 459        n_samples = gt.shape[1]
 460
 461        grouped_output = sample_labels is not None
 462        if sample_labels is None:
 463            labels = np.repeat("__all__", n_samples)
 464        else:
 465            labels = np.asarray(sample_labels)
 466            if labels.ndim != 1:
 467                labels = labels.ravel()
 468            if labels.shape[0] != n_samples:
 469                raise ValueError(
 470                    "'sample_labels' must have length equal to the number of samples in `calldata_gt`."
 471                )
 472
 473        calldata_lai = None
 474        if ancestry is not None:
 475            if self.calldata_lai is not None:
 476                calldata_lai = self.calldata_lai
 477            elif laiobj is not None:
 478                try:
 479                    converted_lai = laiobj.convert_to_snp_level(snpobject=self, lai_format="3D")
 480                    calldata_lai = getattr(converted_lai, "calldata_lai", None)
 481                except Exception:
 482                    calldata_lai = None
 483
 484            if calldata_lai is None:
 485                raise ValueError(
 486                    "Ancestry-specific masking requires SNP-level LAI "
 487                    "(provide a LocalAncestryObject via 'laiobj' or ensure 'self.calldata_lai' is set)."
 488                )
 489
 490        afs, counts, pops = aggregate_pop_allele_freq(
 491            calldata_gt=gt,
 492            sample_labels=labels,
 493            ancestry=ancestry,
 494            calldata_lai=calldata_lai,
 495        )
 496
 497        if grouped_output:
 498            freq_out = afs
 499            count_out = counts
 500            if as_dataframe:
 501                import pandas as pd
 502
 503                freq_out = pd.DataFrame(afs, columns=pops)
 504                count_out = pd.DataFrame(counts, columns=pops)
 505        else:
 506            freq_out = afs[:, 0]
 507            count_out = counts[:, 0]
 508            if as_dataframe:
 509                import pandas as pd
 510
 511                freq_out = pd.DataFrame({"allele_freq": freq_out})
 512                count_out = pd.DataFrame({"called_alleles": count_out})
 513
 514        if return_counts:
 515            return freq_out, count_out
 516        return freq_out
 517
 518    def sum_strands(self, inplace: bool = False) -> Optional['SNPObject']:
 519        """
 520        Sum paternal and maternal strands.
 521
 522        Args:
 523            inplace (bool, default=False): 
 524                If True, modifies `self` in place. If False, returns a new `SNPObject` with the variants 
 525                filtered. Default is False.
 526
 527        Returns:
 528            **Optional[SNPObject]:** 
 529                A new `SNPObject` with summed strands if `inplace=False`. 
 530                If `inplace=True`, modifies `self` in place and returns None.
 531        """
 532        if self.calldata_gt is None:
 533            warnings.warn("Genotype data `calldata_gt` is None.")
 534            return None if not inplace else self
 535
 536        if self.are_strands_summed:
 537            warnings.warn("Genotype data `calldata_gt` is already summed.")
 538            return self if inplace else self.copy()
 539        
 540        if inplace:
 541            self.calldata_gt = self.calldata_gt.sum(axis=2, dtype=np.int8)
 542            return self
 543        else:
 544            snpobj = self.copy()
 545            snpobj.calldata_gt = self.calldata_gt.sum(axis=2, dtype=np.int8)
 546            return snpobj
 547
 548    def filter_variants(
 549            self, 
 550            chrom: Optional[Union[str, Sequence[str], np.ndarray, None]] = None, 
 551            pos: Optional[Union[int, Sequence[int], np.ndarray, None]] = None, 
 552            indexes: Optional[Union[int, Sequence[int], np.ndarray, None]] = None, 
 553            include: bool = True, 
 554            inplace: bool = False
 555        ) -> Optional['SNPObject']:
 556        """
 557        Filter variants based on specified chromosome names, variant positions, or variant indexes.
 558
 559        This method updates the `calldata_gt`, `variants_ref`, `variants_alt`, 
 560        `variants_chrom`, `variants_filter_pass`, `variants_id`, `variants_pos`,  
 561        `variants_qual`, and `lai` attributes to include or exclude the specified variants. The filtering 
 562        criteria can be based on chromosome names, variant positions, or indexes. If multiple 
 563        criteria are provided, their union is used for filtering. The order of the variants is preserved.
 564        
 565        Negative indexes are supported and follow 
 566        [NumPy's indexing conventions](https://numpy.org/doc/stable/user/basics.indexing.html).
 567
 568        Args:
 569            chrom (str or array_like of str, optional): 
 570                Chromosome(s) to filter variants by. Can be a single chromosome as a string or a sequence 
 571                of chromosomes. If both `chrom` and `pos` are provided, they must either have matching lengths 
 572                (pairing each chromosome with a position) or `chrom` should be a single value that applies to 
 573                all positions in `pos`. Default is None. 
 574            pos (int or array_like of int, optional): 
 575                Position(s) to filter variants by. Can be a single position as an integer or a sequence of positions. 
 576                If `chrom` is also provided, `pos` should either match `chrom` in length or `chrom` should be a 
 577                single value. Default is None.
 578            indexes (int or array_like of int, optional): 
 579                Index(es) of the variants to include or exclude. Can be a single index or a sequence
 580                of indexes. Negative indexes are supported. Default is None.
 581            include (bool, default=True): 
 582                If True, includes only the specified variants. If False, excludes the specified
 583                variants. Default is True.
 584            inplace (bool, default=False): 
 585                If True, modifies `self` in place. If False, returns a new `SNPObject` with the variants 
 586                filtered. Default is False.
 587
 588        Returns:
 589            **Optional[SNPObject]:** 
 590                A new `SNPObject` with the specified variants filtered if `inplace=False`. 
 591                If `inplace=True`, modifies `self` in place and returns None.
 592        """
 593        if chrom is None and pos is None and indexes is None:
 594            raise ValueError("At least one of 'chrom', 'pos', or 'indexes' must be provided.")
 595
 596        n_snps = self.n_snps
 597
 598        # Convert inputs to arrays for consistency
 599        chrom = np.atleast_1d(chrom) if chrom is not None else None
 600        pos = np.atleast_1d(pos) if pos is not None else None
 601        indexes = np.atleast_1d(indexes) if indexes is not None else None
 602
 603        # Validate chrom and pos lengths if both are provided
 604        if chrom is not None and pos is not None:
 605            if len(chrom) != len(pos) and len(chrom) > 1:
 606                raise ValueError(
 607                    "When both 'chrom' and 'pos' are provided, they must either be of the same length "
 608                    "or 'chrom' must be a single value."
 609                )
 610
 611        # Create a mask for chromosome and position filtering
 612        mask_combined = np.zeros(n_snps, dtype=bool)
 613        if chrom is not None and pos is not None:
 614            if len(chrom) == 1:
 615                # Apply single chromosome to all positions in `pos`
 616                mask_combined = (self['variants_chrom'] == chrom[0]) & np.isin(self['variants_pos'], pos)
 617            else:
 618                # Vectorized pair matching for chrom and pos
 619                query_pairs = np.array(
 620                    list(zip(chrom, pos)),
 621                    dtype=[
 622                        ('chrom', self['variants_chrom'].dtype),
 623                        ('pos', self['variants_pos'].dtype)
 624                    ]
 625                )
 626                data_pairs = np.array(
 627                    list(zip(self['variants_chrom'], self['variants_pos'])),
 628                    dtype=[
 629                        ('chrom', self['variants_chrom'].dtype),
 630                        ('pos', self['variants_pos'].dtype)
 631                    ]
 632                )
 633                mask_combined = np.isin(data_pairs, query_pairs)
 634
 635        elif chrom is not None:
 636            # Only chromosome filtering
 637            mask_combined = np.isin(self['variants_chrom'], chrom)
 638        elif pos is not None:
 639            # Only position filtering
 640            mask_combined = np.isin(self['variants_pos'], pos)
 641
 642        # Create mask based on indexes if provided
 643        if indexes is not None:
 644            # Validate indexes, allowing negative indexes
 645            out_of_bounds_indexes = indexes[(indexes < -n_snps) | (indexes >= n_snps)]
 646            if out_of_bounds_indexes.size > 0:
 647                raise ValueError(f"One or more sample indexes are out of bounds.")
 648
 649            # Handle negative indexes and check for out-of-bounds indexes
 650            adjusted_indexes = np.mod(indexes, n_snps)
 651
 652            # Create mask for specified indexes
 653            mask_indexes = np.zeros(n_snps, dtype=bool)
 654            mask_indexes[adjusted_indexes] = True
 655
 656            # Combine with `chrom` and `pos` mask using logical OR (union of all specified criteria)
 657            mask_combined = mask_combined | mask_indexes
 658
 659        # Invert mask if `include` is False
 660        if not include:
 661            mask_combined = ~mask_combined
 662
 663        # Define keys to filter
 664        keys = [
 665            'calldata_gt', 'variants_ref', 'variants_alt', 'variants_chrom', 'variants_filter_pass', 
 666            'variants_id', 'variants_pos', 'variants_qual', 'calldata_lai'
 667        ]
 668
 669        # Apply filtering based on inplace parameter
 670        if inplace:
 671            for key in keys:
 672                if self[key] is not None:
 673                    if self[key].ndim > 1:
 674                        self[key] = np.asarray(self[key])[mask_combined, ...]
 675                    else:
 676                        self[key] = np.asarray(self[key])[mask_combined]
 677
 678            return None
 679        else:
 680            # Create A new `SNPObject` with filtered data
 681            snpobj = self.copy()
 682            for key in keys:
 683                if snpobj[key] is not None:
 684                    if snpobj[key].ndim > 1:
 685                        snpobj[key] = np.asarray(snpobj[key])[mask_combined, ...]
 686                    else:
 687                        snpobj[key] = np.asarray(snpobj[key])[mask_combined]
 688
 689            return snpobj
 690
 691    def filter_samples(
 692            self, 
 693            samples: Optional[Union[str, Sequence[str], np.ndarray, None]] = None,
 694            indexes: Optional[Union[int, Sequence[int], np.ndarray, None]] = None,
 695            include: bool = True,
 696            reorder: bool = False,
 697            inplace: bool = False
 698        ) -> Optional['SNPObject']:
 699        """
 700        Filter samples based on specified names or indexes.
 701
 702        This method updates the `samples` and `calldata_gt` attributes to include or exclude the specified 
 703        samples. The order of the samples is preserved. Set `reorder=True` to match the ordering of the
 704        provided `samples` and/or `indexes` lists when including.
 705
 706        If both samples and indexes are provided, any sample matching either a name in samples or an index in 
 707        indexes will be included or excluded.
 708
 709        This method allows inclusion or exclusion of specific samples by their names or 
 710        indexes. When both sample names and indexes are provided, the union of the specified samples 
 711        is used. Negative indexes are supported and follow 
 712        [NumPy's indexing conventions](https://numpy.org/doc/stable/user/basics.indexing.html).
 713
 714        Args:
 715            samples (str or array_like of str, optional): 
 716                 Name(s) of the samples to include or exclude. Can be a single sample name or a
 717                 sequence of sample names. Default is None.
 718            indexes (int or array_like of int, optional):
 719                Index(es) of the samples to include or exclude. Can be a single index or a sequence
 720                of indexes. Negative indexes are supported. Default is None.
 721            include (bool, default=True): 
 722                If True, includes only the specified samples. If False, excludes the specified
 723                samples. Default is True.
 724            inplace (bool, default=False): 
 725                If True, modifies `self` in place. If False, returns a new `SNPObject` with the samples 
 726                filtered. Default is False.
 727
 728        Returns:
 729            **Optional[SNPObject]:** 
 730                A new `SNPObject` with the specified samples filtered if `inplace=False`. 
 731                If `inplace=True`, modifies `self` in place and returns None.
 732        """
 733        if samples is None and indexes is None:
 734            raise ValueError("At least one of 'samples' or 'indexes' must be provided.")
 735
 736        n_samples = self.n_samples
 737        sample_names = np.array(self['samples'])
 738
 739        # Create mask based on sample names
 740        if samples is not None:
 741            samples = np.asarray(samples).ravel()
 742            mask_samples = np.isin(sample_names, samples)
 743            missing_samples = samples[~np.isin(samples, sample_names)]
 744            if missing_samples.size > 0:
 745                raise ValueError(f"The following specified samples were not found: {missing_samples.tolist()}")
 746        else:
 747            mask_samples = np.zeros(n_samples, dtype=bool)
 748
 749        # Create mask based on sample indexes
 750        if indexes is not None:
 751            indexes = np.asarray(indexes).ravel()
 752
 753            # Validate indexes, allowing negative indexes
 754            out_of_bounds_indexes = indexes[(indexes < -n_samples) | (indexes >= n_samples)]
 755            if out_of_bounds_indexes.size > 0:
 756                raise ValueError(f"One or more sample indexes are out of bounds.")
 757            
 758            # Handle negative indexes
 759            adjusted_indexes = np.mod(indexes, n_samples)
 760
 761            mask_indexes = np.zeros(n_samples, dtype=bool)
 762            mask_indexes[adjusted_indexes] = True
 763        else:
 764            mask_indexes = np.zeros(n_samples, dtype=bool)
 765
 766        # Combine masks using logical OR (union of samples)
 767        mask_combined = mask_samples | mask_indexes
 768
 769        if not include:
 770            mask_combined = ~mask_combined
 771
 772        # If requested, compute an ordering of selected samples that follows the provided lists.
 773        ordered_indices = None
 774        if include and reorder:
 775            sel_indices = np.where(mask_combined)[0]
 776            ordered_list: List[int] = []
 777            added = np.zeros(n_samples, dtype=bool)
 778
 779            # Prioritize the order in `samples`
 780            if samples is not None:
 781                name_to_idx = {name: idx for idx, name in enumerate(sample_names)}
 782                for s in samples:
 783                    idx = name_to_idx.get(s)
 784                    if idx is not None and mask_combined[idx] and not added[idx]:
 785                        ordered_list.append(idx)
 786                        added[idx] = True
 787
 788            # Then respect the order in `indexes`
 789            if indexes is not None:
 790                adj_idx = np.mod(np.atleast_1d(indexes), n_samples)
 791                for idx in adj_idx:
 792                    if mask_combined[idx] and not added[idx]:
 793                        ordered_list.append(int(idx))
 794                        added[idx] = True
 795
 796            # Finally, append any remaining selected samples in their original order
 797            for idx in sel_indices:
 798                if not added[idx]:
 799                    ordered_list.append(int(idx))
 800
 801            ordered_indices = np.asarray(ordered_list, dtype=int)
 802
 803        # Define keys to filter
 804        keys = ['samples', 'calldata_gt', 'calldata_lai']
 805
 806        # Apply filtering based on inplace parameter
 807        if inplace:
 808            for key in keys:
 809                if self[key] is not None:
 810                    arr = np.asarray(self[key])
 811                    if ordered_indices is not None:
 812                        if key == 'calldata_lai' and arr.ndim == 2:
 813                            # Haplotype-aware reordering for 2D LAI (n_snps, n_samples*2)
 814                            hap_idx = np.concatenate([2*ordered_indices, 2*ordered_indices + 1])
 815                            self[key] = arr[:, hap_idx]
 816                        elif arr.ndim > 1:
 817                            self[key] = arr[:, ordered_indices, ...]
 818                        else:
 819                            self[key] = arr[ordered_indices]
 820                    else:
 821                        if arr.ndim > 1:
 822                            self[key] = arr[:, mask_combined, ...]
 823                        else:
 824                            self[key] = arr[mask_combined]
 825
 826            return None
 827        else:
 828            # Create A new `SNPObject` with filtered data
 829            snpobj = self.copy()
 830            for key in keys:
 831                if snpobj[key] is not None:
 832                    arr = np.asarray(snpobj[key])
 833                    if ordered_indices is not None:
 834                        if key == 'calldata_lai' and arr.ndim == 2:
 835                            hap_idx = np.concatenate([2*ordered_indices, 2*ordered_indices + 1])
 836                            snpobj[key] = arr[:, hap_idx]
 837                        elif arr.ndim > 1:
 838                            snpobj[key] = arr[:, ordered_indices, ...]
 839                        else:
 840                            snpobj[key] = arr[ordered_indices]
 841                    else:
 842                        if arr.ndim > 1:
 843                            snpobj[key] = arr[:, mask_combined, ...]
 844                        else:
 845                            snpobj[key] = arr[mask_combined]
 846            return snpobj
 847
 848    def detect_chromosome_format(self) -> str:
 849        """
 850        Detect the chromosome naming convention in `variants_chrom` based on the prefix 
 851        of the first chromosome identifier in `unique_chrom`.
 852        
 853        **Recognized formats:**
 854
 855        - `'chr'`: Format with 'chr' prefix, e.g., 'chr1', 'chr2', ..., 'chrX', 'chrY', 'chrM'.
 856        - `'chm'`: Format with 'chm' prefix, e.g., 'chm1', 'chm2', ..., 'chmX', 'chmY', 'chmM'.
 857        - `'chrom'`: Format with 'chrom' prefix, e.g., 'chrom1', 'chrom2', ..., 'chromX', 'chromY', 'chromM'.
 858        - `'plain'`: Plain format without a prefix, e.g., '1', '2', ..., 'X', 'Y', 'M'.
 859        
 860        If the format does not match any recognized pattern, `'Unknown format'` is returned.
 861
 862        Returns:
 863            **str:** 
 864                A string indicating the detected chromosome format (`'chr'`, `'chm'`, `'chrom'`, or `'plain'`).
 865                If no recognized format is matched, returns `'Unknown format'`.
 866        """
 867        # Select the first unique chromosome identifier for format detection
 868        chromosome_str = self.unique_chrom[0]
 869
 870        # Define regular expressions to match each recognized chromosome format
 871        patterns = {
 872            'chr': r'^chr(\d+|X|Y|M)$',    # Matches 'chr' prefixed format
 873            'chm': r'^chm(\d+|X|Y|M)$',    # Matches 'chm' prefixed format
 874            'chrom': r'^chrom(\d+|X|Y|M)$', # Matches 'chrom' prefixed format
 875            'plain': r'^(\d+|X|Y|M)$'       # Matches plain format without prefix
 876        }
 877
 878        # Iterate through the patterns to identify the chromosome format
 879        for prefix, pattern in patterns.items():
 880            if re.match(pattern, chromosome_str):
 881                return prefix  # Return the recognized format prefix
 882
 883        # If no pattern matches, return 'Unknown format'
 884        return 'Unknown format'
 885
 886    def convert_chromosome_format(
 887        self, 
 888        from_format: str, 
 889        to_format: str, 
 890        inplace: bool = False
 891    ) -> Optional['SNPObject']:
 892        """
 893        Convert the chromosome format from one naming convention to another in `variants_chrom`.
 894
 895        **Supported formats:**
 896
 897        - `'chr'`: Format with 'chr' prefix, e.g., 'chr1', 'chr2', ..., 'chrX', 'chrY', 'chrM'.
 898        - `'chm'`: Format with 'chm' prefix, e.g., 'chm1', 'chm2', ..., 'chmX', 'chmY', 'chmM'.
 899        - `'chrom'`: Format with 'chrom' prefix, e.g., 'chrom1', 'chrom2', ..., 'chromX', 'chromY', 'chromM'.
 900        - `'plain'`: Plain format without a prefix, e.g., '1', '2', ..., 'X', 'Y', 'M'.
 901
 902        Args:
 903            from_format (str): 
 904                The current chromosome format. Acceptable values are `'chr'`, `'chm'`, `'chrom'`, or `'plain'`.
 905            to_format (str): 
 906                The target format for chromosome data conversion. Acceptable values match `from_format` options.
 907            inplace (bool, default=False): 
 908                If True, modifies `self` in place. If False, returns a new `SNPObject` with the converted format.
 909                Default is False.
 910
 911        Returns:
 912            **Optional[SNPObject]:** A new `SNPObject` with the converted chromosome format if `inplace=False`. 
 913            If `inplace=True`, modifies `self` in place and returns None.
 914        """
 915        # Define the list of standard chromosome identifiers
 916        chrom_list = [*map(str, range(1, 23)), 'X', 'Y', 'M']  # M for mitochondrial chromosomes
 917
 918        # Format mappings for different chromosome naming conventions
 919        format_mappings = {
 920            'chr': [f'chr{i}' for i in chrom_list],
 921            'chm': [f'chm{i}' for i in chrom_list],
 922            'chrom': [f'chrom{i}' for i in chrom_list],
 923            'plain': chrom_list,
 924        }
 925
 926        # Verify that from_format and to_format are valid naming conventions
 927        if from_format not in format_mappings or to_format not in format_mappings:
 928            raise ValueError(f"Invalid format: {from_format} or {to_format}. Must be one of {list(format_mappings.keys())}.")
 929
 930        # Convert chromosomes to string for consistent comparison
 931        variants_chrom = self['variants_chrom'].astype(str)
 932
 933        # Verify that all chromosomes in the object follow the specified `from_format`
 934        expected_chroms = set(format_mappings[from_format])
 935        mismatched_chroms = set(variants_chrom) - expected_chroms
 936
 937        if mismatched_chroms:
 938            raise ValueError(f"The following chromosomes do not match the `from_format` '{from_format}': {mismatched_chroms}.")
 939
 940        # Create conditions for selecting based on current `from_format` names
 941        conditions = [variants_chrom == chrom for chrom in format_mappings[from_format]]
 942
 943        # Rename chromosomes based on inplace flag
 944        if inplace:
 945            self['variants_chrom'] = np.select(conditions, format_mappings[to_format], default='unknown')
 946            return None
 947        else:
 948            snpobject = self.copy()
 949            snpobject['variants_chrom'] = np.select(conditions, format_mappings[to_format], default='unknown')
 950            return snpobject
 951
 952    def match_chromosome_format(self, snpobj: 'SNPObject', inplace: bool = False) -> Optional['SNPObject']:
 953        """
 954        Convert the chromosome format in `variants_chrom` from `self` to match the format of a reference `snpobj`.
 955
 956        **Recognized formats:**
 957
 958        - `'chr'`: Format with 'chr' prefix, e.g., 'chr1', 'chr2', ..., 'chrX', 'chrY', 'chrM'.
 959        - `'chm'`: Format with 'chm' prefix, e.g., 'chm1', 'chm2', ..., 'chmX', 'chmY', 'chmM'.
 960        - `'chrom'`: Format with 'chrom' prefix, e.g., 'chrom1', 'chrom2', ..., 'chromX', 'chromY', 'chromM'.
 961        - `'plain'`: Plain format without a prefix, e.g., '1', '2', ..., 'X', 'Y', 'M'.
 962
 963        Args:
 964            snpobj (SNPObject): 
 965                The reference SNPObject whose chromosome format will be matched.
 966            inplace (bool, default=False): 
 967                If True, modifies `self` in place. If False, returns a new `SNPObject` with the 
 968                chromosome format matching that of `snpobj`. Default is False.
 969
 970        Returns:
 971            **Optional[SNPObject]:** 
 972                A new `SNPObject` with matched chromosome format if `inplace=False`. 
 973                If `inplace=True`, modifies `self` in place and returns None.
 974        """
 975        # Detect the chromosome naming format of the current SNPObject
 976        fmt1 = self.detect_chromosome_format()
 977        if fmt1 == 'Unknown format':
 978            raise ValueError("The chromosome format of the current SNPObject is unrecognized.")
 979        
 980        # Detect the chromosome naming format of the reference SNPObject
 981        fmt2 = snpobj.detect_chromosome_format()
 982        if fmt2 == 'Unknown format':
 983            raise ValueError("The chromosome format of the reference SNPObject is unrecognized.")
 984
 985        # Convert the current SNPObject's chromosome format to match the reference format
 986        return self.convert_chromosome_format(fmt1, fmt2, inplace=inplace)
 987
 988    def rename_chrom(
 989        self,
 990        to_replace: Union[Dict[str, str], str, List[str]] = {'^([0-9]+)$': r'chr\1', r'^chr([0-9]+)$': r'\1'},
 991        value: Optional[Union[str, List[str]]] = None,
 992        regex: bool = True,
 993        inplace: bool = False
 994    ) -> Optional['SNPObject']:
 995        """
 996        Replace chromosome values in `variants_chrom` using patterns or exact matches.
 997
 998        This method allows flexible chromosome replacements, using regex or exact matches, useful 
 999        for non-standard chromosome formats. For standard conversions (e.g., 'chr1' to '1'), 
1000        consider `convert_chromosome_format`.
1001
1002        Args:
1003            to_replace (dict, str, or list of str): 
1004                Pattern(s) or exact value(s) to be replaced in chromosome names. Default behavior 
1005                transforms `<chrom_num>` to `chr<chrom_num>` or vice versa. Non-matching values 
1006                remain unchanged.
1007                - If str or list of str: Matches will be replaced with `value`.
1008                - If regex (bool), then any regex matches will be replaced with `value`.
1009                - If dict: Keys defines values to replace, with corresponding replacements as values.
1010            value (str or list of str, optional): 
1011                Replacement value(s) if `to_replace` is a string or list. Ignored if `to_replace` 
1012                is a dictionary.
1013            regex (bool, default=True): 
1014                If True, interprets `to_replace` keys as regex patterns.
1015            inplace (bool, default=False): 
1016                If True, modifies `self` in place. If False, returns a new `SNPObject` with the chromosomes
1017                renamed. Default is False.
1018
1019        Returns:
1020            **Optional[SNPObject]:** A new `SNPObject` with the renamed chromosome format if `inplace=False`. 
1021            If `inplace=True`, modifies `self` in place and returns None.
1022        """
1023        # Standardize input format: convert `to_replace` and `value` to a dictionary if needed
1024        if isinstance(to_replace, (str, int)):
1025            to_replace = [to_replace]
1026        if isinstance(value, (str, int)):
1027            value = [value]
1028        if isinstance(to_replace, list) and isinstance(value, list):
1029            dictionary = dict(zip(to_replace, value))
1030        elif isinstance(to_replace, dict) and value is None:
1031            dictionary = to_replace
1032        else:
1033            raise ValueError(
1034            "Invalid input: `to_replace` and `value` must be compatible types (both str, list of str, or dict)."
1035        )
1036
1037        # Vectorized function for replacing values in chromosome array
1038        vec_replace_values = np.vectorize(self._match_to_replace)
1039
1040        # Rename chromosomes based on inplace flag
1041        if inplace:
1042            self.variants_chrom = vec_replace_values(self.variants_chrom, dictionary, regex)
1043            return None
1044        else:
1045            snpobj = self.copy()
1046            snpobj.variants_chrom = vec_replace_values(self.variants_chrom, dictionary, regex)
1047            return snpobj
1048
1049    def rename_missings(
1050        self, 
1051        before: Union[int, float, str] = -1, 
1052        after: Union[int, float, str] = '.', 
1053        inplace: bool = False
1054    ) -> Optional['SNPObject']:
1055        """
1056        Replace missing values in the `calldata_gt` attribute.
1057
1058        This method identifies missing values in 'calldata_gt' and replaces them with a specified 
1059        value. By default, it replaces occurrences of `-1` (often used to signify missing data) with `'.'`.
1060
1061        Args:
1062            before (int, float, or str, default=-1): 
1063                The current representation of missing values in `calldata_gt`. Common values might be -1, '.', or NaN.
1064                Default is -1.
1065            after (int, float, or str, default='.'): 
1066                The value that will replace `before`. Default is '.'.
1067            inplace (bool, default=False): 
1068                If True, modifies `self` in place. If False, returns a new `SNPObject` with the applied 
1069                replacements. Default is False.
1070
1071        Returns:
1072            **Optional[SNPObject]:** 
1073                A new `SNPObject` with the renamed missing values if `inplace=False`. 
1074                If `inplace=True`, modifies `self` in place and returns None.
1075        """
1076        # Rename missing values in the `calldata_gt` attribute based on inplace flag
1077        if inplace:
1078            self['calldata_gt'] = np.where(self['calldata_gt'] == before, after, self['calldata_gt'])
1079            return None
1080        else:
1081            snpobj = self.copy()
1082            snpobj['calldata_gt'] = np.where(snpobj['calldata_gt'] == before, after, snpobj['calldata_gt'])
1083            return snpobj
1084
1085    def get_common_variants_intersection(
1086        self, 
1087        snpobj: 'SNPObject', 
1088        index_by: str = 'pos'
1089    ) -> Tuple[List[str], np.ndarray, np.ndarray]:
1090        """
1091        Identify common variants between `self` and the `snpobj` instance based on the specified `index_by` criterion, 
1092        which may match based on chromosome and position (`variants_chrom`, `variants_pos`), ID (`variants_id`), or both.
1093
1094        This method returns the identifiers of common variants and their corresponding indices in both objects.
1095
1096        Args:
1097            snpobj (SNPObject): 
1098                The reference SNPObject to compare against.
1099            index_by (str, default='pos'): 
1100                Criteria for matching variants. Options:
1101                - `'pos'`: Matches by chromosome and position (`variants_chrom`, `variants_pos`), e.g., 'chr1-12345'.
1102                - `'id'`: Matches by variant ID alone (`variants_id`), e.g., 'rs123'.
1103                - `'pos+id'`: Matches by chromosome, position, and ID (`variants_chrom`, `variants_pos`, `variants_id`), e.g., 'chr1-12345-rs123'.
1104                Default is 'pos'.
1105
1106        Returns:
1107            Tuple containing:
1108            - **list of str:** A list of common variant identifiers (as strings).
1109            - **array:** An array of indices in `self` where common variants are located.
1110            - **array:** An array of indices in `snpobj` where common variants are located.
1111        """
1112        # Create unique identifiers for each variant in both SNPObjects based on the specified criterion
1113        if index_by == 'pos':
1114            query_identifiers = [f"{chrom}-{pos}" for chrom, pos in zip(self['variants_chrom'], self['variants_pos'])]
1115            reference_identifiers = [f"{chrom}-{pos}" for chrom, pos in zip(snpobj['variants_chrom'], snpobj['variants_pos'])]
1116        elif index_by == 'id':
1117            query_identifiers = self['variants_id'].tolist()
1118            reference_identifiers = snpobj['variants_id'].tolist()
1119        elif index_by == 'pos+id':
1120            query_identifiers = [
1121                f"{chrom}-{pos}-{ids}" for chrom, pos, ids in zip(self['variants_chrom'], self['variants_pos'], self['variants_id'])
1122            ]
1123            reference_identifiers = [
1124                f"{chrom}-{pos}-{ids}" for chrom, pos, ids in zip(snpobj['variants_chrom'], snpobj['variants_pos'], snpobj['variants_id'])
1125            ]
1126        else:
1127            raise ValueError("`index_by` must be one of 'pos', 'id', or 'pos+id'.")
1128
1129        # Convert to sets for intersection
1130        common_ids = set(query_identifiers).intersection(reference_identifiers)
1131
1132        # Collect indices for common identifiers
1133        query_idx = [i for i, id in enumerate(query_identifiers) if id in common_ids]
1134        reference_idx = [i for i, id in enumerate(reference_identifiers) if id in common_ids]
1135
1136        return list(common_ids), np.array(query_idx), np.array(reference_idx)
1137
1138    def get_common_markers_intersection(
1139        self, 
1140        snpobj: 'SNPObject'
1141    ) -> Tuple[List[str], np.ndarray, np.ndarray]:
1142        """
1143        Identify common markers between between `self` and the `snpobj` instance. Common markers are identified 
1144        based on matching chromosome (`variants_chrom`), position (`variants_pos`), reference (`variants_ref`), 
1145        and alternate (`variants_alt`) alleles.
1146
1147        This method returns the identifiers of common markers and their corresponding indices in both objects.
1148        
1149        Args:
1150            snpobj (SNPObject): 
1151                The reference SNPObject to compare against.
1152        
1153        Returns:
1154            Tuple containing:
1155            - **list of str:** A list of common variant identifiers (as strings).
1156            - **array:** An array of indices in `self` where common variants are located.
1157            - **array:** An array of indices in `snpobj` where common variants are located.
1158        """
1159        # Generate unique identifiers based on chrom, pos, ref, and alt alleles
1160        query_identifiers = [
1161            f"{chrom}-{pos}-{ref}-{alt}" for chrom, pos, ref, alt in 
1162            zip(self['variants_chrom'], self['variants_pos'], self['variants_ref'], self['variants_alt'])
1163        ]
1164        reference_identifiers = [
1165            f"{chrom}-{pos}-{ref}-{alt}" for chrom, pos, ref, alt in 
1166            zip(snpobj['variants_chrom'], snpobj['variants_pos'], snpobj['variants_ref'], snpobj['variants_alt'])
1167        ]
1168
1169        # Convert to sets for intersection
1170        common_ids = set(query_identifiers).intersection(reference_identifiers)
1171
1172        # Collect indices for common identifiers in both SNPObjects
1173        query_idx = [i for i, id in enumerate(query_identifiers) if id in common_ids]
1174        reference_idx = [i for i, id in enumerate(reference_identifiers) if id in common_ids]
1175
1176        return list(common_ids), np.array(query_idx), np.array(reference_idx)
1177
1178    def subset_to_common_variants(
1179        self, 
1180        snpobj: 'SNPObject', 
1181        index_by: str = 'pos', 
1182        common_variants_intersection: Optional[Tuple[np.ndarray, np.ndarray]] = None,
1183        inplace: bool = False
1184    ) -> Optional['SNPObject']:
1185        """
1186        Subset `self` to include only the common variants with a reference `snpobj` based on 
1187        the specified `index_by` criterion, which may match based on chromosome and position 
1188        (`variants_chrom`, `variants_pos`), ID (`variants_id`), or both.
1189        
1190        Args:
1191            snpobj (SNPObject): 
1192                The reference SNPObject to compare against.
1193            index_by (str, default='pos'): 
1194                Criteria for matching variants. Options:
1195                - `'pos'`: Matches by chromosome and position (`variants_chrom`, `variants_pos`), e.g., 'chr1-12345'.
1196                - `'id'`: Matches by variant ID alone (`variants_id`), e.g., 'rs123'.
1197                - `'pos+id'`: Matches by chromosome, position, and ID (`variants_chrom`, `variants_pos`, `variants_id`), e.g., 'chr1-12345-rs123'.
1198                Default is 'pos'.
1199            common_variants_intersection (Tuple[np.ndarray, np.ndarray], optional): 
1200                Precomputed indices of common variants between `self` and `snpobj`. If None, intersection is 
1201                computed within the function.
1202            inplace (bool, default=False): 
1203                If True, modifies `self` in place. If False, returns a new `SNPObject` with the common variants
1204                subsetted. Default is False.
1205
1206        Returns:
1207            **Optional[SNPObject]:** 
1208                A new `SNPObject` with the common variants subsetted if `inplace=False`. 
1209                If `inplace=True`, modifies `self` in place and returns None.
1210        """
1211        # Get indices of common variants if not provided
1212        if common_variants_intersection is None:
1213            _, query_idx, _ = self.get_common_variants_intersection(snpobj, index_by=index_by)
1214        else:
1215            query_idx, _ = common_variants_intersection
1216
1217        # Use filter_variants method with the identified indices, applying `inplace` as specified
1218        return self.filter_variants(indexes=query_idx, include=True, inplace=inplace)
1219
1220    def subset_to_common_markers(
1221        self, 
1222        snpobj: 'SNPObject', 
1223        common_markers_intersection: Optional[Tuple[np.ndarray, np.ndarray]] = None,
1224        inplace: bool = False
1225    ) -> Optional['SNPObject']:
1226        """
1227        Subset `self` to include only the common markers with a reference `snpobj`. Common markers are identified 
1228        based on matching chromosome (`variants_chrom`), position (`variants_pos`), reference (`variants_ref`), 
1229        and alternate (`variants_alt`) alleles.
1230
1231        Args:
1232            snpobj (SNPObject): 
1233                The reference SNPObject to compare against.
1234            common_markers_intersection (tuple of arrays, optional): 
1235                Precomputed indices of common markers between `self` and `snpobj`. If None, intersection is 
1236                computed within the function.
1237            inplace (bool, default=False): 
1238                If True, modifies `self` in place. If False, returns a new `SNPObject` with the common markers
1239                subsetted. Default is False.
1240
1241        Returns:
1242            **Optional[SNPObject]:** 
1243                A new `SNPObject` with the common markers subsetted if `inplace=False`. 
1244                If `inplace=True`, modifies `self` in place and returns None.
1245        """
1246        # Get indices of common markers if not provided
1247        if common_markers_intersection is None:
1248            _, query_idx, _ = self.get_common_markers_intersection(snpobj)
1249        else:
1250            query_idx, _ = common_markers_intersection
1251
1252        # Use filter_variants method with the identified indices, applying `inplace` as specified
1253        return self.filter_variants(indexes=query_idx, include=True, inplace=inplace)
1254
1255    def merge(
1256            self, 
1257            snpobj: 'SNPObject', 
1258            force_samples: bool = False, 
1259            prefix: str = '2', 
1260            inplace: bool = False
1261        ) -> Optional['SNPObject']:
1262        """
1263        Merge `self` with `snpobj` along the sample axis.
1264
1265        This method expects both SNPObjects to contain the same set of SNPs in the same order, 
1266        then combines their genotype (`calldata_gt`) and LAI (`calldata_lai`) arrays by 
1267        concatenating the sample dimension. Samples from `snpobj` are appended to those in `self`.
1268
1269        Args:
1270            snpobj (SNPObject): 
1271                The SNPObject to merge samples with.
1272            force_samples (bool, default=False): 
1273                If True, duplicate sample names are resolved by prepending the `prefix` to duplicate sample names in 
1274                `snpobj`. Otherwise, merging fails when duplicate sample names are found. Default is False.
1275            prefix (str, default='2'): 
1276                A string prepended to duplicate sample names in `snpobj` when `force_samples=True`. 
1277                Duplicates are renamed from `<sample_name>` to `<prefix>:<sample_name>`. For instance, 
1278                if `prefix='2'` and there is a conflict with a sample called "sample_1", it becomes "2:sample_1".
1279            inplace (bool, default=False): 
1280                If True, modifies `self` in place. If False, returns a new `SNPObject` with the merged samples. 
1281                Default is False.
1282
1283        Returns:
1284            **Optional[SNPObject]**: A new SNPObject containing the merged sample data.
1285        """
1286        # Merge calldata_gt if present and compatible
1287        if self.calldata_gt is not None and snpobj.calldata_gt is not None:
1288            if self.calldata_gt.shape[0] != snpobj.calldata_gt.shape[0]:
1289                raise ValueError(
1290                    f"Cannot merge SNPObjects: Mismatch in the number of SNPs in `calldata_gt`.\n"
1291                    f"`self.calldata_gt` has {self.calldata_gt.shape[0]} SNPs, "
1292                    f"while `snpobj.calldata_gt` has {snpobj.calldata_gt.shape[0]} SNPs."
1293                )
1294            if self.are_strands_summed and not snpobj.are_strands_summed:
1295                raise ValueError(
1296                    "Cannot merge SNPObjects: `self` has summed strands, but `snpobj` does not.\n"
1297                    "Ensure both objects have the same genotype summation state before merging."
1298                )
1299            if not self.are_strands_summed and snpobj.are_strands_summed:
1300                raise ValueError(
1301                    "Cannot merge SNPObjects: `snpobj` has summed strands, but `self` does not.\n"
1302                    "Ensure both objects have the same genotype summation state before merging."
1303                )
1304            calldata_gt = np.concatenate([self.calldata_gt, snpobj.calldata_gt], axis=1)
1305        else:
1306            calldata_gt = None
1307
1308        # Merge samples if present and compatible, handling duplicates if `force_samples=True`
1309        if self.samples is not None and snpobj.samples is not None:
1310            overlapping_samples = set(self.samples).intersection(set(snpobj.samples))
1311            if overlapping_samples:
1312                if not force_samples:
1313                    raise ValueError(
1314                        f"Cannot merge SNPObjects: Found overlapping sample names {overlapping_samples}.\n"
1315                        "Samples must be strictly non-overlapping. To allow merging with renaming, set `force_samples=True`."
1316                    )
1317                else:
1318                    # Rename duplicate samples by prepending the file index
1319                    renamed_samples = [f"{prefix}:{sample}" if sample in overlapping_samples else sample for sample in snpobj.samples]
1320                    samples = np.concatenate([self.samples, renamed_samples], axis=0)
1321            else:
1322                samples = np.concatenate([self.samples, snpobj.samples], axis=0)
1323        else:
1324            samples = None
1325
1326        # Merge LAI data if present and compatible
1327        if self.calldata_lai is not None and snpobj.calldata_lai is not None:
1328            if self.calldata_lai.ndim != snpobj.calldata_lai.ndim:
1329                raise ValueError(
1330                    f"Cannot merge SNPObjects: Mismatch in `calldata_lai` dimensions.\n"
1331                    f"`self.calldata_lai` has {self.calldata_lai.ndim} dimensions, "
1332                    f"while `snpobj.calldata_lai` has {snpobj.calldata_lai.ndim} dimensions."
1333                )
1334            if self.calldata_lai.shape[0] != snpobj.calldata_lai.shape[0]:
1335                raise ValueError(
1336                    f"Cannot merge SNPObjects: Mismatch in the number of SNPs in `calldata_lai`.\n"
1337                    f"`self.calldata_lai` has {self.calldata_lai.shape[0]} SNPs, "
1338                    f"while `snpobj.calldata_lai` has {snpobj.calldata_lai.shape[0]} SNPs."
1339                )
1340            calldata_lai = np.concatenate([self.calldata_lai, snpobj.calldata_lai], axis=1)
1341        else:
1342            calldata_lai = None
1343
1344        if inplace:
1345            self.calldata_gt = calldata_gt
1346            self.calldata_lai = calldata_lai
1347            self.samples = samples
1348            return self
1349
1350        # Create and return a new SNPObject containing the merged samples
1351        return SNPObject(
1352            calldata_gt=calldata_gt,
1353            samples=samples,
1354            variants_ref=self.variants_ref,
1355            variants_alt=self.variants_alt,
1356            variants_chrom=self.variants_chrom,
1357            variants_filter_pass=self.variants_filter_pass,
1358            variants_id=self.variants_id,
1359            variants_pos=self.variants_pos,
1360            variants_qual=self.variants_qual,
1361            calldata_lai=calldata_lai,
1362            ancestry_map=self.ancestry_map
1363        )
1364    
1365    def concat(
1366        self,
1367        snpobj: 'SNPObject', 
1368        inplace: bool = False
1369    ) -> Optional['SNPObject']:
1370        """
1371        Concatenate self with snpobj along the SNP axis.
1372
1373        This method expects both SNPObjects to contain the same set of samples in the same order, 
1374        and that the chromosome(s) in snpobj follow (i.e. have higher numeric identifiers than) 
1375        those in self.
1376
1377        Args:
1378            snpobj (SNPObject):
1379                The SNPObject to concatenate SNPs with.
1380            inplace (bool, default=False):
1381                If True, modifies `self` in place. If False, returns a new `SNPObject` with the concatenated SNPs. 
1382                Default is False.
1383        
1384        Returns:
1385            **Optional[SNPObject]**: A new SNPObject containing the concatenated SNP data.
1386        """
1387        # Merge calldata_gt if present and compatible
1388        if self.calldata_gt is not None and snpobj.calldata_gt is not None:
1389            if self.calldata_gt.shape[1] != snpobj.calldata_gt.shape[1]:
1390                raise ValueError(
1391                    f"Cannot merge SNPObjects: Mismatch in the number of samples in `calldata_gt`.\n"
1392                    f"`self.calldata_gt` has {self.calldata_gt.shape[1]} samples, "
1393                    f"while `snpobj.calldata_gt` has {snpobj.calldata_gt.shape[1]} samples."
1394                )
1395            if self.are_strands_summed and not snpobj.are_strands_summed:
1396                raise ValueError(
1397                    "Cannot merge SNPObjects: `self` has summed strands, but `snpobj` does not.\n"
1398                    "Ensure both objects have the same genotype summation state before merging."
1399                )
1400            if not self.are_strands_summed and snpobj.are_strands_summed:
1401                raise ValueError(
1402                    "Cannot merge SNPObjects: `snpobj` has summed strands, but `self` does not.\n"
1403                    "Ensure both objects have the same genotype summation state before merging."
1404                )
1405            calldata_gt = np.concatenate([self.calldata_gt, snpobj.calldata_gt], axis=0)
1406        else:
1407            calldata_gt = None
1408
1409        # Merge SNP-related attributes if present
1410        attributes = [
1411            'variants_ref', 'variants_alt', 'variants_chrom', 'variants_filter_pass', 'variants_id', 'variants_pos', 'variants_qual'
1412        ]
1413        merged_attrs = {}
1414        for attr in attributes:
1415            self_attr = getattr(self, attr, None)
1416            obj_attr = getattr(snpobj, attr, None)
1417
1418            # Concatenate if both present
1419            if self_attr is not None and obj_attr is not None:
1420                merged_attrs[attr] = np.concatenate([self_attr, obj_attr], axis=0)
1421            else:
1422                # If either is None, store None
1423                merged_attrs[attr] = None
1424
1425        # Merge LAI data if present and compatible
1426        if self.calldata_lai is not None and snpobj.calldata_lai is not None:
1427            if self.calldata_lai.ndim != snpobj.calldata_lai.ndim:
1428                raise ValueError(
1429                    f"Cannot merge SNPObjects: Mismatch in `calldata_lai` dimensions.\n"
1430                    f"`self.calldata_lai` has {self.calldata_lai.ndim} dimensions, "
1431                    f"while `snpobj.calldata_lai` has {snpobj.calldata_lai.ndim} dimensions."
1432                )
1433            if self.calldata_lai.shape[1] != snpobj.calldata_lai.shape[1]:
1434                raise ValueError(
1435                    f"Cannot merge SNPObjects: Mismatch in the number of samples in `calldata_lai`.\n"
1436                    f"`self.calldata_lai` has {self.calldata_lai.shape[1]} samples, "
1437                    f"while `snpobj.calldata_lai` has {snpobj.calldata_lai.shape[1]} samples."
1438                )
1439            calldata_lai = np.concatenate([self.calldata_lai, snpobj.calldata_lai], axis=0)
1440        else:
1441            calldata_lai = None
1442        
1443        if inplace:
1444            self.calldata_gt = calldata_gt
1445            self.calldata_lai = calldata_lai
1446            for attr in attributes:
1447                self[attr] = merged_attrs[attr]
1448            return self
1449        
1450        # Create and return a new SNPObject containing the concatenated SNPs
1451        return SNPObject(
1452            calldata_gt=calldata_gt,
1453            calldata_lai=calldata_lai,
1454            samples=self.samples,
1455            variants_ref=merged_attrs['variants_ref'],
1456            variants_alt=merged_attrs['variants_alt'],
1457            variants_chrom=merged_attrs['variants_chrom'],
1458            variants_id=merged_attrs['variants_id'],
1459            variants_pos=merged_attrs['variants_pos'],
1460            variants_qual=merged_attrs['variants_qual'],
1461            variants_filter_pass=merged_attrs['variants_filter_pass'],
1462            ancestry_map=self.ancestry_map
1463        )
1464
1465    def remove_strand_ambiguous_variants(self, inplace: bool = False) -> Optional['SNPObject']:
1466        """
1467        A strand-ambiguous variant has reference (`variants_ref`) and alternate (`variants_alt`) alleles 
1468        in the pairs A/T, T/A, C/G, or G/C, where both alleles are complementary and thus indistinguishable 
1469        in terms of strand orientation.
1470
1471        Args:
1472            inplace (bool, default=False): 
1473                If True, modifies `self` in place. If False, returns a new `SNPObject` with the 
1474                strand-ambiguous variants removed. Default is False.
1475
1476        Returns:
1477            **Optional[SNPObject]:** A new `SNPObject` with non-ambiguous variants only if `inplace=False`. 
1478            If `inplace=True`, modifies `self` in place and returns None.
1479        """
1480        # Identify strand-ambiguous SNPs using vectorized comparisons
1481        is_AT = (self['variants_ref'] == 'A') & (self['variants_alt'] == 'T')
1482        is_TA = (self['variants_ref'] == 'T') & (self['variants_alt'] == 'A')
1483        is_CG = (self['variants_ref'] == 'C') & (self['variants_alt'] == 'G')
1484        is_GC = (self['variants_ref'] == 'G') & (self['variants_alt'] == 'C')
1485
1486        # Create a combined mask for all ambiguous variants
1487        ambiguous_mask = is_AT | is_TA | is_CG | is_GC
1488        non_ambiguous_idx = np.where(~ambiguous_mask)[0]
1489
1490        # Count each type of ambiguity using numpy's sum on boolean arrays
1491        A_T_count = np.sum(is_AT)
1492        T_A_count = np.sum(is_TA)
1493        C_G_count = np.sum(is_CG)
1494        G_C_count = np.sum(is_GC)
1495
1496        # Log the counts of each type of strand-ambiguous variants
1497        total_ambiguous = A_T_count + T_A_count + C_G_count + G_C_count
1498        log.info(f'{A_T_count} ambiguities of A-T type.')
1499        log.info(f'{T_A_count} ambiguities of T-A type.')
1500        log.info(f'{C_G_count} ambiguities of C-G type.')
1501        log.info(f'{G_C_count} ambiguities of G-C type.')
1502
1503        # Filter out ambiguous variants and keep non-ambiguous ones
1504        log.debug(f'Removing {total_ambiguous} strand-ambiguous variants...')
1505        return self.filter_variants(indexes=non_ambiguous_idx, include=True, inplace=inplace)
1506
1507    def correct_flipped_variants(
1508        self, 
1509        snpobj: 'SNPObject', 
1510        check_complement: bool = True, 
1511        index_by: str = 'pos', 
1512        common_variants_intersection: Optional[Tuple[np.ndarray, np.ndarray]] = None,
1513        log_stats: bool = True,
1514        inplace: bool = False
1515    ) -> Optional['SNPObject']:
1516        """
1517        Correct flipped variants between between `self` and a reference `snpobj`, where reference (`variants_ref`) 
1518        and alternate (`variants_alt`) alleles are swapped.
1519
1520        **Flip Detection Based on `check_complement`:**
1521
1522        - If `check_complement=False`, only direct allele swaps are considered:
1523            1. **Direct Swap:** `self.variants_ref == snpobj.variants_alt` and `self.variants_alt == snpobj.variants_ref`.
1524
1525        - If `check_complement=True`, both direct and complementary swaps are considered, with four possible cases:
1526            1. **Direct Swap:** `self.variants_ref == snpobj.variants_alt` and `self.variants_alt == snpobj.variants_ref`.
1527            2. **Complement Swap of Ref:** `complement(self.variants_ref) == snpobj.variants_alt` and `self.variants_alt == snpobj.variants_ref`.
1528            3. **Complement Swap of Alt:** `self.variants_ref == snpobj.variants_alt` and `complement(self.variants_alt) == snpobj.variants_ref`.
1529            4. **Complement Swap of both Ref and Alt:** `complement(self.variants_ref) == snpobj.variants_alt` and `complement(self.variants_alt) == snpobj.variants_ref`.
1530
1531        **Note:** Variants where `self.variants_ref == self.variants_alt` are ignored as they are ambiguous.
1532
1533        **Correction Process:** 
1534        - Swaps `variants_ref` and `variants_alt` alleles in `self` to align with `snpobj`.
1535        - Flips `calldata_gt` values (0 becomes 1, and 1 becomes 0) to match the updated allele configuration.
1536
1537        Args:
1538            snpobj (SNPObject): 
1539                The reference SNPObject to compare against.
1540            check_complement (bool, default=True): 
1541                If True, also checks for complementary base pairs (A/T, T/A, C/G, and G/C) when identifying swapped variants.
1542                Default is True.
1543            index_by (str, default='pos'): 
1544                Criteria for matching variants. Options:
1545                - `'pos'`: Matches by chromosome and position (`variants_chrom`, `variants_pos`), e.g., 'chr1-12345'.
1546                - `'id'`: Matches by variant ID alone (`variants_id`), e.g., 'rs123'.
1547                - `'pos+id'`: Matches by chromosome, position, and ID (`variants_chrom`, `variants_pos`, `variants_id`), e.g., 'chr1-12345-rs123'.
1548                Default is 'pos'.
1549            common_variants_intersection (tuple of arrays, optional): 
1550                Precomputed indices of common variants between `self` and `snpobj`. If None, intersection is 
1551                computed within the function.
1552            log_stats (bool, default=True): 
1553                If True, logs statistical information about matching and ambiguous alleles. Default is True.
1554            inplace (bool, default=False): 
1555                If True, modifies `self` in place. If False, returns a new `SNPObject` with corrected 
1556                flips. Default is False.
1557
1558        Returns:
1559            **Optional[SNPObject]**: 
1560                A new `SNPObject` with corrected flips if `inplace=False`. 
1561                If `inplace=True`, modifies `self` in place and returns None.
1562        """
1563        # Define complement mappings for nucleotides
1564        complement_map = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
1565
1566        # Helper function to get the complement of a base
1567        def get_complement(base: str) -> str:
1568            return complement_map.get(base, base)
1569
1570        # Get common variant indices if not provided
1571        if common_variants_intersection != None:
1572            query_idx, reference_idx = common_variants_intersection
1573        else:
1574            _, query_idx, reference_idx = self.get_common_variants_intersection(snpobj, index_by=index_by)
1575
1576        # Log statistics on matching alleles if enabled
1577        if log_stats:
1578            matching_ref = np.sum(self['variants_ref'][query_idx] == snpobj['variants_ref'][reference_idx])
1579            matching_alt = np.sum(self['variants_alt'][query_idx] == snpobj['variants_alt'][reference_idx])
1580            ambiguous = np.sum(self['variants_ref'][query_idx] == self['variants_alt'][query_idx])
1581            log.info(f"Matching reference alleles (ref=ref'): {matching_ref}, Matching alternate alleles (alt=alt'): {matching_alt}.")
1582            log.info(f"Number of ambiguous alleles (ref=alt): {ambiguous}.")
1583
1584        # Identify indices where `ref` and `alt` alleles are swapped
1585        if not check_complement:
1586            # Simple exact match for swapped alleles
1587            swapped_ref = (self['variants_ref'][query_idx] == snpobj['variants_alt'][reference_idx])
1588            swapped_alt = (self['variants_alt'][query_idx] == snpobj['variants_ref'][reference_idx])
1589        else:
1590            # Check for swapped or complementary-swapped alleles
1591            swapped_ref = (
1592                (self['variants_ref'][query_idx] == snpobj['variants_alt'][reference_idx]) |
1593                (np.vectorize(get_complement)(self['variants_ref'][query_idx]) == snpobj['variants_alt'][reference_idx])
1594            )
1595            swapped_alt = (
1596                (self['variants_alt'][query_idx] == snpobj['variants_ref'][reference_idx]) |
1597                (np.vectorize(get_complement)(self['variants_alt'][query_idx]) == snpobj['variants_ref'][reference_idx])
1598            )
1599
1600        # Filter out ambiguous variants where `ref` and `alt` alleles match (ref=alt)
1601        not_ambiguous = (self['variants_ref'][query_idx] != self['variants_alt'][query_idx])
1602
1603        # Indices in `self` of flipped variants
1604        flip_idx_query = query_idx[swapped_ref & swapped_alt & not_ambiguous]
1605
1606        # Correct the identified variant flips
1607        if len(flip_idx_query) > 0:
1608            log.info(f'Correcting {len(flip_idx_query)} variant flips...')
1609
1610            temp_alts = self['variants_alt'][flip_idx_query]
1611            temp_refs = self['variants_ref'][flip_idx_query]
1612
1613            # Correct the variant flips based on whether the operation is in-place or not
1614            if inplace:
1615                self['variants_alt'][flip_idx_query] = temp_refs
1616                self['variants_ref'][flip_idx_query] = temp_alts
1617                self['calldata_gt'][flip_idx_query] = 1 - self['calldata_gt'][flip_idx_query]
1618                return None
1619            else:
1620                snpobj = self.copy()
1621                snpobj['variants_alt'][flip_idx_query] = temp_refs
1622                snpobj['variants_ref'][flip_idx_query] = temp_alts
1623                snpobj['calldata_gt'][flip_idx_query] = 1 - snpobj['calldata_gt'][flip_idx_query]
1624                return snpobj
1625        else:
1626            log.info('No variant flips found to correct.')
1627            return self if not inplace else None
1628
1629    def remove_mismatching_variants(
1630        self, 
1631        snpobj: 'SNPObject', 
1632        index_by: str = 'pos', 
1633        common_variants_intersection: Optional[Tuple[np.ndarray, np.ndarray]] = None,
1634        inplace: bool = False
1635    ) -> Optional['SNPObject']:
1636        """
1637        Remove variants from `self`, where reference (`variants_ref`) and/or alternate (`variants_alt`) alleles 
1638        do not match with a reference `snpobj`.
1639
1640        Args:
1641            snpobj (SNPObject): 
1642                The reference SNPObject to compare against.
1643            index_by (str, default='pos'): 
1644                Criteria for matching variants. Options:
1645                - `'pos'`: Matches by chromosome and position (`variants_chrom`, `variants_pos`), e.g., 'chr1-12345'.
1646                - `'id'`: Matches by variant ID alone (`variants_id`), e.g., 'rs123'.
1647                - `'pos+id'`: Matches by chromosome, position, and ID (`variants_chrom`, `variants_pos`, `variants_id`), e.g., 'chr1-12345-rs123'.
1648                Default is 'pos'.
1649            common_variants_intersection (tuple of arrays, optional): 
1650                Precomputed indices of common variants between `self` and the reference `snpobj`.
1651                If None, the intersection is computed within the function.
1652            inplace (bool, default=False): 
1653                If True, modifies `self` in place. If False, returns a new `SNPObject` without 
1654                mismatching variants. Default is False.
1655
1656        Returns:
1657            **Optional[SNPObject]:** 
1658                A new `SNPObject` without mismatching variants if `inplace=False`. 
1659                If `inplace=True`, modifies `self` in place and returns None.
1660        """
1661        # Get common variant indices if not provided
1662        if common_variants_intersection is not None:
1663            query_idx, reference_idx = common_variants_intersection
1664        else:
1665            _, query_idx, reference_idx = self.get_common_variants_intersection(snpobj, index_by=index_by)
1666
1667        # Vectorized comparison of `ref` and `alt` alleles
1668        ref_mismatch = self['variants_ref'][query_idx] != snpobj['variants_ref'][reference_idx]
1669        alt_mismatch = self['variants_alt'][query_idx] != snpobj['variants_alt'][reference_idx]
1670        mismatch_mask = ref_mismatch | alt_mismatch
1671
1672        # Identify indices in `self` of mismatching variants
1673        mismatch_idx = query_idx[mismatch_mask]
1674
1675        # Compute total number of variant mismatches
1676        total_mismatches = np.sum(mismatch_mask)
1677
1678        # Filter out mismatching variants
1679        log.debug(f'Removing {total_mismatches} mismatching variants...')
1680        return self.filter_variants(indexes=mismatch_idx, include=True, inplace=inplace)
1681
1682    def shuffle_variants(self, inplace: bool = False) -> Optional['SNPObject']:
1683        """
1684        Randomly shuffle the positions of variants in the SNPObject, ensuring that all associated 
1685        data (e.g., `calldata_gt` and variant-specific attributes) remain aligned.
1686
1687        Args:
1688            inplace (bool, default=False): 
1689                If True, modifies `self` in place. If False, returns a new `SNPObject` with 
1690                shuffled variants. Default is False.
1691
1692        Returns:
1693            **Optional[SNPObject]:** 
1694                A new `SNPObject` without shuffled variant positions if `inplace=False`. 
1695                If `inplace=True`, modifies `self` in place and returns None.
1696        """
1697        # Generate a random permutation index for shuffling variant positions
1698        shuffle_index = np.random.permutation(self.n_snps)
1699
1700        # Apply shuffling to all relevant attributes using the class's dictionary-like interface
1701        if inplace:
1702            for key in self.keys():
1703                if self[key] is not None:
1704                    if key == 'calldata_gt':
1705                        # `calldata_gt`` has a different shape, so it's shuffled along axis 0
1706                        self[key] = self[key][shuffle_index, ...]
1707                    elif 'variant' in key:
1708                        # snpobj attributes are 1D arrays
1709                        self[key] = np.asarray(self[key])[shuffle_index]
1710            return None
1711        else:
1712            shuffled_snpobj = self.copy()
1713            for key in shuffled_snpobj.keys():
1714                if shuffled_snpobj[key] is not None:
1715                    if key == 'calldata_gt':
1716                        shuffled_snpobj[key] = shuffled_snpobj[key][shuffle_index, ...]
1717                    elif 'variant' in key:
1718                        shuffled_snpobj[key] = np.asarray(shuffled_snpobj[key])[shuffle_index]
1719            return shuffled_snpobj
1720
1721    def set_empty_to_missing(self, inplace: bool = False) -> Optional['SNPObject']:
1722        """
1723        Replace empty strings `''` with missing values `'.'` in attributes of `self`.
1724
1725        Args:
1726            inplace (bool, default=False): 
1727                If True, modifies `self` in place. If False, returns a new `SNPObject` with empty 
1728                strings `''` replaced by missing values `'.'`. Default is False.
1729
1730        Returns:
1731            **Optional[SNPObject]:** 
1732                A new `SNPObject` with empty strings replaced if `inplace=False`. 
1733                If `inplace=True`, modifies `self` in place and returns None.
1734        """
1735        if inplace:
1736            if self.variants_alt is not None:
1737                self.variants_alt[self.variants_alt == ''] = '.'
1738            if self.variants_ref is not None:
1739                self.variants_ref[self.variants_ref == ''] = '.'
1740            if self.variants_qual is not None:
1741                self.variants_qual = self.variants_qual.astype(str)
1742                self.variants_qual[(self.variants_qual == '') | (self.variants_qual == 'nan')] = '.'
1743            if self.variants_chrom is not None:
1744                self.variants_chrom = self.variants_chrom.astype(str)
1745                self.variants_chrom[self.variants_chrom == ''] = '.'
1746            if self.variants_filter_pass is not None:
1747                self.variants_filter_pass[self.variants_filter_pass == ''] = '.'
1748            if self.variants_id is not None:
1749                self.variants_id[self.variants_id == ''] = '.'
1750            return self
1751        else:
1752            snpobj = self.copy()
1753            if snpobj.variants_alt is not None:
1754                snpobj.variants_alt[snpobj.variants_alt == ''] = '.'
1755            if snpobj.variants_ref is not None:
1756                snpobj.variants_ref[snpobj.variants_ref == ''] = '.'
1757            if snpobj.variants_qual is not None:
1758                snpobj.variants_qual = snpobj.variants_qual.astype(str)
1759                snpobj.variants_qual[(snpobj.variants_qual == '') | (snpobj.variants_qual == 'nan')] = '.'
1760            if snpobj.variants_chrom is not None:
1761                snpobj.variants_chrom[snpobj.variants_chrom == ''] = '.'
1762            if snpobj.variants_filter_pass is not None:
1763                snpobj.variants_filter_pass[snpobj.variants_filter_pass == ''] = '.'
1764            if snpobj.variants_id is not None:
1765                snpobj.variants_id[snpobj.variants_id == ''] = '.'
1766            return snpobj
1767
1768    def convert_to_window_level(
1769        self,
1770        window_size: Optional[int] = None,
1771        physical_pos: Optional[np.ndarray] = None,
1772        chromosomes: Optional[np.ndarray] = None,
1773        window_sizes: Optional[np.ndarray] = None,
1774        laiobj: Optional['LocalAncestryObject'] = None
1775    ) -> 'LocalAncestryObject':
1776        """
1777        Aggregate the `calldata_lai` attribute into genomic windows within a 
1778        `snputils.ancestry.genobj.LocalAncestryObject`.
1779
1780        **Options for defining windows (in order of precedence):**
1781
1782        1. **Fixed window size**:
1783        - Use `window_size` to specify how many SNPs go into each window. The last window on each 
1784        chromosome may be larger if SNPs are not evenly divisible by the size.
1785
1786        2. **Custom start and end positions**:
1787        - Provide `physical_pos` (2D array of shape (n_windows, 2)) as the [start, end] base-pair 
1788         coordinates for each window. 
1789        - If `chromosomes` is not provided and `self` has exactly one chromosome, all windows are 
1790        assumed to belong to that chromosome. 
1791        - If multiple chromosomes exist but `chromosomes` is missing, an error will be raised.
1792        - Optionally, provide `window_sizes` to store the SNP count per-window.
1793
1794        3. **Matching existing windows**:
1795        - Reuse window definitions (`physical_pos`, `chromosomes`, `window_sizes`) from an existing `laiobj`.
1796
1797        Args:
1798            window_size (int, optional): 
1799                Number of SNPs in each window if defining fixed-size windows. If the total number of 
1800                SNPs in a chromosome is not evenly divisible by the window size, the last window on that 
1801                chromosome will include all remaining SNPs and therefore be larger than the specified size.
1802            physical_pos (array of shape (n_windows, 2), optional): 
1803                A 2D array containing the start and end physical positions for each window.
1804            chromosomes (array of shape (n_windows,), optional): 
1805                An array with chromosome numbers corresponding to each genomic window.
1806            window_sizes (array of shape (n_windows,), optional): 
1807                An array specifying the number of SNPs in each genomic window.
1808            laiobj (LocalAncestryObject, optional): 
1809                A reference `LocalAncestryObject` from which to copy existing window definitions.
1810
1811        Returns:
1812            **LocalAncestryObject:** 
1813                A LocalAncestryObject containing window-level ancestry data.
1814        """
1815        from snputils.ancestry.genobj.local import LocalAncestryObject
1816
1817        if window_size is None and physical_pos is None and laiobj is None:
1818            raise ValueError("One of `window_size`, `physical_pos`, or `laiobj` must be provided.")
1819        
1820        # Fixed window size
1821        if window_size is not None:
1822            physical_pos = []   # Boundaries [start, end] of each window
1823            chromosomes = []    # Chromosome for each window
1824            window_sizes = []   # Number of SNPs for each window
1825            for chrom in self.unique_chrom:
1826                # Extract indices corresponding to this chromosome
1827                mask_chrom = (self.variants_chrom == chrom)
1828                # Subset to this chromosome
1829                pos_chrom = self.variants_pos[mask_chrom]
1830                # Number of SNPs for this chromosome
1831                n_snps_chrom = pos_chrom.size
1832                
1833                # Initialize the start of the first window with the position of the first SNP
1834                current_start = self.variants_pos[0]
1835
1836                # Number of full windows with exactly `window_size` SNPs
1837                n_full_windows = n_snps_chrom // window_size
1838
1839                # Build all but the last window
1840                for i in range(n_full_windows-1):
1841                    current_end = self.variants_pos[(i+1) * window_size - 1]
1842                    physical_pos.append([current_start, current_end])
1843                    chromosomes.append(chrom)
1844                    window_sizes.append(window_size)
1845                    current_start = self.variants_pos[(i+1) * window_size]
1846                
1847                # Build the last window
1848                current_end = self.variants_pos[-1]
1849                physical_pos.append([current_start, current_end])
1850                chromosomes.append(chrom)
1851                window_sizes.append(n_snps_chrom - ((n_full_windows - 1) * window_size))
1852                
1853            physical_pos = np.array(physical_pos)
1854            chromosomes = np.array(chromosomes)
1855            window_sizes = np.array(window_sizes)
1856        
1857        # Custom start and end positions
1858        elif physical_pos is not None:
1859            # Check if there is exactly one chromosome
1860            if chromosomes is None:
1861                unique_chrom = self.unique_chrom
1862                if len(unique_chrom) == 1:
1863                    # We assume all windows belong to this single chromosome
1864                    single_chrom = unique_chrom[0]
1865                    chromosomes = np.array([single_chrom] * physical_pos.shape[0])
1866                else:
1867                    raise ValueError("Multiple chromosomes detected, but `chromosomes` was not provided.")
1868
1869        # Match existing windows to a reference laiobj
1870        elif laiobj is not None:
1871            physical_pos = laiobj.physical_pos
1872            chromosomes = laiobj.chromosomes
1873            window_sizes = laiobj.window_sizes
1874
1875        # Allocate an output LAI array
1876        n_windows = physical_pos.shape[0]
1877        n_samples = self.n_samples
1878        if self.calldata_lai.ndim == 3:
1879            lai = np.zeros((n_windows, n_samples, 2))
1880        else:
1881            lai = np.zeros((n_windows, n_samples*2))
1882
1883        # For each window, find the relevant SNPs and compute the mode of the ancestries
1884        for i, ((start, end), chrom) in enumerate(zip(physical_pos, chromosomes)):
1885            snps_mask = (
1886                (self.variants_chrom == chrom) &
1887                (self.variants_pos >= start) &
1888                (self.variants_pos <= end)
1889            )
1890            if np.any(snps_mask):
1891                lai_mask = self.calldata_lai[snps_mask, ...]
1892                mode_ancestries = mode(lai_mask, axis=0, nan_policy='omit').mode
1893                lai[i] = mode_ancestries
1894            else:
1895                lai[i] = np.nan
1896
1897        # Generate haplotype labels, e.g. "Sample1.0", "Sample1.1"
1898        haplotypes = [f"{sample}.{i}" for sample in self.samples for i in range(2)]
1899
1900        # If original data was (n_snps, n_samples, 2), flatten to (n_windows, n_samples*2)
1901        if self.calldata_lai.ndim == 3:
1902            lai = lai.reshape(n_windows, -1)
1903
1904        # Aggregate into a LocalAncestryObject
1905        return LocalAncestryObject(
1906            haplotypes=haplotypes,
1907            lai=lai,
1908            samples=self.samples,
1909            ancestry_map=self.ancestry_map,
1910            window_sizes=window_sizes,
1911            physical_pos=physical_pos,
1912            chromosomes=chromosomes
1913        )
1914
1915    def save(self, file: Union[str, Path]) -> None:
1916        """
1917        Save the data stored in `self` to a specified file.
1918
1919        The format of the saved file is determined by the file extension provided in the `file` 
1920        argument. 
1921        
1922        **Supported formats:**
1923        
1924        - `.bed`: Binary PED (Plink) format.
1925        - `.pgen`: Plink2 binary genotype format.
1926        - `.vcf`: Variant Call Format.
1927        - `.pkl`: Pickle format for saving `self` in serialized form.
1928
1929        Args:
1930            file (str or pathlib.Path): 
1931                Path to the file where the data will be saved. The extension of the file determines the save format. 
1932                Supported extensions: `.bed`, `.pgen`, `.vcf`, `.pkl`.
1933        """
1934        ext = Path(file).suffix.lower()
1935        if ext == '.bed':
1936            self.save_bed(file)
1937        elif ext == '.pgen':
1938            self.save_pgen(file)
1939        elif ext == '.vcf':
1940            self.save_vcf(file)
1941        elif ext == '.pkl':
1942            self.save_pickle(file)
1943        else:
1944            raise ValueError(f"Unsupported file extension: {ext}")
1945
1946    def save_bed(self, file: Union[str, Path]) -> None:
1947        """
1948        Save the data stored in `self` to a `.bed` file.
1949
1950        Args:
1951            file (str or pathlib.Path): 
1952                Path to the file where the data will be saved. It should end with `.bed`. 
1953                If the provided path does not have this extension, it will be appended.
1954        """
1955        from snputils.snp.io.write.bed import BEDWriter
1956        writer = BEDWriter(snpobj=self, filename=file)
1957        writer.write()
1958
1959    def save_pgen(self, file: Union[str, Path]) -> None:
1960        """
1961        Save the data stored in `self` to a `.pgen` file.
1962
1963        Args:
1964            file (str or pathlib.Path): 
1965                Path to the file where the data will be saved. It should end with `.pgen`. 
1966                If the provided path does not have this extension, it will be appended.
1967        """
1968        from snputils.snp.io.write.pgen import PGENWriter
1969        writer = PGENWriter(snpobj=self, filename=file)
1970        writer.write()
1971
1972    def save_vcf(self, file: Union[str, Path]) -> None:
1973        """
1974        Save the data stored in `self` to a `.vcf` file.
1975
1976        Args:
1977            file (str or pathlib.Path): 
1978                Path to the file where the data will be saved. It should end with `.vcf`. 
1979                If the provided path does not have this extension, it will be appended.
1980        """
1981        from snputils.snp.io.write.vcf import VCFWriter
1982        writer = VCFWriter(snpobj=self, filename=file)
1983        writer.write()
1984
1985    def save_pickle(self, file: Union[str, Path]) -> None:
1986        """
1987        Save `self` in serialized form to a `.pkl` file.
1988
1989        Args:
1990            file (str or pathlib.Path): 
1991                Path to the file where the data will be saved. It should end with `.pkl`. 
1992                If the provided path does not have this extension, it will be appended.
1993        """
1994        import pickle
1995        with open(file, 'wb') as file:
1996            pickle.dump(self, file)
1997
1998    @staticmethod
1999    def _match_to_replace(val: Union[str, int, float], dictionary: Dict[Any, Any], regex: bool = True) -> Union[str, int, float]:
2000        """
2001        Find a matching key in the provided dictionary for the given value `val`
2002        and replace it with the corresponding value.
2003
2004        Args:
2005            val (str, int, or float): 
2006                The value to be matched and potentially replaced.
2007            dictionary (Dict): 
2008                A dictionary containing keys and values for matching and replacement.
2009                The keys should match the data type of `val`.
2010            regex (bool): 
2011                If True, interprets keys in `dictionary` as regular expressions.
2012                Default is True.
2013
2014        Returns:
2015            str, int, or float: 
2016                The replacement value from `dictionary` if a match is found; otherwise, the original `val`.
2017        """
2018        if regex:
2019            # Use regular expression matching to find replacements
2020            for key, value in dictionary.items():
2021                if isinstance(key, str):
2022                    match = re.match(key, val)
2023                    if match:
2024                        # Replace using the first matching regex pattern
2025                        return re.sub(key, value, val)
2026            # Return the original value if no regex match is found
2027            return val
2028        else:
2029            # Return the value for `val` if present in `dictionary`; otherwise, return `val`
2030            return dictionary.get(val, val)
2031
2032    @staticmethod
2033    def _get_chromosome_number(chrom_string: str) -> Union[int, str]:
2034        """
2035        Extracts the chromosome number from the given chromosome string.
2036
2037        Args:
2038            chrom_string (str): 
2039                The chromosome identifier.
2040
2041        Returns:
2042            int or str: 
2043                The numeric representation of the chromosome if detected. 
2044                Returns 10001 for 'X' or 'chrX', 10002 for 'Y' or 'chrY', 
2045                and the original `chrom_string` if unrecognized.
2046        """
2047        if chrom_string.isdigit():
2048            return int(chrom_string)
2049        else:
2050            chrom_num = re.search(r'\d+', chrom_string)
2051            if chrom_num:
2052                return int(chrom_num.group())
2053            elif chrom_string.lower() in ['x', 'chrx']:
2054                return 10001
2055            elif chrom_string.lower() in ['y', 'chry']:
2056                return 10002
2057            else:
2058                log.warning(f"Chromosome nomenclature not standard. Chromosome: {chrom_string}")
2059                return chrom_string
2060
2061    def _sanity_check(self) -> None:
2062        """
2063        Perform sanity checks to ensure LAI and ancestry map consistency.
2064
2065        This method checks that all unique ancestries in the LAI data are represented 
2066        in the ancestry map if it is provided.
2067        """
2068        if self.__calldata_lai is not None and self.__ancestry_map is not None:
2069            unique_ancestries = np.unique(self.__calldata_lai)
2070            missing_ancestries = [anc for anc in unique_ancestries if str(anc) not in self.__ancestry_map]
2071            if missing_ancestries:
2072                warnings.warn(f"Missing ancestries in ancestry_map: {missing_ancestries}")