scmkl.estimate_sigma

  1import numpy as np
  2import scipy
  3from sklearn.decomposition import TruncatedSVD, PCA
  4
  5from scmkl.calculate_z import _process_data
  6from scmkl.tfidf_normalize import _tfidf
  7
  8def estimate_sigma(adata, n_features = 5000):
  9    '''
 10    Calculate kernel widths to inform distribution for projection of 
 11    Fourier Features. Calculates one sigma per group of features.
 12
 13    Parameters
 14    ----------
 15    **adata** : *AnnData* 
 16        > Created by `create_adata`.
 17    
 18    **n_features** : *int*  
 19        > Number of random features to include when estimating sigma. 
 20        Will be scaled for the whole pathway set according to a 
 21        heuristic. Used for scalability.
 22    
 23    Returns
 24    -------
 25    **adata** : *AnnData*
 26        > Key added `adata.uns['sigma']`.
 27
 28    Examples
 29    --------
 30    >>> adata = scmkl.estimate_sigma(adata)
 31    >>> adata.uns['sigma']
 32    array([10.4640895 , 10.82011454,  6.16769438,  9.86156855, ...])
 33    '''
 34    sigma_list = []
 35
 36    # Loop over every group in group_dict
 37    for m, group_features in enumerate(adata.uns['group_dict'].values()):
 38
 39        # Select only features in that group and downsample for scalability
 40        num_group_features = len(group_features)
 41        group_array = np.array(list(group_features))
 42        n_feats = min([n_features, num_group_features])
 43        group_features = adata.uns['seed_obj'].choice(group_array, n_feats, 
 44                                                      replace = False) 
 45
 46        # Use on the train data to estimate sigma
 47        X_train = adata[adata.uns['train_indices'], group_features].X
 48
 49        
 50        # Sample cells for scalability
 51        sample_idx = np.arange(X_train.shape[0])
 52        n_samples = np.min((2000, X_train.shape[0]))
 53        distance_indices = adata.uns['seed_obj'].choice(sample_idx, n_samples, 
 54                                                        replace = False)
 55
 56        X_train = _process_data(X_train = X_train, 
 57                                 scale_data = adata.uns['scale_data'], 
 58                                 return_dense = True)
 59        
 60
 61        if adata.uns['tfidf']:
 62            X_train = _tfidf(X_train, mode = 'normalize')
 63
 64        if adata.uns['reduction'].lower() == 'svd':
 65
 66            SVD_func = TruncatedSVD(n_components = np.min([50, X_train.shape[1]]), random_state = 1)
 67            X_train = SVD_func.fit_transform(scipy.sparse.csr_array(X_train[distance_indices,:]))[:,1:]
 68
 69        elif adata.uns['reduction'].lower() == 'pca':
 70            PCA_func = PCA(n_components = np.min([50, X_train.shape[1]]), random_state = 1)
 71            X_train = PCA_func.fit_transform(np.asarray(X_train[distance_indices,:]))
 72
 73        elif adata.uns['reduction'].lower() == 'linear':
 74
 75            X_train = X_train[distance_indices] @ adata.uns['seed_obj'].choice([0,1], [0.02, 0.98]).reshape((len(distance_indices), 50))
 76
 77        else:
 78            X_train = X_train[distance_indices, :]
 79
 80        if scipy.sparse.issparse(X_train):
 81            X_train = X_train.toarray()
 82
 83        # Calculate Distance Matrix with specified metric
 84        sigma = scipy.spatial.distance.cdist(X_train, 
 85                                             X_train, 
 86                                             adata.uns['distance_metric'])
 87        sigma = np.mean(sigma)
 88
 89        # sigma = 0 is numerically unusable in later steps
 90        # Using such a small sigma will result in wide distribution, and 
 91        # typically a non-predictive Z
 92        if sigma == 0:
 93            sigma += 1e-5
 94
 95        if n_features < num_group_features:
 96            # Heuristic we calculated to account for fewer features used in 
 97            # distance calculation
 98            sigma = sigma * num_group_features / n_features 
 99
100        sigma_list.append(sigma)
101    
102    adata.uns['sigma'] = np.array(sigma_list)
103        
104    return adata