Repository: rmarkello/abagen
Pull request #119 (2019-11-18 13:56:46)

[REF,ENH] Modifies and adds normalization procedures
Adds new normalization functions and modifies when/where normalization occurs in primary workflow.

To Do:
- [x] Add sample normalization option to `get_expression_data()` workflow
- [x] Add tests for new normalization procedures

coveralls: ## Pull Request Test Coverage Report for [Build 334](https://coveralls.io/builds/27058868)

* **0** of **0**   changed or added relevant lines in **0** files are covered.
* **6** unchanged lines in **2** files lost coverage.
* Overall coverage decreased (**-0.03%**) to **91.353%**

---


|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [/home/travis/virtualenv/python3.6.7/lib/python3.6/site-packages/abagen/allen.py](https://coveralls.io/builds/27058868/source?filename=%2Fhome%2Ftravis%2Fvirtualenv%2Fpython3.6.7%2Flib%2Fpython3.6%2Fsite-packages%2Fabagen%2Fallen.py#L362) | 2 | 97.98% |
| [/home/travis/virtualenv/python3.6.7/lib/python3.6/site-packages/abagen/correct.py](https://coveralls.io/builds/27058868/source?filename=%2Fhome%2Ftravis%2Fvirtualenv%2Fpython3.6.7%2Flib%2Fpython3.6%2Fsite-packages%2Fabagen%2Fcorrect.py#L239) | 4 | 96.64% |
<!-- | **Total:** | **6** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/27058868/badge)](https://coveralls.io/builds/27058868) |
| :-- | --: |
| Change from base [Build 330](https://coveralls.io/builds/26212965): |  -0.03% |
| Covered Lines: | 17558 |
| Relevant Lines: | 19220 |

---
##### 💛  - [Coveralls](https://coveralls.io)

codecov[bot]: # [Codecov](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=h1) Report
> Merging [#119](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=desc) into [master](https://codecov.io/gh/rmarkello/abagen/commit/9137748d7e2896d445ab956487f0a4765a7d7e51?src=pr&el=desc) will **decrease** coverage by `0.03%`.
> The diff coverage is `92.68%`.

[![Impacted file tree graph](https://codecov.io/gh/rmarkello/abagen/pull/119/graphs/tree.svg?width=650&token=jX4aAk8SDv&height=150&src=pr)](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #119      +/-   ##
==========================================
- Coverage   91.45%   91.41%   -0.04%     
==========================================
  Files          33       33              
  Lines        1895     1922      +27     
==========================================
+ Hits         1733     1757      +24     
- Misses        162      165       +3
```


| [Impacted Files](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [abagen/allen.py](https://codecov.io/gh/rmarkello/abagen/pull/119/diff?src=pr&el=tree#diff-YWJhZ2VuL2FsbGVuLnB5) | `97.97% <100%> (+0.17%)` | :arrow_up: |
| [abagen/tests/test\_correct.py](https://codecov.io/gh/rmarkello/abagen/pull/119/diff?src=pr&el=tree#diff-YWJhZ2VuL3Rlc3RzL3Rlc3RfY29ycmVjdC5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [abagen/correct.py](https://codecov.io/gh/rmarkello/abagen/pull/119/diff?src=pr&el=tree#diff-YWJhZ2VuL2NvcnJlY3QucHk=) | `96.63% <88.88%> (-2.38%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=footer). Last update [9137748...efceb34](https://codecov.io/gh/rmarkello/abagen/pull/119?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/abagen/allen.py b/abagen/allen.py
index af6aa95..4d58c08 100644
--- a/abagen/allen.py
+++ b/abagen/allen.py
@@ -56,17 +56,20 @@ def groupby_label(microarray, sample_labels, labels=None, metric='mean'):
                                .groupby('label')
                                .aggregate(metric)
                                .append(labels)
-                               .drop([0])
                                .sort_index()
                                .rename_axis('label'))
 
+    # remove "zero" label (if it exists)
+    if 0 in gene_by_label.index:
+        gene_by_label = gene_by_label.drop([0], axis=0)
+
     return gene_by_label
 
 
 def get_expression_data(atlas, atlas_info=None, *, exact=True,
                         tolerance=2, metric='mean', ibf_threshold=0.5,
-                        probe_selection='diff_stability',
-                        lr_mirror=False, donor_norm='srs',
+                        probe_selection='diff_stability', lr_mirror=False,
+                        sample_norm='srs', donor_norm='srs',
                         corrected_mni=True, reannotated=True,
                         return_counts=False, return_donors=False,
                         donors='all', data_dir=None, verbose=1, n_proc=1):
@@ -110,9 +113,9 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
     parameter; see the parameter description for more information.
 
     Once all samples have been matched to parcels for all supplied donors, the
-    microarray expression data are optionally normalized within-donor via the
-    provided `donor_norm` function before being combined across donors via the
-    supplied `metric`.
+    microarray expression data are optionally normalized via the provided
+    `sample_norm` and `donor_norm` functions before being combined within
+    parcels and across donors via the supplied `metric.
 
     Parameters
     ----------
@@ -160,12 +163,18 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         increase spatial coverage. This will duplicate samples across both
         hemispheres (i.e., L->R and R->L), approximately doubling the number of
         available samples. Default: False
-    donor_norm : {'srs', 'zscore', 'batch', None}, optional
+    sample_norm : {'rs', 'srs', 'minmax', 'center', 'zscore', None}, optional
+        Method by which to normalize microarray expression values for each
+        sample. Expression values are normalized separately for each sample for
+        each donor across all genes; see Notes for more information on
+        different methods. If None is specified then no normalization is
+        performed. Default: 'srs'
+    donor_norm : {'rs', 'srs', 'minmax', 'center', 'zscore', None}, optional
         Method by which to normalize microarray expression values for each
         donor. Expression values are normalized separately for each gene for
         each donor across all regions in `atlas`; see Notes for more
-        information on different methods. If not specified no normalization
-        is performed. Default: 'srs'
+        information on different methods. If None is specified then no
+        normalization is performed. Default: 'srs'
     corrected_mni : bool, optional
         Whether to use the "corrected" MNI coordinates shipped with the
         `alleninf` package instead of the coordinates provided with the AHBA
@@ -257,25 +266,26 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
     The following methods can be used for normalizing microarray expression
     values for each donor prior to aggregating:
 
-        1. ``donor_norm='srs'``
+        1. ``norm=='rs'``
+
+        Uses a robust sigmoid function as in [A3]_ to normalize values
+
+        2. ``norm='srs'``
 
-        Uses a scaled robust sigmoid function as in [A3]_ to normalize
-        expression values for each gene across regions to within the unit
-        normal (i.e., in the range 0-1).
+        Same as 'rs' but scales output to the unit normal (i.e., range 0-1)
 
-        2. ``donor_norm='zscore'``
+        3. ``norm='minmax'``
 
-        Applies a basic z-score (subtract mean, divide by standard deviation)
-        to expression values for each gene across regions. Uses degrees of
-        freedom equal to one for standard deviation calculation.
+        Scales data in each column to the unit normal (i.e., range 0-1)
 
-        3. ``donor_norm='batch'``
+        4. ``norm='center'``
 
-        Uses a linear model to remove donor effects from expression values.
-        Differs from other methods in that all donors are simultaneously fit
-        to the same model and expression values are residualized based on
-        estimated betas. Linear model includes the intercept but the
-        residualization does not remove it.
+        Removes the mean of expression values
+
+        5. ``norm='zscore'``
+
+        Applies a basic z-score (subtract mean, divide by standard deviation);
+        uses degrees of freedom equal to one for standard deviation
 
     References
     ----------
@@ -388,13 +398,33 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         annotation[subj] = samples.drop_mismatch_samples(annotation[subj],
                                                          ontology[subj])
 
-        # subset representative probes + samples from microarray data
-        microarray[subj] = microarray[subj].loc[annotation[subj].index]
-
-        # assign samples to regions and aggregate samples w/i the same region
+        # assign samples to regions
         labels = samples.label_samples(annotation[subj], atlas,
                                        atlas_info=atlas_info,
                                        tolerance=tolerance)
+
+        # subset representative probes + samples from microarray data
+        microarray[subj] = microarray[subj].loc[annotation[subj].index]
+
+        if exact:  # remove all samples not assigned a label before norming
+            nz = np.asarray(labels != 0).squeeze()
+            microarray[subj] = microarray[subj].loc[nz]
+            labels = labels[nz]
+
+        if sample_norm is not None:
+            if subj == 0:  # only provide norm log once
+                lgr.info('Normalizing expression across genes with function: '
+                         '{}'.format(sample_norm))
+            microarray[subj] = correct.normalize_expression(microarray[subj].T,
+                                                            norm=sample_norm).T
+        if donor_norm is not None:
+            if subj == 0:  # only provide norm log once
+                lgr.info('Normalizing expression across sample with function: '
+                         '{}'.format(donor_norm))
+            microarray[subj] = correct.normalize_expression(microarray[subj],
+                                                            norm=donor_norm)
+
+        # aggregate normalized samples within regions
         expression += [groupby_label(microarray[subj], labels,
                                      all_labels, metric=metric)]
         lgr.info('{:>3} / {} samples matched to regions for donor #{}'
@@ -412,7 +442,7 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         if not exact:
             cols = ['mni_x', 'mni_y', 'mni_z']
             coords = utils.xyz_to_ijk(annotation[subj][cols], atlas.affine)
-            empty = ~np.in1d(all_labels, labs)
+            empty = np.logical_not(np.in1d(all_labels, labs))
             idx, dist = utils.closest_centroid(coords, centroids[empty],
                                                return_dist=True)
             idx = microarray[subj].loc[annotation[subj].iloc[idx].index]
@@ -436,17 +466,14 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
             expression[ind].loc[roi] = missing[ind][0].loc[roi]
             counts.loc[roi, ind] += 1
 
-    # normalize data with SRS
-    if donor_norm is not None:
-        lgr.info('Normalizing donor expression data with function: "{}"'
-                 .format(donor_norm))
-        expression = correct.normalize_expression(expression, norm=donor_norm)
-
     # aggregate across donors if individual donor dataframes not requested
     if not return_donors:
         lgr.info('Aggregating donor expression data with function: "{}"'
                  .format(metric))
         expression = pd.concat(expression).groupby('label').aggregate(metric)
+        # some genes may have been poorly normalized; remove these
+        drop = expression.dropna(axis=0, how='all').isna().any(axis=0)
+        expression = expression.drop(drop[drop].index, axis=1)
 
     # drop the "zero" label from the counts dataframe (this is background)
     if return_counts:
diff --git a/abagen/correct.py b/abagen/correct.py
index 0fad225..9f7d7ba 100644
--- a/abagen/correct.py
+++ b/abagen/correct.py
@@ -4,6 +4,7 @@ Functions for processing and correcting gene expression data
 """
 
 import itertools
+import warnings
 
 import numpy as np
 import pandas as pd
@@ -13,6 +14,24 @@ from scipy.spatial.distance import cdist
 from . import utils
 
 
+def _unpack_tuple(var):
+    """
+    Returns first entry of var if there is only one entry
+
+    Parameters
+    ----------
+    var : tuple
+
+    Returns
+    -------
+    var
+    """
+
+    if len(var) == 1:
+        return var[0]
+    return var
+
+
 def _batch_correct(data):
     """
     Performs batch correction on `data`
@@ -52,14 +71,44 @@ def _batch_correct(data):
     return residualized
 
 
-def _srs(data, axis=0):
+def _rescale(data, low=0, high=1, axis=0):
     """
-    Normalizes `data` with a scaled robust sigmoid function
+    Rescales `data` to range [`low`, `high`]
+
+    Parameters
+    ----------
+    data : array_like
+        Input data array to be rescaled
+    low : float, optional
+        Lower bound for rescaling. Default: 0
+    high : float, optional
+        Upper bound for rescaling. Default: 1
+    axis : int, optional
+        Axis of `data` to scale along. Default: 0
+
+    Returns
+    -------
+    rescaled : np.ndarray
+        Rescaled data
+    """
+
+    data = np.asarray(data)
+
+    dmin = data.min(axis=axis, keepdims=True)
+    dmax = data.max(axis=axis, keepdims=True)
+    rescaled = ((data - dmin) / (dmax - dmin)) * (high - low) + low
+
+    return rescaled
+
+
+def _rs(data, axis=0):
+    """
+    Normalizes `data` with a robust sigmoid function
 
     Parameters
     ----------
     data : array_like
-        Input data to be passed through SRS function
+        Input data array to be transformed
     axis : int, optional
         Axis of `data` to be normalized
 
@@ -74,12 +123,35 @@ def _srs(data, axis=0):
     # calculate sigmoid normalization
     med = np.median(data, axis=axis, keepdims=True)
     iqr = sstats.iqr(data, axis=axis, scale='normal', keepdims=True)
-    srs = 1 / (1 + np.exp(-(data - med) / iqr))
+    with np.errstate(divide='ignore', invalid='ignore'):
+        rs = 1 / (1 + np.exp(-(data - med) / iqr))
+
+    return rs
+
+
+def _srs(data, low=0, high=1, axis=0):
+    """
+    Normalizes `data` with a scaled robust sigmoid function
+
+    Parameters
+    ----------
+    data : array_like
+        Input data array to be transformed
+    low : float, optional
+        Lower bound for rescaling. Default: 0
+    high : float, optional
+        Upper bound for rescaling. Default: 1
+    axis : int, optional
+        Axis of `data` to be normalized
 
-    # rescale normalized values to a unit interval
-    srs_min = srs.min(axis=axis, keepdims=True)
-    srs_max = srs.max(axis=axis, keepdims=True)
-    normed = (srs - srs_min) / (srs_max - srs_min)
+    Returns
+    -------
+    normed : array_like
+        Normalized input `data`
+    """
+
+    srs = _rs(np.asarray(data), axis=axis)  # robust sigmoid transform
+    normed = _rescale(srs, low=low, high=high, axis=axis)  # scales data
 
     return normed
 
@@ -93,8 +165,8 @@ def normalize_expression(expression, norm='srs'):
     expression : list of (S, G) pandas.DataFrame
         Microarray expression data, where `S` is samples (or regions) and `G`
         is genes
-    norm : {'srs', 'zscore', 'batch'}, optional
-        Function by which to normalize expression data; see Notes for more
+    norm : {'rs', 'srs', 'center', 'zscore', 'minmax', 'batch'}, optional
+        Function with which to normalize expression data. See Notes for more
         information. Default: 'srs'
 
     Returns
@@ -105,36 +177,48 @@ def normalize_expression(expression, norm='srs'):
     Notes
     -----
     The following methods can be used for normalizing gene expression values
-    for each donor:
+    for each donor (adapted from [PC2]_):
+
+    1. ``norm=='rs'``
 
-    1. ``norm='srs'``
+    Uses a robust sigmoid function ([PC1]_) to normalize data in each column
 
-    Uses a scaled robust sigmoid function as in [PC1]_ to normalize expression
-    values for each gene across regions to within the unit normal (i.e., in the
-    range 0-1).
+    2. ``norm='srs'``
 
-    2. ``norm='zscore'``
+    Same as 'rs' but scales output to the unit normal (i.e., range 0-1)
+
+    3. ``norm='center'``
+
+    Removes the mean of data in each column
+
+    4. ``norm='zscore'``
 
     Applies a basic z-score (subtract mean, divide by standard deviation) to
-    expression values for each gene across regions. Uses degrees of freedom
-    equal to one for standard deviation calculation.
+    each column; uses degrees of freedom equal to one for standard deviation
 
-    3. ``norm='batch'``
+    5. ``norm='minmax'``
 
-    Uses a linear model to remove donor effects from expression values. Differs
-    from other methods in that all donors are simultaneously fit to the same
-    model and expression values are residualized based on estimated betas.
-    Linear model includes the intercept but the residualization does not remove
-    it.
+    Scales data in each column to the unit normal (i.e., range 0-1)
+
+    6. ``norm='batch'``
+
+    Uses a linear model to remove donor effects from data. Differs from other
+    methods in that all donors are simultaneously fit to the same model and
+    data are residualized based on estimated betas. Linear model includes the
+    intercept but it is not removed during residualization
 
     References
     ----------
     .. [PC1] Fulcher, B. D., & Fornito, A. (2016). A transcriptional signature
        of hub connectivity in the mouse connectome. Proceedings of the National
        Academy of Sciences, 113(5), 1435-1440.
+    .. [PC2] Fulcher, B. D., Little, M. A., & Jones, N. S. (2013). Highly
+       comparative time-series analysis: the empirical structure of time series
+       and their methods. Journal of the Royal Society Interface, 10(83),
+       20130048
     """
 
-    norms = ['srs', 'zscore', 'batch']
+    norms = ['rs', 'srs', 'center', 'zscore', 'batch', 'minmax']
     if norm not in norms:
         raise ValueError('Provided value for `norm` not recognized. Must be '
                          'one of {}. Received: {}'.format(norms, norm))
@@ -149,13 +233,19 @@ def normalize_expression(expression, norm='srs'):
     normexp = []
     for n, exp in enumerate(expression):
         # get non-NaN values
-        notna = exp.notna().all(axis=1)
+        notna = np.logical_not(exp.isna().all(axis=1))
         data = np.asarray(exp)[notna]
 
+        if norm == 'rs':
+            normed = _rs(data, axis=0)
         if norm == 'srs':
-            normed = _srs(data, axis=0)
+            normed = _srs(data, low=0, high=1, axis=0)
+        elif norm == 'center':
+            normed = data - data.mean(axis=0, keepdims=True)
         elif norm == 'zscore':
             normed = sstats.zscore(data, axis=0, ddof=1)
+        elif norm == 'minmax':
+            normed = _rescale(data, low=0, high=1, axis=0)
         elif norm == 'batch':
             normed = corrected[n]
 
@@ -164,7 +254,7 @@ def normalize_expression(expression, norm='srs'):
         normalized.loc[notna] = normed
         normexp.append(normalized)
 
-    return normexp
+    return _unpack_tuple(normexp)
 
 
 def remove_distance(coexpression, atlas, atlas_info=None, labels=None):
@@ -326,8 +416,12 @@ def keep_stable_genes(expression, threshold=0.9, percentile=True, rank=True):
         gene_corrs[:, n] = utils.efficient_corr(for_corr[s1].loc[regions],
                                                 for_corr[s2].loc[regions])
 
-    # average similarity across donors
-    gene_corrs = gene_corrs.mean(axis=1)
+    # average similarity across donors (ignore NaNs)
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', category=RuntimeWarning,
+                                message='Mean of empty slice')
+        gene_corrs = np.nan_to_num(np.nanmean(gene_corrs, axis=1))
+
     # calculate absolute threshold if percentile is desired
     if percentile:
         threshold = np.percentile(gene_corrs, threshold * 100)
diff --git a/docs/user_guide/expression.rst b/docs/user_guide/expression.rst
index ee85161..8df7e3b 100644
--- a/docs/user_guide/expression.rst
+++ b/docs/user_guide/expression.rst
@@ -82,15 +82,15 @@ expression data normalized and aggregated across donors:
     >>> print(expression)
     gene_symbol      A1BG  A1BG-AS1       A2M  ...       ZYX     ZZEF1      ZZZ3
     label                                      ...
-    1            0.533876  0.764106  0.437582  ...  0.515566  0.211259  0.391993
-    2            0.426346  0.690736  0.438414  ...  0.418789  0.106776  0.546173
+    1            0.584804  0.680186  0.407087  ...  0.628665  0.371389  0.488942
+    2            0.518708  0.647558  0.418932  ...  0.553781  0.290940  0.564853
     3                 NaN       NaN       NaN  ...       NaN       NaN       NaN
     ...               ...       ...       ...  ...       ...       ...       ...
-    81           0.385569  0.050425  0.409866  ...  0.087828  0.607935  0.367567
-    82           0.822947  0.154579  0.371812  ...  0.140344  0.288985  0.708077
-    83           0.356487  0.027261  0.804702  ...  0.282122  0.882226  0.458286
+    81           0.495329  0.341752  0.508363  ...  0.266736  0.573451  0.490982
+    82           0.756134  0.391646  0.421420  ...  0.292576  0.379677  0.674880
+    83           0.428424  0.222963  0.660437  ...  0.411786  0.661593  0.439605
     <BLANKLINE>
-    [83 rows x 15656 columns]
+    [83 rows x 15655 columns]
 
 
 By default the data are normalized using a scaled robust sigmoid function such
@@ -159,15 +159,15 @@ will return a dense matrix (at the expense of some anatomical precision):
     >>> print(exp_exact)
     gene_symbol      A1BG  A1BG-AS1       A2M  ...       ZYX     ZZEF1      ZZZ3
     label                                      ...
-    1            0.524360  0.790998  0.445545  ...  0.529348  0.224699  0.391277
-    2            0.415996  0.715841  0.462067  ...  0.424328  0.117925  0.553627
-    3            0.648735  0.994571  0.205451  ...  0.488812  0.024174  0.718533
+    1            0.578008  0.688770  0.406520  ...  0.556295  0.374005  0.497330
+    2            0.510950  0.655648  0.420663  ...  0.481134  0.296428  0.569186
+    3            0.664877  0.818036  0.319912  ...  0.538502  0.147872  0.624574
     ...               ...       ...       ...  ...       ...       ...       ...
-    81           0.387863  0.045774  0.473089  ...  0.093229  0.596697  0.410640
-    82           0.792881  0.133295  0.372055  ...  0.143174  0.279238  0.715495
-    83           0.355565  0.022527  0.802119  ...  0.282611  0.886450  0.469348
+    81           0.486826  0.354870  0.517001  ...  0.251989  0.566561  0.498105
+    82           0.749263  0.412805  0.440311  ...  0.270990  0.376689  0.667931
+    83           0.416738  0.239668  0.661322  ...  0.389969  0.654765  0.448254
     <BLANKLINE>
-    [83 rows x 15656 columns]
+    [83 rows x 15655 columns]
 
 
 .. _usage_expression_lrmirror:
@@ -197,15 +197,15 @@ but it will dramatically increase the likelihood that this will happen:
     >>> print(exp_mirror)
     gene_symbol      A1BG  A1BG-AS1       A2M  ...       ZYX     ZZEF1      ZZZ3
     label                                      ...
-    1            0.652426  0.774026  0.359871  ...  0.517450  0.311367  0.490487
-    2            0.536944  0.616062  0.381963  ...  0.492545  0.162643  0.463119
-    3            0.917679  0.370264  0.317911  ...  0.357439  0.359224  0.935759
+    1            0.633077  0.680281  0.447366  ...  0.606646  0.378124  0.543360
+    2            0.578483  0.609816  0.478465  ...  0.597044  0.294332  0.540773
+    3            0.828558  0.643389  0.431069  ...  0.584006  0.436933  0.801877
     ...               ...       ...       ...  ...       ...       ...       ...
-    81           0.339675  0.044840  0.285005  ...  0.080648  0.608044  0.381425
-    82           0.828870  0.168068  0.418436  ...  0.103946  0.364523  0.744152
-    83           0.338746  0.021799  0.824666  ...  0.279235  0.904594  0.476128
+    81           0.478099  0.335455  0.427816  ...  0.262703  0.577574  0.492169
+    82           0.762843  0.409607  0.423332  ...  0.249643  0.389806  0.686732
+    83           0.417837  0.219333  0.676730  ...  0.404555  0.675958  0.458144
     <BLANKLINE>
-    [83 rows x 15656 columns]
+    [83 rows x 15655 columns]
 
 
 Note that since this effectively duplicates the number of tissue samples the
diff --git a/setup.cfg b/setup.cfg
index c3a3707..3a81f6b 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -61,7 +61,7 @@ console_scripts =
 
 [coverage:run]
 omit =
-    abagen/_version.py
+    */abagen/_version.py
 
 [flake8]
 doctests = True


Test patch:
diff --git a/abagen/tests/test_correct.py b/abagen/tests/test_correct.py
index a309bd3..c20d9ba 100644
--- a/abagen/tests/test_correct.py
+++ b/abagen/tests/test_correct.py
@@ -20,6 +20,15 @@ def donor_expression(testfiles, atlas):
                                      donors=['12876', '15496'])
 
 
+def test__unpack_tuple():
+    assert correct._unpack_tuple((3,)) == 3
+    assert correct._unpack_tuple((3, 3)) == (3, 3)
+    assert correct._unpack_tuple([2]) == 2
+    assert correct._unpack_tuple([2, 4]) == [2, 4]
+    assert correct._unpack_tuple(np.array([3])) == 3
+    assert np.all(correct._unpack_tuple(np.array([3, 3])) == [3, 3])
+
+
 def test__batch():
     rs = np.random.RandomState(1234)
     # p-values for ANOVA should all be ~0 (large group differences) before
@@ -39,26 +48,69 @@ def test__batch():
         correct._batch_correct([y[0]])
 
 
-def test__srs():
+def test__rescale():
+    rs = np.random.RandomState(1234)
+    y = rs.normal(size=(100, 1000)) + 10
+    out = correct._rescale(y)
+
+    # default max = 1, min =0
+    assert np.allclose(out.max(axis=0), 1) and np.allclose(out.min(axis=0), 0)
+
+    # can specify alternative min/max
+    out = correct._rescale(y, low=5, high=6)
+    assert np.allclose(out.max(axis=0), 6) and np.allclose(out.min(axis=0), 5)
+
+    # different axis works, too!
+    out = correct._rescale(y, axis=1)
+    assert np.allclose(out.max(axis=1), 1) and np.allclose(out.min(axis=1), 0)
+
+
+@pytest.mark.parametrize('a', [0, 1])
+def test__rs(a):
     rs = np.random.RandomState(1234)
 
     # create an array with a pretty ridiculous outlier effect to try and fix
     y = rs.normal(size=(100, 1000))
     y[0] += 1000
-    out = correct._srs(y)
+    y[:, 0] += 1000
+    out = correct._rs(y, axis=a)
 
-    # basic check for scaling
-    assert np.allclose(out.max(0), 1) and np.allclose(out.min(), 0)
+    # max will always be less than one, min will always be greater than zero
+    assert np.all(out.max(axis=a) <= 1) and np.all(out.min(axis=a) >= 0)
 
     # we should have reduced skewness / kurtosis compared to the original
-    assert np.all(sstats.skew(out) < sstats.skew(y))
-    assert np.all(sstats.kurtosis(out) < sstats.kurtosis(y))
+    assert np.all(sstats.skew(out, axis=a) < sstats.skew(y, axis=a))
+    assert np.all(sstats.kurtosis(out, axis=a) < sstats.kurtosis(y, axis=a))
 
     # this is a weird test; we're gonna bin the data at 0.2 intervals and make
     # sure no bins are empty. if one is something probably went wrong, right?
     for low in np.arange(0, 1, 0.2):
         hi = low + 0.2 + np.spacing(1)  # include 1
-        assert np.all(np.sum(np.logical_and(out >= low, out < hi), axis=0) > 0)
+        assert np.all(np.sum(np.logical_and(out >= low, out < hi), axis=a) > 0)
+
+
+@pytest.mark.parametrize('a', [0, 1])
+def test__srs(a):
+    rs = np.random.RandomState(1234)
+
+    # create an array with a pretty ridiculous outlier effect to try and fix
+    y = rs.normal(size=(100, 1000))
+    y[0] += 1000
+    y[:, 0] += 1000
+    out = correct._srs(y, axis=a)
+
+    # max will always be one, min will always be zero
+    assert np.allclose(out.max(axis=a), 1) and np.allclose(out.min(axis=a), 0)
+
+    # we should have reduced skewness / kurtosis compared to the original
+    assert np.all(sstats.skew(out, axis=a) < sstats.skew(y, axis=a))
+    assert np.all(sstats.kurtosis(out, axis=a) < sstats.kurtosis(y, axis=a))
+
+    # this is a weird test; we're gonna bin the data at 0.2 intervals and make
+    # sure no bins are empty. if one is something probably went wrong, right?
+    for low in np.arange(0, 1, 0.2):
+        hi = low + 0.2 + np.spacing(1)  # include 1
+        assert np.all(np.sum(np.logical_and(out >= low, out < hi), axis=a) > 0)
 
 
 def test_normalize_expression_real(testfiles):
@@ -68,27 +120,55 @@ def test_normalize_expression_real(testfiles):
     for n, idx in enumerate(inds):
         micro[n].iloc[idx] = np.nan
 
-    # min-max scaling (with some extra pizzazz)
-    srs = correct.normalize_expression(micro, norm='srs')
-    for exp, idx in zip(srs, inds):
+    # min-max scaling
+    out = correct.normalize_expression(micro, norm='minmax')
+    for exp, idx in zip(out, inds):
         assert np.all(np.isnan(exp.iloc[idx]))
         exp = exp.dropna(axis=1, how='all')
         assert np.allclose(exp.max(axis=0), 1)
         assert np.allclose(exp.min(axis=0), 0)
+    del out
+
+    # robust sigmoid
+    out = correct.normalize_expression(micro, norm='rs')
+    for exp, idx in zip(out, inds):
+        assert np.all(np.isnan(exp.iloc[idx]))
+        exp = exp.dropna(axis=1, how='all')
+        assert np.all(exp.max(axis=0) <= 1)
+        assert np.all(exp.min(axis=0) >= 0)
+    del out
+
+    # scaled robust sigmoid
+    out = correct.normalize_expression(micro, norm='srs')
+    for exp, idx in zip(out, inds):
+        assert np.all(np.isnan(exp.iloc[idx]))
+        exp = exp.dropna(axis=1, how='all')
+        assert np.allclose(exp.max(axis=0), 1)
+        assert np.allclose(exp.min(axis=0), 0)
+    del out
+
+    # centering
+    out = correct.normalize_expression(micro, norm='center')
+    for exp, idx in zip(out, inds):
+        assert np.all(np.isnan(exp.iloc[idx]))
+        exp = exp.dropna(axis=1, how='all')
+        assert np.allclose(exp.mean(axis=0), 0)
+    del out
 
     # z-scoring: mean = 0, std = 1
-    zscore = correct.normalize_expression(micro, norm='zscore')
-    for exp, idx in zip(zscore, inds):
+    out = correct.normalize_expression(micro, norm='zscore')
+    for exp, idx in zip(out, inds):
         assert np.all(np.isnan(exp.iloc[idx]))
         exp = exp.dropna(axis=1, how='all')
         assert np.allclose(exp.mean(axis=0), 0)
         assert np.allclose(exp.std(axis=0, ddof=1), 1)
+    del out
 
     # batch correct: force means identical
-    batch = correct.normalize_expression(micro, norm='batch')
-    assert np.allclose(*[e.mean(axis=0, skipna=True) for e in batch])
+    out = correct.normalize_expression(micro, norm='batch')
+    assert np.allclose(*[e.mean(axis=0, skipna=True) for e in out])
     # the NaN values should still be there, though
-    for exp, idx in zip(batch, inds):
+    for exp, idx in zip(out, inds):
         assert np.all(np.isnan(exp.iloc[idx]))
 
     # invalid norm parameter
@@ -102,6 +182,7 @@ def test_normalize_expression_real(testfiles):
 
 def test_remove_distance(donor_expression, atlas):
     expr = pd.concat(donor_expression).groupby('label').aggregate(np.mean)
+    expr = expr.dropna(axis=1, how='any')
     coexpr = np.corrcoef(expr)
     for atlas_info in [None, atlas['info']]:
         out = correct.remove_distance(coexpr, atlas['image'], atlas_info)