Repository: rmarkello/abagen
Pull request #131 (2019-11-22 21:39:05)

[ENH,REF] Adds new `region_agg` parameter
Related to #125 (insomuch as it temporarily removes `batch` normalization as an option)

This adds the `region_agg` parameter (and renames `metric` to `agg_metric`) in the `abagen.get_expression_data()` function. This parameter controls how samples belonging to the same region are combined (i.e., [1] aggregated separately for each donor and then combined across donors or [2] aggregated across all donors simultaneously).

Also some minor updates to various typos and doc-strings throughout the package.

To Do:
- [x] Update *Donor normalization* doc page to highlight `sample_norm` vs `donor_norm` differences
- [x] Add new doc page to highlight `region_agg` and `agg_metric` options
- [x] Ensure all documentation between `abagen.get_expression_data()` and `abagen` CLI is up-to-date
- [x] Add tests for new functions in `abagen.samples`
- [ ] Add tests for new parameters to `abagen.get_expression_data` and `abagen` CLI


Patch:
diff --git a/abagen/allen.py b/abagen/allen.py
index 1291ea0..7d08e59 100644
--- a/abagen/allen.py
+++ b/abagen/allen.py
@@ -16,60 +16,11 @@ lgr = logging.getLogger('abagen')
 lgr_levels = dict(zip(range(3), [40, 20, 10]))
 
 
-def groupby_label(microarray, sample_labels, labels=None, metric='mean'):
-    """
-    Averages expression data in `microarray` over samples with same label
-
-    Parameters
-    ----------
-    microarray : (S, G) pandas.DataFrame
-        Microarray expression data, where `S` is samples and `G` is genes
-    sample_labels : (S, 1) pandas.DataFrame
-        Parcel labels for `S` samples, as returned by e.g., `label_samples()`
-    labels : (L,) array_like, optional
-        All possible labels for parcellation (to account for possibility that
-        some parcels have NO expression data). Default: None
-    metric : str or func, optional
-        Mechanism by which to collapse across samples within a parcel. If a
-        str, should be in ['mean', 'median']; if a function, should be able to
-        accept an `N`-dimensional input and the `axis` keyword argument and
-        return an `N-1`-dimensional output. Default: 'mean'
-
-    Returns
-    -------
-    gene_by_label : (L, G) pandas.DataFrame
-        Microarray expression data
-    """
-
-    # get combination function
-    metric = utils.check_metric(metric)
-
-    # get missing labels
-    if labels is not None:
-        missing = np.setdiff1d(labels, sample_labels)
-        labels = pd.DataFrame(columns=microarray.columns,
-                              index=pd.Series(missing, name='label'))
-
-    gene_by_label = (microarray.merge(sample_labels,
-                                      left_index=True,
-                                      right_index=True)
-                               .groupby('label')
-                               .aggregate(metric)
-                               .append(labels)
-                               .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,
-                        sample_norm='srs', donor_norm='srs',
+def get_expression_data(atlas, atlas_info=None, *,
+                        ibf_threshold=0.5, probe_selection='diff_stability',
+                        lr_mirror=False, exact=True, tolerance=2,
+                        sample_norm='srs', gene_norm='srs',
+                        region_agg='donors', agg_metric='mean',
                         corrected_mni=True, reannotated=True,
                         return_counts=False, return_donors=False,
                         donors='all', data_dir=None, verbose=1, n_proc=1):
@@ -114,8 +65,8 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
 
     Once all samples have been matched to parcels for all supplied donors, the
     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.
+    `sample_norm` and `gene_norm` functions before being combined within
+    parcels and across donors via the supplied `agg_metric`.
 
     Parameters
     ----------
@@ -127,54 +78,63 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         `atlas`. Must have at least columns 'id', 'hemisphere', and 'structure'
         containing information mapping atlas IDs to hemisphere (i.e, "L", "R")
         and broad structural class (i.e., "cortex", "subcortex", "cerebellum").
-        Default: None
-    exact : bool, optional
-        Whether to use exact matching of donor tissue samples to parcels in
-        `atlas`. If True, this function will match tissue samples to parcels
-        within `threshold` mm of the sample; any samples that are beyond
-        `threshold` mm of a parcel will be discarded. This may result in some
-        parcels having no assigned sample / expression data. If False, the
-        default matching procedure will be performed and followed by a check
-        for parcels with no assigned samples; any such parcels will be matched
-        to the nearest sample (nearest defined as the sample with the closest
-        Euclidean distance to the parcel centroid). Default: True
-    tolerance : int, optional
-        Distance (in mm) that a sample must be from a parcel for it to be
-        matched to that parcel. This is only considered if the sample is not
-        directly within a parcel. Default: 2
-    metric : {'mean', 'median'} or callable, optional
-        Mechanism by which to reduce donor-level expression data into a single
-        dataframe. If a callable, should be able to accept an `N`-dimensional
-        input and the `axis` keyword argument and return an `N-1`-dimensional
-        output. Default: 'mean'
+        If provided, this will constrain matching of tissue samples to regions
+        in `atlas`. Default: None
     ibf_threshold : [0, 1] float, optional
-        Threshold for intensity-based filtering specifying. This number should
-        specify the ratio of samples, across all supplied donors, for which a
-        probe must have signal above background noise in order to be retained.
-        Default: 0.5
+        Threshold for intensity-based filtering. This number specifies the
+        ratio of samples, across all supplied donors, for which a probe must
+        have signal significantly greater background noise in order to be
+        retained. Default: 0.5
     probe_selection : str, optional
         Selection method for subsetting (or collapsing across) probes that
         index the same gene. Must be one of 'average', 'max_intensity',
         'max_variance', 'pc_loading', 'corr_variance', 'corr_intensity', or
-        'diff_stability'; see Notes for more information. Default:
-        'diff_stability'
+        'diff_stability'; see Notes for more information on different options.
+        Default: 'diff_stability'
     lr_mirror : bool, optional
         Whether to mirror microarray expression samples across hemispheres to
         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
+    exact : bool, optional
+        Whether to use exact matching of donor tissue samples to parcels in
+        `atlas`. If True, this function will ONLY match tissue samples to
+        parcels within `threshold` mm of the sample; any samples that are
+        beyond `threshold` mm of a parcel will be discarded. This may result
+        in some parcels having no assigned sample / expression data. If False,
+        the default matching procedure will be performed and followed by a
+        check for parcels with no assigned samples; any such parcels will be
+        matched to the nearest sample (defined as the sample with the closest
+        Euclidean distance to the parcel centroid). Default: True
+    tolerance : int, optional
+        Distance (in mm) that a sample must be from a parcel for it to be
+        matched to that parcel. This is only considered if the sample is not
+        directly within a parcel. Default: 2
     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
+        sample. Expression values are normalized separately for each sample and
+        donor across all genes; see Notes for more information on different
+        methods. If None is specified then no normalization is performed.
+        Default: 'srs'
+    gene_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 None is specified then no
-        normalization is performed. Default: 'srs'
+        donor. Expression values are normalized separately for each gene and
+        donor across all samples; see Notes for more information on different
+        methods. If None is specified then no normalization is performed.
+        Default: 'srs'
+    region_agg : {'samples', 'donors'}, optional
+        When multiple samples are identified as belonging to a region in
+        `atlas` this determines how they are aggegated. If 'samples',
+        expression data from all samples for all donors assigned to a given
+        region are combined. If 'donors', expression values for all samples
+        assigned to a given region are combined independently for each donor
+        before being combined across donors. See `agg_metric` for mechanism by
+        which samples are combined. Default: 'donors'
+    agg_metric : {'mean', 'median'} or callable, optional
+        Mechanism by which to reduce sample-level expression data into region-
+        level expression (see `region_agg`). If a callable, should be able to
+        accept an `N`-dimensional input and the `axis` keyword argument and
+        return an `N-1`-dimensional output. Default: 'mean'
     corrected_mni : bool, optional
         Whether to use the "corrected" MNI coordinates shipped with the
         `alleninf` package instead of the coordinates provided with the AHBA
@@ -185,11 +145,11 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         reannotated information will discard probes that could not be reliably
         matched to genes. Default: True
     return_counts : bool, optional
-        Whether to return how many samples were assigned to each parcel in
-        `atlas` for each donor. Default: False
+        Whether to return dataframe containing information on how many samples
+        were assigned to each parcel in `atlas` for each donor. Default: False
     return_donors : bool, optional
         Whether to return donor-level expression arrays instead of aggregating
-        expression across donors with provided `metric`. Default: False
+        expression across donors with provided `agg_metric`. Default: False
     donors : list, optional
         List of donors to use as sources of expression data. Can be either
         donor numbers or UID. If not specified will use all available donors.
@@ -302,13 +262,13 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
     # set logging verbosity level
     lgr.setLevel(lgr_levels.get(int(verbose), 1))
 
-    # load atlas_info, if provided
+    # load atlas and atlas_info, if provided
     atlas = utils.check_img(atlas)
     if atlas_info is not None:
         atlas_info = utils.check_atlas_info(atlas, atlas_info)
 
     # get combination functions
-    metric = utils.check_metric(metric)
+    agg_metric = utils.check_metric(agg_metric)
 
     # check probe_selection input
     if probe_selection not in probes.SELECTION_METHODS:
@@ -320,13 +280,8 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
     # fetch files (downloading if necessary) and unpack to variables
     files = datasets.fetch_microarray(data_dir=data_dir, donors=donors,
                                       verbose=verbose, n_proc=n_proc)
-    microarray = files['microarray']
-    annotation = files['annotation']
-    pacall = files['pacall']
-    ontology = files['ontology']
-    probe_info = files['probes'][0]
 
-    if probe_selection == 'diff_stability' and len(microarray) == 1:
+    if probe_selection == 'diff_stability' and len(files['microarray']) == 1:
         raise ValueError('Cannot use diff_stability for probe_selection with '
                          'only one donor. Please specify a different probe_'
                          'selection method or use more donors.')
@@ -336,41 +291,38 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
     if not exact:
         lgr.info('Pre-calculating centroids for {} regions in provided atlas'
                  .format(len(all_labels)))
-        centroids = utils.get_centroids(atlas, labels=all_labels)
+        centroids = utils.get_centroids(atlas, labels=all_labels,
+                                        image_space=True)
 
-    # update the annotation "files"
-    # this handles updating the MNI coordinates, dropping mistmatched samples
-    # (where MNI coordinates don't match ontology), and mirror samples across
-    # hemisphere, as applicable.
-    for n, (annot, ont) in enumerate(zip(annotation, ontology)):
-        annotation[n] = samples.update_samples(annot, ont, lr_mirror=lr_mirror,
-                                               corrected_mni=corrected_mni)
+    # update the annotation "files". this handles updating the MNI coordinates,
+    # dropping mistmatched samples (where MNI coordinates don't match the
+    # provided ontology), and mirroring samples across hemispheres, if desired
+    annotation = files['annotation']
+    ontology = files['ontology']
+    for n, (annot, ontol) in enumerate(zip(annotation, ontology)):
+        if corrected_mni:
+            annot = samples.update_mni_coords(annot)
+        annot = samples.drop_mismatch_samples(annot, ontol)
+        if lr_mirror:
+            annot = samples.mirror_samples(annot, ontol)
+        annotation[n] = annot
 
     # get dataframe of probe information (reannotated or otherwise)
-    probe_info = io.read_probes(probe_info, copy=True)
+    probe_info = io.read_probes(files['probes'][0])
     if reannotated:
-        lgr.info('Reannotating probes with information from Arnatkevic̆iūtė '
-                 'et al., 2019, NeuroImage')
         probe_info = probes.reannotate_probes(probe_info)
 
     # intensity-based filtering of probes
-    probe_info = probes.filter_probes(pacall, annotation, probe_info,
+    probe_info = probes.filter_probes(files['pacall'], annotation, probe_info,
                                       threshold=ibf_threshold)
-    lgr.info('{} probes survive intensity-based filtering with threshold of {}'
-             .format(len(probe_info), ibf_threshold))
 
     # get probe-reduced microarray expression data for all donors based on
     # selection method; this will be a list of gene x sample dataframes (one
     # for each donor)
-    lgr.info('Reducing probes indexing same gene with provided method: "{}"'
-             .format(probe_selection))
-    microarray = probes.collapse_probes(microarray, annotation,
-                                        probe_info, method=probe_selection,
-                                        inplace=True)
-    lgr.info('{} genes remain after probe filtering and selection'
-             .format(microarray[0].shape[-1]))
-
-    expression, missing = [], []
+    microarray = probes.collapse_probes(files['microarray'], annotation,
+                                        probe_info, method=probe_selection)
+
+    missing = []
     counts = pd.DataFrame(np.zeros((len(all_labels) + 1, len(microarray)),
                                    dtype=int),
                           index=np.append([0], all_labels))
@@ -382,80 +334,63 @@ def get_expression_data(atlas, atlas_info=None, *, exact=True,
         # assign samples to regions
         labels = samples.label_samples(annotation[subj], atlas,
                                        atlas_info, tolerance=tolerance)
-
         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))
+        if gene_norm is not None:
             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 #{}'
-                 .format(np.sum(np.asarray(labels) != 0),
-                         len(annotation[subj]),
-                         datasets.WELL_KNOWN_IDS.value_set('subj')[subj]))
+                                                            norm=gene_norm)
 
         # get counts of samples collapsed into each ROI
         labs, num = np.unique(labels, return_counts=True)
         counts.loc[labs, subj] = num
+        lgr.info('{:>3} / {} samples matched to regions for donor #{}'
+                 .format(counts.loc[:, subj].sum(), len(annotation[subj]),
+                         datasets.WELL_KNOWN_IDS.value_set('subj')[subj]))
 
         # if we don't want to do exact matching then cache which parcels are
         # missing data and the expression data for the closest sample to that
         # parcel; we'll use this once we've iterated through all donors
         if not exact:
-            cols = ['mni_x', 'mni_y', 'mni_z']
-            coords = utils.xyz_to_ijk(annotation[subj][cols], atlas.affine)
             empty = np.logical_not(np.in1d(all_labels, labs))
-            idx, dist = utils.closest_centroid(coords, centroids[empty],
+            cols = ['mni_x', 'mni_y', 'mni_z']
+            idx, dist = utils.closest_centroid(annotation[subj][cols],
+                                               centroids[empty],
                                                return_dist=True)
             idx = microarray[subj].loc[annotation[subj].iloc[idx].index]
             empty = all_labels[empty]
             idx.index = pd.Series(empty, name='label')
             missing += [(idx, dict(zip(empty, np.diag(dist))))]
 
-    # check for missing ROIs and fill in, as needed
-    if not exact:
-        # find labels that are missing across all donors
+        microarray[subj].index = labels['label']
+
+    if not exact:  # check for missing ROIs and fill in, as needed
+        # labels that are missing across all donors
         empty = reduce(set.intersection, [set(f.index) for f, d in missing])
         lgr.info('Matching {} regions with no data to nearest samples'
                  .format(len(empty)))
         for roi in empty:
             # find donor with sample closest to centroid of empty parcel
-            ind = np.argmin([d.get(roi) for f, d in missing])
+            ind = np.argmin([dist.get(roi) for micro, dist in missing])
             donor = datasets.WELL_KNOWN_IDS.value_set("subj")[ind]
             lgr.debug('Assigning sample from donor {} to region id #{}'
                       .format(donor, roi))
             # assign expression data from that sample and add to count
-            expression[ind].loc[roi] = missing[ind][0].loc[roi]
+            microarray[ind] = microarray[ind].append(missing[ind][0].loc[roi])
             counts.loc[roi, ind] += 1
 
-    # 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)
-        lgr.info('Dropping {} gene from concatenated expression data due to '
-                 'poor normalization'.format(drop.sum()))
-        expression = expression.drop(drop[drop].index, axis=1)
+    microarray = samples.aggregate_samples(microarray, labels=all_labels,
+                                           region_agg=region_agg,
+                                           agg_metric=agg_metric,
+                                           return_donors=return_donors)
 
     # drop the "zero" label from the counts dataframe (this is background)
     if return_counts:
-        return expression, counts.iloc[1:]
+        return microarray, counts.iloc[1:]
 
-    return expression
+    return microarray
diff --git a/abagen/cli/run.py b/abagen/cli/run.py
index c63a887..c9963c9 100644
--- a/abagen/cli/run.py
+++ b/abagen/cli/run.py
@@ -45,6 +45,8 @@ def get_parser():
     """
 
     from .. import __version__
+    from ..correct import NORMALIZATION_METHODS
+    from ..probes import SELECTION_METHODS
 
     verstr = 'abagen {}'.format(__version__)
     parser = argparse.ArgumentParser(
@@ -82,14 +84,14 @@ process:
 
 If at any step a sample can be assigned to a parcel the matching process is
 terminated. If multiple sample are assigned to the same parcel they are
-aggregated with the metric specified via the `metric` parameter. More control
-over the sample matching can be obtained by setting the `inexact` parameter;
-see the parameter description for more information.
+aggregated with the metric specified via the `agg_metric` parameter. More
+control over the sample matching can be obtained by setting the `inexact`
+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 `gene_norm` functions before being combined within parcels
+and across donors via the supplied `agg_metric`.
 """
     )
 
@@ -175,8 +177,23 @@ supplied `metric`.
                              'samples, across all supplied donors, for which '
                              'a probe must have signal above background noise '
                              'in order to be retained. Default: 0.5')
-    w_data.add_argument('--metric', action='store', default='mean',
-                        metavar='METHOD', choices=['mean', 'median'],
+    w_data.add_argument('--region_agg', '--region-agg', action='store',
+                        default='donors', metavar='METHOD',
+                        choices=['donors', 'samples'],
+                        help='When multiple samples are identified as '
+                             'belonging to a region in `atlas` this '
+                             'determines how they are aggegated. If '
+                             '\'samples\', expression data from all samples '
+                             'for all donors assigned to a given region are '
+                             'combined. If \'donors\', expression values for '
+                             'all samples assigned to a given region are '
+                             'combined independently for each donor before '
+                             'being combined across donors. See `agg_metric` '
+                             'for mechanism by which samples are combined. '
+                             'Default: \'donors\'')
+    w_data.add_argument('--agg_metric', '--agg-metric', action='store',
+                        default='mean', metavar='METHOD',
+                        choices=['mean', 'median'],
                         help='Mechanism by which to (1) reduce expression '
                              'data of multiple samples in the same `atlas` '
                              'region, and (2) reduce donor-level expression '
@@ -185,10 +202,7 @@ supplied `metric`.
                              'Default: "mean"')
     w_data.add_argument('--probe_selection', '--probe-selection',
                         action='store', default='diff_stability',
-                        metavar='METHOD',
-                        choices=['average', 'mean', 'max_intensity',
-                                 'max_variance', 'pc_loading', 'corr_variance',
-                                 'corr_intensity', 'diff_stability'],
+                        metavar='METHOD', choices=sorted(SELECTION_METHODS),
                         help='Selection method for subsetting (or collapsing '
                              'across) probes that index the same gene. Must '
                              'be one of {"average", "mean", "max_intensity", '
@@ -202,15 +216,25 @@ supplied `metric`.
                              'both hemispheres (i.e., L->R and R->L), '
                              'approximately doubling the number of available '
                              'samples. Default: False (i.e., no mirroring)')
-    w_data.add_argument('--donor_norm', '--donor-norm', action='store',
+    w_data.add_argument('--gene_norm', '--gene-norm', action='store',
                         default='srs', metavar='METHOD', type=_resolve_none,
-                        choices=['srs', 'zscore', 'batch', None],
+                        choices=sorted(NORMALIZATION_METHODS) + ['None'],
                         help='Method by which to normalize microarray '
                              'expression values for each donor prior to '
                              'collapsing across donors. Expression values are '
                              'normalized separately for each gene for each '
-                             'donor across all regions in `atlas`. Must be '
-                             'one of {"srs", "zscore", "batch", None}. '
+                             'donor across all expression samples. If None is '
+                             'specified then no normalization is performed. '
+                             'Default: "srs"')
+    w_data.add_argument('--sample_norm', '--sample-norm', action='store',
+                        default='srs', metavar='METHOD', type=_resolve_none,
+                        choices=sorted(NORMALIZATION_METHODS) + ['None'],
+                        help='Method by which to normalize microarray '
+                             'expression values for each sample prior to '
+                             'collapsing into regions in `atlas`. Expression '
+                             'values are normalized separately for each '
+                             'sample and donor across genes. If None is '
+                             'specified then no normalization is performed. '
                              'Default: "srs"')
 
     p_data = parser.add_argument_group('Options to modify the AHBA data used')
@@ -233,7 +257,7 @@ supplied `metric`.
                              'to anatomical regions. Default: False (i.e., '
                              'use corrected coordinates)')
 
-    o_data = parser.add_argument_group('Options to modify how data is output')
+    o_data = parser.add_argument_group('Options to modify how data are output')
     o_data.add_argument('--stdout', action='store_true',
                         help='Generated region x gene dataframes will be '
                              'printed to stdout for piping to other things. '
@@ -258,7 +282,7 @@ supplied `metric`.
     o_data.add_argument('--save-donors', '--save_donors', action='store_true',
                         help='Whether to save donor-level expression '
                              'dataframes instead of aggregating expression '
-                             'across donosr with provided `metric`. If '
+                             'across donors with provided `agg_metric`. If '
                              'specified, dataframes will be saved to path '
                              'specified by `output-file`, appending donor IDs '
                              'to the end of the filename. Default: False')
@@ -289,11 +313,13 @@ def main(args=None):
                                      atlas_info=opts.atlas_info,
                                      exact=opts.exact,
                                      tolerance=opts.tolerance,
-                                     metric=opts.metric,
+                                     region_agg=opts.region_agg,
+                                     agg_metric=opts.agg_metric,
                                      ibf_threshold=opts.ibf_threshold,
                                      probe_selection=opts.probe_selection,
                                      lr_mirror=opts.lr_mirror,
-                                     donor_norm=opts.donor_norm,
+                                     gene_norm=opts.gene_norm,
+                                     sample_norm=opts.sample_norm,
                                      corrected_mni=opts.corrected_mni,
                                      reannotated=opts.reannotated,
                                      return_counts=opts.save_counts,
diff --git a/abagen/correct.py b/abagen/correct.py
index 9f7d7ba..0a8b3b5 100644
--- a/abagen/correct.py
+++ b/abagen/correct.py
@@ -3,6 +3,7 @@
 Functions for processing and correcting gene expression data
 """
 
+import functools
 import itertools
 import warnings
 
@@ -156,6 +157,15 @@ def _srs(data, low=0, high=1, axis=0):
     return normed
 
 
+NORMALIZATION_METHODS = dict(
+    rs=functools.partial(_rs, axis=0),
+    srs=functools.partial(_srs, low=0, high=1, axis=0),
+    center=lambda x: x - x.mean(axis=0, keepdims=True),
+    zscore=functools.partial(sstats.zscore, axis=0, ddof=1),
+    minmax=functools.partial(_rescale, low=0, high=1, axis=0)
+)
+
+
 def normalize_expression(expression, norm='srs'):
     """
     Performs normalization on `expression` data
@@ -218,10 +228,12 @@ def normalize_expression(expression, norm='srs'):
        20130048
     """
 
-    norms = ['rs', 'srs', 'center', 'zscore', 'batch', 'minmax']
-    if norm not in norms:
+    try:
+        normfunc = NORMALIZATION_METHODS[norm]
+    except KeyError:
         raise ValueError('Provided value for `norm` not recognized. Must be '
-                         'one of {}. Received: {}'.format(norms, norm))
+                         'one of {}. Received: {}'
+                         .format(list(NORMALIZATION_METHODS), norm))
 
     # FIXME: I hate having to do this...
     if isinstance(expression, pd.DataFrame):
@@ -236,18 +248,8 @@ def normalize_expression(expression, norm='srs'):
         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, 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]
+        # normalize the data (however was specified)
+        normed = normfunc(data) if norm != 'batch' else corrected[n]
 
         # recreate dataframe and fill non-NaN values
         normalized = pd.DataFrame(np.nan, columns=exp.columns, index=exp.index)
diff --git a/abagen/io.py b/abagen/io.py
index e2b52ed..204cf6c 100644
--- a/abagen/io.py
+++ b/abagen/io.py
@@ -138,7 +138,7 @@ def read_ontology(fname, copy=False):
         data = pd.read_csv(fname)
     except ValueError:
         if not isinstance(fname, pd.DataFrame):
-            raise TypeError('Provided fname must be filepath to Ontology.csv'
+            raise TypeError('Provided fname must be filepath to Ontology.csv '
                             'file from Allen Human Brain Atlas.')
         data = fname.copy() if copy else fname
 
@@ -197,7 +197,7 @@ def read_pacall(fname, copy=False, parquet=True):
         data.columns = pd.Series(range(len(data.columns)), name='sample_id')
     except (AttributeError, ValueError):
         if not isinstance(fname, pd.DataFrame):
-            raise TypeError('Provided fname must be filepath to PACall.csv'
+            raise TypeError('Provided fname must be filepath to PACall.csv '
                             'file from Allen Human Brain Atlas.')
         data = fname.copy() if copy else fname
 
@@ -236,7 +236,7 @@ def read_probes(fname, copy=False):
         data = pd.read_csv(fname, index_col=0)
     except ValueError:
         if not isinstance(fname, pd.DataFrame):
-            raise TypeError('Provided fname must be filepath to Probes.csv'
+            raise TypeError('Provided fname must be filepath to Probes.csv '
                             'file from Allen Human Brain Atlas.')
         data = fname.copy() if copy else fname
 
@@ -278,7 +278,7 @@ def read_annotation(fname, copy=False):
         data.index.name = 'sample_id'
     except ValueError:
         if not isinstance(fname, pd.DataFrame):
-            raise TypeError('Provided fname must be filepath to Annotation.'
+            raise TypeError('Provided fname must be filepath to Annotation'
                             '.csv file from Allen Human Brain Atlas.')
         data = fname.copy() if copy else fname
 
diff --git a/abagen/probes.py b/abagen/probes.py
index ddc96bb..6d75576 100644
--- a/abagen/probes.py
+++ b/abagen/probes.py
@@ -7,6 +7,7 @@ import functools
 import gzip
 from io import StringIO
 import itertools
+import logging
 from pkg_resources import resource_filename
 
 import numpy as np
@@ -14,6 +15,8 @@ import pandas as pd
 
 from . import io, utils
 
+lgr = logging.getLogger('abagen')
+
 
 def reannotate_probes(probes):
     """
@@ -41,6 +44,9 @@ def reannotate_probes(probes):
        data. NeuroImage, 189, 353-367.
     """
 
+    lgr.info('Reannotating probes with information from Arnatkevic̆iūtė '
+             'et al., 2019, NeuroImage')
+
     # load in reannotated probes
     reannot = resource_filename('abagen', 'data/reannotated.csv.gz')
     with gzip.open(reannot, 'r') as src:
@@ -91,6 +97,9 @@ def filter_probes(pacall, samples, probes, threshold=0.5):
 
     threshold = np.clip(threshold, 0.0, 1.0)
 
+    lgr.info('Filtering probes with intensity-based threshold of {}'
+             .format(threshold))
+
     probes = io.read_probes(probes)
     signal, n_samp = np.zeros(len(probes), dtype=int), 0
     for pa, annot in zip(pacall, samples):
@@ -103,6 +112,8 @@ def filter_probes(pacall, samples, probes, threshold=0.5):
     # calculate proportion of signal to noise for given probe across samples
     keep = (signal / n_samp) >= threshold
 
+    lgr.info('{} probes survive intensity-based filtering'.format(keep.sum()))
+
     return probes[keep]
 
 
@@ -429,8 +440,7 @@ SELECTION_METHODS = dict(
 )
 
 
-def collapse_probes(microarray, annotation, probes, method='diff_stability',
-                    inplace=False):
+def collapse_probes(microarray, annotation, probes, method='diff_stability'):
     """
     Reduces `microarray` to a sample x gene expression dataframe
 
@@ -460,9 +470,6 @@ def collapse_probes(microarray, annotation, probes, method='diff_stability',
         same gene. Must be one of 'average', 'intensity', 'variance',
         'pc_loading', 'corr_variance', 'corr_intensity', or 'diff_stability';
         see Notes for more information. Default: 'diff_stability'
-    inplace : bool, optional
-        Whether to conserve memory by editing dataframes in-place instead of
-        returning edited copies. Default: False
 
     Returns
     -------
@@ -608,12 +615,18 @@ def collapse_probes(microarray, annotation, probes, method='diff_stability',
         raise ValueError('Provided method must be one of {}, not: {}'
                          .format(list(SELECTION_METHODS), method))
 
+    lgr.info('Reducing probes indexing same gene with method: "{}"'
+             .format(method))
+
     # read in microarray data for all subjects; this can be quite slow...
     probes = io.read_probes(probes)
     exp = [
-        (io.read_microarray(m, copy=not inplace)
-           .loc[probes.index, io.read_annotation(a, copy=not inplace).index])
+        io.read_microarray(m).loc[probes.index, io.read_annotation(a).index]
         for m, a in zip(microarray, annotation)
     ]
+    exp = [e.T for e in collfunc(exp, probes, annotation)]
+
+    lgr.info('{} genes remain after probe filtering and selection'
+             .format(exp[0].shape[-1]))
 
-    return [e.T for e in collfunc(exp, probes, annotation)]
+    return exp
diff --git a/abagen/samples.py b/abagen/samples.py
index ce7764b..804f972 100644
--- a/abagen/samples.py
+++ b/abagen/samples.py
@@ -3,6 +3,7 @@
 Functions for cleaning and processing the AHBA microarray dataset
 """
 
+import logging
 from pkg_resources import resource_filename
 
 from nibabel.volumeutils import Recoder
@@ -12,6 +13,8 @@ from scipy.spatial.distance import cdist
 
 from . import io, utils
 
+lgr = logging.getLogger('abagen')
+
 # AHBA structure IDs corresponding to different brain parts
 ONTOLOGY = Recoder(
     (('4008', 'cerebral cortex', 'cortex'),
@@ -279,48 +282,6 @@ def mirror_samples(annotation, ontology):
     return annotation
 
 
-def update_samples(annotation, ontology, lr_mirror=False, corrected_mni=True):
-    """
-    Updates samples in provided `annotation` file
-
-    Sample MNI coordinates are updated (if `corrected_mni`; see
-    :func:`~.update_mni_coords`), mismatched samples are dropped (where MNI
-    x-coordinate does not match hemisphere designation of structure in
-    `ontology`; see :func:`~.drop_mismatch_samples`), and samples are mirrored
-    across hemispheres (if `lr_mirror`; see :func:`mirror_samples`)
-
-    Parameters
-    ----------
-    annotation : str or pandas.DataFrame
-        List of filepaths to annotation files from Allen Brain Institute (i.e.,
-        as obtained by calling :func:`abagen.fetch_microarray` and accessing
-        the `annotation` attribute on the resulting object).
-    ontology : str or pandas.DataFrame
-        List of filepaths to ontology files from Allen Brain Institute (i.e.,
-        as obtained by calling :func:`abagen.fetch_microarray` and accessing
-        the `ontology` attribute on the resulting object).
-    lr_mirror : bool, optional
-    corrected_mni : bool, optional
-
-    Returns
-    -------
-    annotation : pandas.DataFrame
-    """
-
-    # update MNI coordinates first; this will affect all later stages
-    if corrected_mni:
-        annotation = update_mni_coords(annotation)
-
-    # drop samples where (new?) MNI coords don't match hemisphere in ontology
-    annotation = drop_mismatch_samples(annotation, ontology)
-
-    # optionally duplicate samples across hemisphere
-    if lr_mirror:
-        annotation = mirror_samples(annotation, ontology)
-
-    return annotation
-
-
 def _assign_sample(sample, atlas, sample_info=None, atlas_info=None,
                    tolerance=2):
     """
@@ -352,7 +313,7 @@ def _assign_sample(sample, atlas, sample_info=None, atlas_info=None,
     """
 
     # pull relevant info from atlas
-    label_data = np.asarray(utils.check_img(atlas).dataobj)
+    atlas = utils.check_img(atlas)
 
     # expand provided coordinates to include those w/i `tolerance` of `coords`
     # set a hard euclidean distance limit to account for different voxel sizes
@@ -360,7 +321,7 @@ def _assign_sample(sample, atlas, sample_info=None, atlas_info=None,
     coords = coords[cdist(sample, coords).squeeze() < tolerance]
 
     # grab non-zero labels for expanded coordinates
-    possible_labels = label_data[coords[:, 0], coords[:, 1], coords[:, 2]]
+    possible_labels = atlas.dataobj[coords[:, 0], coords[:, 1], coords[:, 2]]
     nz_labels = possible_labels[possible_labels.nonzero()]
     labels, counts = np.unique(nz_labels, return_counts=True)
 
@@ -418,8 +379,8 @@ def _check_label(label, sample_info, atlas_info):
 
     if label != 0:
         sample_info = sample_info[cols]
-        atlas_info = atlas_info.loc[label][cols]
-        if not np.all(sample_info.values == atlas_info.values):
+        atlas_info = atlas_info.loc[label, cols]
+        if not np.all(np.asarray(sample_info) == np.asarray(atlas_info)):
             label = 0
 
     return label
@@ -499,3 +460,119 @@ def label_samples(annotation, atlas, atlas_info=None, tolerance=2):
 
     return pd.DataFrame(labelled, dtype=int,
                         columns=['label'], index=annotation.index)
+
+
+def groupby_index(microarray, labels=None, metric='mean'):
+    """
+    Averages expression data in `microarray` over samples with same label
+
+    Parameters
+    ----------
+    microarray : (S, G) pandas.DataFrame
+        Microarray expression data, where `S` is samples and `G` is genes
+    sample_labels : (S, 1) pandas.DataFrame
+        Parcel labels for `S` samples, as returned by e.g., `label_samples()`
+    labels : (L,) array_like, optional
+        All possible labels for parcellation (to account for possibility that
+        some parcels have NO expression data). Default: None
+    metric : str or func, optional
+        Mechanism by which to collapse across samples within a parcel. If a
+        str, should be in ['mean', 'median']; if a function, should be able to
+        accept an `N`-dimensional input and the `axis` keyword argument and
+        return an `N-1`-dimensional output. Default: 'mean'
+
+    Returns
+    -------
+    gene_by_label : (L, G) pandas.DataFrame
+        Microarray expression data
+    """
+
+    # get combination function
+    metric = utils.check_metric(metric)
+
+    # get missing labels
+    if labels is not None:
+        missing = np.setdiff1d(labels, np.unique(microarray.index))
+        labels = pd.DataFrame(columns=microarray.columns,
+                              index=pd.Series(missing, name='label'))
+
+    gene_by_label = (microarray.groupby('label')
+                               .aggregate(metric)
+                               .append(labels)
+                               .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 aggregate_samples(microarray, labels=None, region_agg='donors',
+                      agg_metric='mean', return_donors=False):
+    """
+    Aggregates samples in `microarray` belonging to same regions
+
+    Parameters
+    ----------
+    microarray : list of (S, G) pandas.DataFrame
+        Microarray expression data for multiple donor, where `S` is samples and
+        `G` is genes. Index of dataframes should identify to which region each
+        sample was assigned
+    labels : array_like
+
+    region_agg : {'samples', 'donors'}, optional
+        When multiple samples are identified as belonging to a region in
+        `atlas` this determines how they are aggegated. If 'samples',
+        expression data from all samples for all donors assigned to a given
+        region are combined. If 'donors', expression values for all samples
+        assigned to a given region are combined independently for each donor
+        before being combined across donors. See `agg_metric` for mechanism by
+        which samples are combined. Default: 'donors'
+    agg_metric : {'mean', 'median'} or callable, optional
+        Mechanism by which to reduce sample-level expression data into region-
+        level expression (see `region_agg`). If a callable, should be able to
+        accept an `N`-dimensional input and the `axis` keyword argument and
+        return an `N-1`-dimensional output. Default: 'mean'
+    return_donors : bool, optional
+        Whether to return donor-level expression arrays instead of aggregating
+        expression across donors with provided `agg_metric`. Default: False
+
+    Returns
+    -------
+    microarray : (R, G) pandas.DataFrame or list of (R, G) pandas.DataFrame
+        Microarray expression data aggregated across samples within each of `R`
+        regions. If `return_donors=True` then a list of dataframes is returned
+    """
+
+    region_aggs = ['donors', 'samples']
+    if region_agg not in region_aggs:
+        raise ValueError('Provided `region_agg` {} invalid. Must be one of {}'
+                         .format(region_agg, region_aggs))
+    if region_agg == 'samples' and return_donors:
+        raise ValueError('Cannot use region_agg=\'samples\' when '
+                         '`return_donors=True`.')
+    metric = utils.check_metric(agg_metric)
+
+    lgr.info('Aggregating samples to regions with provided region_agg: {}'
+             .format(region_agg))
+
+    if region_agg == 'donors':
+        microarray = [
+            groupby_index(e, labels=labels, metric=metric) for e in microarray
+        ]
+    elif region_agg == 'samples':
+        microarray = [
+            groupby_index(pd.concat(microarray), labels=labels, metric=metric)
+        ]
+
+    if not return_donors:
+        microarray = pd.concat(microarray).groupby('label').aggregate(metric)
+        # some genes may have been poorly normalized; remove these
+        drop = microarray.dropna(axis=0, how='all').isna().any(axis=0)
+        lgr.info('Dropping {} gene from concatenated expression data due to '
+                 'poor normalization'.format(drop.sum()))
+        microarray = microarray.drop(drop[drop].index, axis=1)
+
+    return microarray
diff --git a/docs/usage.rst b/docs/usage.rst
index 079ed58..15b972c 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -17,11 +17,12 @@ abagen/issues>`_.
 
 .. toctree::
    :caption: Table of Contents
-   :maxdepth: 5
+   :maxdepth: 2
    :numbered:
 
    user_guide/download.rst
    user_guide/parcellations.rst
    user_guide/expression.rst
    user_guide/probes.rst
-   user_guide/donors.rst
+   user_guide/normalization.rst
+   user_guide/aggregation.rst
diff --git a/docs/user_guide/aggregation.rst b/docs/user_guide/aggregation.rst
new file mode 100644
index 0000000..1b9304b
--- /dev/null
+++ b/docs/user_guide/aggregation.rst
@@ -0,0 +1,38 @@
+.. _usage_aggregation:
+
+Sample aggregation options
+==========================
+
+The primary goal of :func:`abagen.get_expression_data` is to allow users to
+aggregate the ~3,700 disparate tissue samples from the Allen Human Brain Atlas
+into regions of interest defined by an atlas or parcellation file. However,
+there exist several options for exactly *how* to aggregate samples within each
+region of the specified atlas.
+
+These options are controlled via two parameters to
+:func:`abagen.get_expression_data`: ``region_agg`` and ``agg_metric``. We
+discuss both parameters and the different options available to each below.
+
+.. _usage_region_agg:
+
+The `region_agg` parameter
+--------------------------
+
+This parameter determines how samples are aggregated together to generate the
+expression values for a region. It can take two values: `'donors'` or
+`'samples'`.
+
+If set to `'donors'`, expression values for all samples assigned to a region
+are aggregated independently for each donor and *then* aggregated across
+donors. If set to `'samples'`, expression values for all samples for all
+donors assigned to a region are aggregated simultaneously.
+
+.. _usage_agg_metric:
+
+The `agg_metric` parameter
+--------------------------
+
+This parameter determines the actual metric used for aggregating samples into
+regional expression values. It can be set to any callable function (as long as
+that function accepts the keyword ``axis`` argument), but generally either
+`'mean'` (the default) or `'median'` will suffice.
diff --git a/docs/user_guide/donors.rst b/docs/user_guide/donors.rst
deleted file mode 100644
index f3f687a..0000000
--- a/docs/user_guide/donors.rst
+++ /dev/null
@@ -1,140 +0,0 @@
-.. _usage_donor_norm:
-
-Donor normalization options
-===========================
-
-The microarray expression data provided by the AHBA has been subjected to some
-`normalization procedures <help.brain-map.org/display/humanbrain/
-Documentation>`_ designed to mitigate potential differences in expression
-values between donors due to "batch effects." Despite these procedures, there
-are still some notable differences between donors present in the downloadable
-data.
-
-By default, :func:`abagen.get_expression_data` aggregates expression data
-across donors (though this can be prevented via the ``return_donors``
-parameter). Prior to aggregation, the function performs a within-donor
-normalization procedure to attempt to mitigate donor-specific effects; however,
-there are a number of ways to do that.
-
-Currently, ``abagen`` supports three options for normalizing donor data:
-
-    1. :ref:`usage_donors_srs`,
-    2. :ref:`usage_donors_zscore`, and
-    3. :ref:`usage_donors_batch`
-
-All the options have been used at various points throughout the published
-record, so while there is no "right" choice we do encourage using the default
-option (:ref:`scaled robust sigmoid <usage_donors_srs>`) due to recent work by
-Arnatkevičiūte et al., 2019 showing that it is---as the name would
-suggest---robust to outlier effects commonly observed in microarray data.
-
-We describe all the methods in detail here; these can be implemented by passing
-the ``donor_norm`` keyword argument to :func:`abagen.get_expression_data`. For
-a selection of references to published works that have used these different
-methods please see the documentation of :func:`abagen.normalize_expression`.
-
-.. _usage_donors_srs:
-
-Scaled robust sigmoid
----------------------
-
-.. code-block:: python
-
-    >>> abagen.get_expression_data(atlas['image'], donor_norm='srs')
-
-Microarray values are separately normalized for each gene across brain regions
-using a robust sigmoid function:
-
-.. math::
-
-   x_{y} = \frac{1}
-                {1 + \exp \left( \frac{-(x_{y} - \langle x \rangle)}
-                                      {\text{IQR}_{x}}
-                          \right)}
-
-where :math:`\langle x \rangle` is the median and :math:`\text{IQR}_{x}` is the
-normalized interquartile range of the microarray expression values given as:
-
-.. math::
-
-   \DeclareMathOperator\erf{erf}
-   \text{IQR}_{x} = \frac{Q_{3} - Q{1}}
-                         {2 \cdot \sqrt{2} \cdot \erf^{-1}\left(\frac{1}{2}\right)}
-            \approx \frac{Q_{3} - Q_{1}}
-                         {1.35}
-
-Once the expression values have been normalized they are rescaled to the unit
-interval with:
-
-.. math::
-
-   x_{norm} = \frac{x_{y} - \text{min}(x)}
-                   {\text{max}(x) - \text{min}(x)}
-
-Normalization is performed separately for each donor.
-
-.. _usage_donors_zscore:
-
-Z-score
--------
-
-.. code-block:: python
-
-    >>> abagen.get_expression_data(atlas['image'], donor_norm='zscore')
-
-Microarray values are separately normalized for each gene across brain regions
-using a basic z-score function:
-
-.. math::
-
-    x_{norm} = \frac{x_{y} - \bar{x}}
-                    {\sigma_{x}}
-
-where :math:`\bar{x}` is the mean and :math:`\sigma_{x}` is the sample standard
-deviation of the microarray expression values.
-
-Normalization is performed separately for each donor.
-
-.. _usage_donors_batch:
-
-Batch correction
-----------------
-
-.. code-block:: python
-
-    >>> abagen.get_expression_data(atlas['image'], donor_norm='batch')
-
-Region by gene expression matrices for each donor are vertically concatenated
-(across donors) and donor-specific indicator variables are fit to the resulting
-expression data matrix with a simple linear regression. Beta estimates for
-donor effects are estimated independently for each gene:
-
-.. math::
-
-    x = \beta_{0} + \beta_{1} I_{1} + \beta_{1} I_{1} + \ldots + \beta_{n} I_{n} \epsilon
-
-where :math:`\beta_{0}` is the intercept and :math:`I_{1}` is the indicator
-variable for a given donor. Concatenated expression data are residualized based
-on the regression fit and then unstacked into individual donor expression
-matrices:
-
-.. math::
-
-   x_{norm} = x - (\beta_{1} I_{1} + \beta_{2} I_{2} + \ldots + \beta_{n} I_{n})
-
-Note that the linear model fit includes the intercept but the intercept is not
-removed during the residualization process.
-
-Normalization is performed simultaneously for all donors.
-
-No normalization
-----------------
-
-.. code-block:: python
-
-    >>> abagen.get_expression_data(atlas['image'], donor_norm=None)
-
-Providing ``None`` to the ``donor_norm`` parameter will prevent any
-normalization procedure from being performed on the data. If you use this it
-is **strongly** encouraged that you also specify ``return_donors=True`` so that
-expression data are not aggregated across donors.
diff --git a/docs/user_guide/expression.rst b/docs/user_guide/expression.rst
index 9e70e16..aba38c4 100644
--- a/docs/user_guide/expression.rst
+++ b/docs/user_guide/expression.rst
@@ -52,17 +52,21 @@ about what's happening as it goes. However, briefly the function:
     6. Selects a representative probe amongst those probes indexing the same
        gene. This occurs by default by selecting the probe with the highest
        differential stability amongst donors; refer to parameter
-       ``probe_selection`` for more info (or see :ref:`usage_probe_selection`).
-    7. Matches tissue samples to regions in the user-specified ``atlas`` and
-       aggregates samples within each region separately for each donor. Refer
-       to parameters ``atlas``, ``atlas_info``, ``exact``, ``tolerance``, and
-       ``metric`` for more info (or see :ref:`usage_expression_exact`).
-    8. Normalizes expression values for each gene across brain regions for each
+       ``probe_selection`` for more info (or see :ref:`usage_probes`).
+    7. Matches tissue samples to regions in the user-specified ``atlas``. Refer
+       to parameters ``atlas``, ``atlas_info``, ``exact``, and ``tolerance``
+       for more info (or see :ref:`usage_expression_exact`).
+    8. Normalizes expression values for each sample across genes for each
        donor. This occurs by default using a scaled robust sigmoid
-       normalization function; refer to parameter ``donor_norm`` for more info.
-    9. Aggregates expression data across all donors. This occurs by default
-       by averaging the expression values; refer to parameters ``metric`` and
-       ``return_donors`` for more info.
+       normalization function; refere to parameter ``sample_norm`` for more
+       info.
+    9. Normalizes expression values for each gene across samples for each
+       donor. This occurs by default using a scaled robust sigmoid
+       normalization function; refer to parameter ``gene_norm`` for more info.
+    10. Aggregates samples within regions in the user-specified ``atlas`` based
+        on matches made in Step 7. By default, samples are averaged separately
+        for each donor and then averaged across donors. Refer to parameters
+        ``region_agg``, ``agg_metric``, and ``return_donors`` for more info.
 
 You can investigate all these parameters and options for modifying how the
 ``expression`` array is generated by looking at the :ref:`api_ref`.
@@ -208,7 +212,6 @@ but it will dramatically increase the likelihood that this will happen:
     [83 rows x 15661 columns]
 
 
-
 Note that since this effectively duplicates the number of tissue samples the
 function runtime will increase somewhat. Also notice how the ``lr_mirror``
 parameter changes the expression values for all the regions more dramatically
diff --git a/docs/user_guide/normalization.rst b/docs/user_guide/normalization.rst
new file mode 100644
index 0000000..5a30f59
--- /dev/null
+++ b/docs/user_guide/normalization.rst
@@ -0,0 +1,162 @@
+.. _usage_normalization:
+
+Data normalization options
+==========================
+
+The microarray expression data provided by the AHBA has been subjected to some
+`normalization procedures <help.brain-map.org/display/humanbrain/
+Documentation>`_ designed to mitigate potential differences in expression
+values between donors due to "batch effects." Despite these procedures, there
+are still some notable differences between donors present in the downloadable
+data.
+
+By default, :func:`abagen.get_expression_data` aggregates expression data
+across donors (though this can be prevented via the ``return_donors``
+parameter). Prior to aggregation, the function performs a within-donor
+normalization procedure to attempt to mitigate donor-specific effects; however,
+there are a number of ways to achieve this.
+
+Currently, ``abagen`` supports five options for normalizing data:
+
+    1. :ref:`usage_norm_center`,
+    2. :ref:`usage_norm_zscore`,
+    3. :ref:`usage_norm_minmax`,
+    4. :ref:`usage_norm_rs`, and
+    5. :ref:`usage_norm_srs`
+
+All the options have been used at various points throughout the published
+record, so while there is no "right" choice we do encourage using the default
+option (:ref:`scaled robust sigmoid <usage_norm_srs>`) due to recent work by
+Arnatkevičiūte et al., 2019 showing that it is---as the name might
+suggest---robust to outlier effects commonly observed in microarray data.
+
+We describe all the methods in detail here; these can be implemented by passing
+the ``sample_norm`` and ``gene_norm`` keyword arguments to
+:func:`abagen.get_expression_data`. For a selection of references to published
+works that have used these different methods please see the documentation of
+:func:`abagen.normalize_expression`.
+
+.. _usage_norm_sampledonor:
+
+``sample_norm`` vs ``gene_norm``
+--------------------------------
+
+Microarray expression data can be normalized in two directions:
+
+    1. Each sample can be normalized across all genes, or
+    2. Each gene can be normalized across all samples
+
+These different forms of normalization are controlled by two parameters in the
+:func:`abagen.get_expression_data` function: ``sample_norm`` and ``gene_norm``.
+Note that normalization of each sample across all genes occurs before
+normalization of each gene across all samples.
+
+Both parameters can accept the same five arguments (detailed below), and both
+are turned on by default.
+
+.. _usage_norm_center:
+
+Centering
+---------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm='center', gene_norm='center')
+
+Microarray values are centered with:
+
+.. math::
+
+    x_{norm} = x_{y} - \bar{x}
+
+where :math:`\bar{x}` is the mean of the microarray expression values.
+
+.. _usage_norm_zscore:
+
+Z-score
+-------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm='zscore', gene_norm='zscore')
+
+Microarray values are normalized using a basic z-score function:
+
+.. math::
+
+    x_{norm} = \frac{x_{y} - \bar{x}}
+                    {\sigma_{x}}
+
+where :math:`\bar{x}` is the mean and :math:`\sigma_{x}` is the sample standard
+deviation of the microarray expression values.
+
+.. _usage_norm_minmax:
+
+Min-max
+-------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm='minmax', gene_norm='minmax')
+
+Microarray values are rescaled to the unit interval with:
+
+.. math::
+
+   x_{norm} = \frac{x_{y} - \text{min}(x)}
+                   {\text{max}(x) - \text{min}(x)}
+
+.. _usage_norm_rs:
+
+Robust sigmoid
+--------------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm='rs', gene_norm='rs')
+
+Microarray values are normalized using a robust sigmoid function:
+
+.. math::
+
+   x_{y} = \frac{1}
+                {1 + \exp \left( \frac{-(x_{y} - \langle x \rangle)}
+                                      {\text{IQR}_{x}}
+                          \right)}
+
+where :math:`\langle x \rangle` is the median and :math:`\text{IQR}_{x}` is the
+normalized interquartile range of the microarray expression values given as:
+
+.. math::
+
+   \DeclareMathOperator\erf{erf}
+   \text{IQR}_{x} = \frac{Q_{3} - Q{1}}
+                         {2 \cdot \sqrt{2} \cdot \erf^{-1}\left(\frac{1}{2}\right)}
+            \approx \frac{Q_{3} - Q_{1}}
+                         {1.35}
+
+.. _usage_norm_srs:
+
+Scaled robust sigmoid
+---------------------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm='rs', gene_norm='srs')
+
+Microarray values are processed with the :ref:`robust sigmoid <usage_norm_rs>`
+function and then rescaled to the unit interval with the :ref:`min-max
+<usage_norm_minmax>` function.
+
+.. _usage_norm_none:
+
+No normalization
+----------------
+
+.. code-block:: python
+
+    >>> abagen.get_expression_data(atlas['image'], sample_norm=None, gene_norm=None)
+
+Providing ``None`` to the ``sample_norm`` and ``gene_norm`` parameters will
+prevent any normalization procedure from being performed on the data. Use this
+with caution!
diff --git a/docs/user_guide/probes.rst b/docs/user_guide/probes.rst
index 06f7cc6..8bd193f 100644
--- a/docs/user_guide/probes.rst
+++ b/docs/user_guide/probes.rst
@@ -1,4 +1,4 @@
-.. _usage_probe_selection:
+.. _usage_probes:
 
 Probe selection options
 =======================
@@ -7,7 +7,7 @@ The probes used to measure microarray expression levels in the AHBA data are
 often redundant; that is, there are frequently several probes indexing the same
 gene. Since the output of the :func:`abagen.get_expression_data` workflow is a
 region by gene dataframe, at some point we need to transition from indexing
-probe expression level to indexing gene expression levels. Effectively, this
+probe expression levels to indexing gene expression levels. Effectively, this
 means we need to select from or condense the redundant probes for each gene;
 however, there are a number of ways to do that.
 


Test patch:
diff --git a/abagen/tests/cli/test_run.py b/abagen/tests/cli/test_run.py
index fbd9268..c6b6803 100644
--- a/abagen/tests/cli/test_run.py
+++ b/abagen/tests/cli/test_run.py
@@ -44,7 +44,10 @@ def test_run_get_parser(capsys, atlas, datadir):
         parser.parse_args(['--probe-selection', 'notamethod', atlas['image']])
     assert "invalid choice: 'notamethod'" in capsys.readouterr().err
     with pytest.raises(SystemExit):
-        parser.parse_args(['--donor-norm', 'notamethod', atlas['image']])
+        parser.parse_args(['--gene-norm', 'notamethod', atlas['image']])
+    assert "invalid choice: 'notamethod'" in capsys.readouterr().err
+    with pytest.raises(SystemExit):
+        parser.parse_args(['--sample-norm', 'notamethod', atlas['image']])
     assert "invalid choice: 'notamethod'" in capsys.readouterr().err
 
     # just test every option
@@ -58,10 +61,12 @@ def test_run_get_parser(capsys, atlas, datadir):
         '--inexact',
         '--tol', '5',
         '--ibf-threshold', '0.6',
-        '--metric', 'median',
+        '--region-agg', 'donors',
+        '--agg-metric', 'median',
         '--probe-selection', 'average',
         '--lr_mirror',
-        '--donor-norm', 'srs',
+        '--gene-norm', 'srs',
+        '--sample-norm', 'srs',
         '--no-reannotated', '--no-corrected-mni',
         '--output-file', 'test.csv',
         '--save-counts', '--save-donors',
diff --git a/abagen/tests/test_correct.py b/abagen/tests/test_correct.py
index c20d9ba..4555ec3 100644
--- a/abagen/tests/test_correct.py
+++ b/abagen/tests/test_correct.py
@@ -164,20 +164,20 @@ def test_normalize_expression_real(testfiles):
         assert np.allclose(exp.std(axis=0, ddof=1), 1)
     del out
 
-    # batch correct: force means identical
-    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(out, inds):
-        assert np.all(np.isnan(exp.iloc[idx]))
+    # # batch correct: force means identical
+    # 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(out, inds):
+    #     assert np.all(np.isnan(exp.iloc[idx]))
 
     # invalid norm parameter
     with pytest.raises(ValueError):
         correct.normalize_expression(micro, norm='notanorm')
 
-    # can't do batch correction with only one donor
-    with pytest.raises(ValueError):
-        correct.normalize_expression(micro[0], norm='batch')
+    # # can't do batch correction with only one donor
+    # with pytest.raises(ValueError):
+    #     correct.normalize_expression(micro[0], norm='batch')
 
 
 def test_remove_distance(donor_expression, atlas):
diff --git a/abagen/tests/test_samples.py b/abagen/tests/test_samples.py
index 9b368aa..67192b3 100644
--- a/abagen/tests/test_samples.py
+++ b/abagen/tests/test_samples.py
@@ -138,19 +138,20 @@ def test_drop_mismatch_samples(mm_annotation, ontology):
     pd.testing.assert_frame_equal(out, expected, check_like=True)
 
 
-@pytest.mark.xfail
-def test__assign_sample():
-    assert False
-
-
-@pytest.mark.xfail
-def test__check_label():
-    assert False
-
+def test__mirror_ontology(annotation, ontology):
+    aexp = pd.DataFrame(dict(mni_x=[-10, 30, 0],
+                             structure_acronym=['S', 'Cl', 'CC'],
+                             structure_id=[4260, 4322, 9422],
+                             structure_name=['subiculum, right',
+                                             'claustrum, left',
+                                             'central canal']),
+                        index=pd.Series(range(3), name='sample_id'))
 
-@pytest.mark.xfail
-def test_label_samples():
-    assert False
+    # this function doesn't touch mni_x -- it just assumes that all the
+    # hemisphere designation are incorrect and updates the structure id, name,
+    # and acronyms accordingly
+    a = samples._mirror_ontology(annotation, ontology)
+    pd.testing.assert_frame_equal(a, aexp, check_like=True)
 
 
 def test_mirror_samples(annotation, ontology):
@@ -170,20 +171,102 @@ def test_mirror_samples(annotation, ontology):
     pd.testing.assert_frame_equal(a, aexp, check_like=True)
 
 
-def test__mirror_ontology(annotation, ontology):
-    aexp = pd.DataFrame(dict(mni_x=[-10, 30, 0],
-                             structure_acronym=['S', 'Cl', 'CC'],
-                             structure_id=[4260, 4322, 9422],
-                             structure_name=['subiculum, right',
-                                             'claustrum, left',
-                                             'central canal']),
-                        index=pd.Series(range(3), name='sample_id'))
+@pytest.mark.xfail
+def test__assign_sample(atlas):
+    atlas = atlas['image']
+    assert samples._assign_sample([[0, 0, 0]], atlas, tolerance=0) == 0
+    assert samples._assign_sample([[26, 96, 66]], atlas, tolerance=0) == 71
+    assert False
 
-    # this function doesn't touch mni_x -- it just assumes that all the
-    # hemisphere designation are incorrect and updates the structure id, name,
-    # and acronyms accordingly
-    a = samples._mirror_ontology(annotation, ontology)
-    pd.testing.assert_frame_equal(a, aexp, check_like=True)
+
+@pytest.mark.xfail
+def test__check_label():
+    assert False
+
+
+@pytest.mark.xfail
+def test_label_samples():
+    assert False
+
+
+def test_groupby_index():
+    # default usage (no params)
+    microarray = pd.DataFrame([[0., 1.], [1., 2.], [5., 6.], [0., 1.]],
+                              index=pd.Series([1, 1, 1, 2], name='label'))
+    expected = pd.DataFrame([[2., 3.], [0., 1.]],
+                            index=pd.Series([1, 2], name='label'))
+    out = samples.groupby_index(microarray)
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+    # supplying `labels` appends NaN rows to output
+    expected = pd.DataFrame([[2., 3.], [0., 1.], [np.nan, np.nan]],
+                            index=pd.Series([1, 2, 3], name='label'))
+    out = samples.groupby_index(microarray, labels=[1, 2, 3])
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+    # changing `metric` works
+    expected = pd.DataFrame([[1., 2.], [0., 1.]],
+                            index=pd.Series([1, 2], name='label'))
+    out = samples.groupby_index(microarray, metric='median')
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+
+def test_aggregate_samples():
+    m1 = pd.DataFrame([[0., 1.], [1., 2.], [5., 6.], [0., 1.]],
+                      index=pd.Series([1, 1, 1, 2], name='label'))
+    m2 = pd.DataFrame([[0., 1.], [1., 2.], [5., 6.], [10., 11.]],
+                      index=pd.Series([1, 1, 2, 3], name='label'))
+
+    # region_agg='donors'
+    expected = pd.DataFrame([[1.25, 2.25],
+                             [2.5, 3.5],
+                             [10., 11.],
+                             [np.nan, np.nan]],
+                            index=pd.Series([1, 2, 3, 4], name='label'))
+    out = samples.aggregate_samples([m1, m2], labels=[1, 2, 3, 4],
+                                    region_agg='donors', agg_metric='mean')
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+    # region_agg = 'samples'
+    expected = pd.DataFrame([[1.4, 2.4],
+                             [2.5, 3.5],
+                             [10., 11.],
+                             [np.nan, np.nan]],
+                            index=pd.Series([1, 2, 3, 4], name='label'))
+    out = samples.aggregate_samples([m1, m2], labels=[1, 2, 3, 4],
+                                    region_agg='samples', agg_metric='mean')
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+    # return_donors=True, agg_metric='median'
+    expected = [
+        pd.DataFrame([[1., 2.], [0., 1.], [np.nan, np.nan], [np.nan, np.nan]],
+                     index=pd.Series([1, 2, 3, 4], name='label')),
+        pd.DataFrame([[0.5, 1.5], [5, 6], [10, 11], [np.nan, np.nan]],
+                     index=pd.Series([1, 2, 3, 4], name='label'))
+    ]
+    out = samples.aggregate_samples([m1, m2], labels=[1, 2, 3, 4],
+                                    return_donors=True, agg_metric='median')
+    for e, o in zip(expected, out):
+        pd.testing.assert_frame_equal(o, e, check_like=True)
+
+    # check that poorly normalized genes are removed
+    m1.loc[2, 1], m2.loc[2, 1] = np.nan, np.nan
+    expected = pd.DataFrame([[1.25],
+                             [2.5],
+                             [10],
+                             [np.nan]],
+                            index=pd.Series([1, 2, 3, 4], name='label'))
+    out = samples.aggregate_samples([m1, m2], labels=[1, 2, 3, 4],
+                                    region_agg='donors', agg_metric='mean')
+    pd.testing.assert_frame_equal(out, expected, check_like=True)
+
+    # invalid method for region_agg
+    with pytest.raises(ValueError):
+        samples.aggregate_samples([m1, m2], region_agg='notamethod')
+    # region_agg='samples' incompatible with return_donors=True
+    with pytest.raises(ValueError):
+        samples.aggregate_samples([m1, m2], region_agg='samples',
+                                  return_donors=True)
 
 
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #