Repository: rmarkello/abagen
Pull request #177 (2021-03-12 21:32:03)

[ENH] Add Methods auto-generation
Adds a `return_report` option (and the `abagen.reporting` module) to auto-generate Methods sections detailing the processing performed with `abagen.get_expression_data()`.

Also adds some smoke tests and some documentation to describe the new functionality.


Patch:
diff --git a/abagen/allen.py b/abagen/allen.py
index 8cbcb44..b4bd6f5 100644
--- a/abagen/allen.py
+++ b/abagen/allen.py
@@ -4,21 +4,19 @@ Functions for mapping AHBA microarray dataset to atlases and and parcellations
 """
 
 from functools import reduce
+import logging
 import warnings
 
 import nibabel as nib
 import numpy as np
 import pandas as pd
 
-from . import correct, datasets, images, io, matching, probes_, samples_, utils
-from .datasets import WELL_KNOWN_IDS
+from . import (correct, datasets, images, io, matching, probes_, reporting,
+               samples_, utils)
 from .transforms import xyz_to_ijk
 from .utils import first_entry, flatten_dict
 
-import logging
-logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
-lgr = logging.getLogger('abagen')
-lgr_levels = dict(zip(range(3), [40, 20, 10]))
+LGR = logging.getLogger('abagen')
 
 
 def get_expression_data(atlas,
@@ -40,6 +38,7 @@ def get_expression_data(atlas,
                         reannotated=True,
                         return_counts=False,
                         return_donors=False,
+                        return_report=False,
                         donors='all',
                         data_dir=None,
                         verbose=0,
@@ -107,11 +106,12 @@ def get_expression_data(atlas,
         "B") and broad structural class (i.e., "cortex", "subcortex/brainstem",
         "cerebellum"). If provided, this will constrain matching of tissue
         samples to regions in `atlas`. If `atlas` is a tuple of GIFTI images
-        with valid label tables this will be intuited. Default: None
+        with valid label tables this will be intuited from the data. Default:
+        None
     ibf_threshold : [0, 1] float, optional
         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
+        have signal significantly greater than background noise in order to be
         retained. Default: 0.5
     probe_selection : str, optional
         Selection method for subsetting (or collapsing across) probes that
@@ -119,7 +119,7 @@ def get_expression_data(atlas,
         'max_variance', 'pc_loading', 'corr_variance', 'corr_intensity', or
         'diff_stability', 'rnaseq'; see Notes for more information on different
         options. Default: 'diff_stability'
-    donor_probes : str, optional
+    donor_probes : {'aggregate', 'independent', 'common'}, optional
         Whether specified `probe_selection` method should be performed with
         microarray data from all donors ('aggregate'), independently for each
         donor ('independent'), or based on the most common selected probe
@@ -196,6 +196,10 @@ def get_expression_data(atlas,
     return_donors : bool, optional
         Whether to return donor-level expression arrays instead of aggregating
         expression across donors with provided `agg_metric`. Default: False
+    return_report : bool, optional
+        Whether to return a string containing longform text describing the
+        processing procedures used to generate the `expression` DataFrames
+        returned by this function. 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.
@@ -226,6 +230,10 @@ def get_expression_data(atlas,
         Number of samples assigned to each of `R` regions in `atlas` for each
         of `D` donors (if multiple donors were specified); only returned if
         ``return_counts=True``.
+    report : str
+        Methods describing processing procedures implemented to generate
+        `expression`, suitable to be used in a manuscript Methods section. Only
+        returned if ``return_report=True``.
 
     Notes
     -----
@@ -321,10 +329,10 @@ def get_expression_data(atlas,
     """
 
     # set logging verbosity level
-    lgr.setLevel(lgr_levels.get(int(verbose), 1))
+    LGR.setLevel(dict(zip(range(3), [40, 20, 10])).get(int(verbose), 2))
 
     # load atlas and atlas_info, if provided, and coerce to dict
-    atlas, group_atlas = coerce_atlas_to_dict(atlas, donors, atlas_info)
+    atlas, group_atlas = images.coerce_atlas_to_dict(atlas, donors, atlas_info)
 
     # get combination functions
     agg_metric = utils.check_metric(agg_metric)
@@ -352,6 +360,10 @@ def get_expression_data(atlas,
         raise ValueError('Provided lr_mirror method is invalid, must be one '
                          f'of {mirror_opts}. Received value: \'{lr_mirror}\'')
 
+    if return_donors and region_agg == 'samples':
+        raise ValueError('Cannot return donor-level expresison data when '
+                         'region_agg parameter is set to \'samples\'.')
+
     # 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)
@@ -368,8 +380,7 @@ def get_expression_data(atlas,
     all_labels = utils.first_entry(atlas).labels
     n_gb = (8 * len(all_labels) * 30000) / (1024 ** 3)
     if n_gb > 1:
-        lgr.warning(f'Output region x gene matrix may require up to {n_gb:.2f}'
-                    'GB RAM.')
+        warnings.warn(f'Output matrix may require up to {n_gb:.2f} GB RAM')
 
     # update the annotation "files". this handles updating the MNI coordinates,
     # dropping mistmatched samples (where MNI coordinates don't match the
@@ -429,8 +440,7 @@ def get_expression_data(atlas,
         # otherwise, we'll remove the non-labelled samples after normalization
         nz = np.asarray(labels != 0).squeeze()
         if nz.sum() == 0:
-            warnings.warn(f'No samples matched to atlas for donor {subj}',
-                          stacklevel=2)
+            warnings.warn(f'No samples matched to atlas for donor {subj}')
             microarray[subj].index = labels['label']
             if not exact:
                 missing += [(pd.DataFrame(), {})]
@@ -455,7 +465,7 @@ def get_expression_data(atlas,
         # get counts of samples collapsed into each ROI
         labs, num = np.unique(labels, return_counts=True)
         counts.loc[labs, subj] = num
-        lgr.info(f'{counts.iloc[1:][subj].sum():>3} / {len(nz)} '
+        LGR.info(f'{counts.iloc[1:][subj].sum():>3} / {len(nz)} '
                  f'samples matched to regions for donor #{subj}')
 
         # if we don't want to do exact matching then cache which parcels are
@@ -491,13 +501,13 @@ def get_expression_data(atlas,
     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(f'Matching {len(empty)} region(s) with no data to the '
+        LGR.info(f'Matching {len(empty)} region(s) with no data to the '
                  'nearest tissue sample(s)')
         for roi in empty:
             # find donor with sample closest to centroid of empty parcel
             ind = np.argmin([dist.get(roi) for micro, dist in missing])
             subj = list(microarray.keys())[ind]
-            lgr.debug(f'Assigning sample from donor {subj} to region #{roi}')
+            LGR.debug(f'Assigning sample from donor {subj} to region #{roi}')
             # assign expression data from that sample and add to count
             exp = missing[ind][0].loc[roi]
             microarray[subj] = microarray[subj].append(exp)
@@ -521,11 +531,35 @@ def get_expression_data(atlas,
                                             agg_metric=agg_metric,
                                             return_donors=return_donors)
 
-    # drop the "zero" label from the counts dataframe (this is background)
+    if return_report:  # generate report
+        report = reporting.Report(atlas, atlas_info=atlas[subj].atlas_info,
+                                  ibf_threshold=ibf_threshold,
+                                  probe_selection=probe_selection,
+                                  donor_probes=donor_probes,
+                                  lr_mirror=lr_mirror, exact=exact,
+                                  tolerance=tolerance, sample_norm=sample_norm,
+                                  gene_norm=gene_norm,
+                                  norm_matched=norm_matched,
+                                  norm_structures=norm_structures,
+                                  region_agg=region_agg, agg_metric=agg_metric,
+                                  corrected_mni=corrected_mni,
+                                  reannotated=reannotated, donors=donors,
+                                  return_donors=return_donors,
+                                  data_dir=data_dir).body
+        report = report.format(n_probes=len(probe_info),
+                               n_genes=(microarray[0].shape[1] if return_donors
+                                        else microarray.shape[1]))
+
+    # pack outputs
+    out = (microarray,)
     if return_counts:
-        return microarray, counts.iloc[1:]
+        out += (counts.drop([0], axis=0),)
+    if return_report:
+        out += (report,)
+    if len(out) == 1:
+        out = out[0]
 
-    return microarray
+    return out
 
 
 def get_samples_in_mask(mask=None, **kwargs):
@@ -608,65 +642,3 @@ def get_samples_in_mask(mask=None, **kwargs):
     coords = coords[np.isin(well_id, exp.index)]
 
     return exp, coords
-
-
-def coerce_atlas_to_dict(atlas, donors, atlas_info=None, data_dir=None):
-    """
-    Coerces `atlas` to dict with keys `donors`
-
-    If already a dictionary, confirms that `atlas` has entries for all values
-    in `donors`
-
-    Parameters
-    ----------
-    atlas : niimg-like object
-        A parcellation image in MNI space, where each parcel is identified by a
-        unique integer ID
-    donors : array_like
-        Donors that should have entries in returned `atlas` dictionary
-    atlas_info : os.PathLike or pandas.DataFrame, optional
-        Filepath to or pre-loaded dataframe containing information about
-        `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/brainstem",
-        "cerebellum"). If provided, this will constrain matching of tissue
-        samples to regions in `atlas`. Default: None
-    data_dir : str, optional
-        Directory where data should be downloaded and unpacked. Only used if
-        provided `atlas` is a dictionary of surface files. Default: $HOME/
-        abagen-data
-
-    Returns
-    -------
-    atlas : dict
-        Dict where keys are `donors` and values are `atlas`. If a dict was
-        provided it is checked to ensure
-    group_atlas : bool
-        Whether one atlas was provided for all donors (True) instead of
-        donor-specific atlases (False)
-    """
-
-    donors = datasets.fetchers.check_donors(donors)
-    group_atlas = True
-
-    # FIXME: so that we're not depending on type checks so much :grimacing:
-    if isinstance(atlas, dict):
-        atlas = {
-            WELL_KNOWN_IDS.subj[donor]: images.check_atlas(atl, atlas_info,
-                                                           donor, data_dir)
-            for donor, atl in atlas.items()
-        }
-        group_atlas = False
-        missing = set(donors) - set(atlas)
-        if len(missing) > 0:
-            raise ValueError('Provided `atlas` does not have entry for all '
-                             f'requested donors. Missing donors: {donors}.')
-        lgr.info('Donor-specific atlases provided; using native coords for '
-                 'tissue samples')
-    else:
-        atlas = images.check_atlas(atlas, atlas_info)
-        atlas = {donor: atlas for donor in donors}
-        lgr.info('Group-level atlas provided; using MNI coords for '
-                 'tissue samples')
-
-    return atlas, group_atlas
diff --git a/abagen/cli/run.py b/abagen/cli/run.py
index b67b88c..33bdfbf 100644
--- a/abagen/cli/run.py
+++ b/abagen/cli/run.py
@@ -7,7 +7,7 @@ from pathlib import Path
 from typing import Iterable
 import sys
 
-lgr = logging.getLogger('abagen')
+LGR = logging.getLogger('abagen')
 
 
 def isiterable(val):
@@ -388,7 +388,7 @@ def main(args=None):
     if opts.save_counts:
         expression, counts = expression
         counts_fname = os.path.join(output_path, fname_pref + '_counts.csv')
-        lgr.info('Saving samples counts to {}'.format(counts_fname))
+        LGR.info('Saving samples counts to {}'.format(counts_fname))
         counts.to_csv(counts_fname)
 
     # determine how best to save expression output files
@@ -402,7 +402,7 @@ def main(args=None):
         for donor, exp in zip(donors, expression):
             exp_fname = os.path.join(output_path,
                                      fname_pref + '_{}.csv'.format(donor))
-            lgr.info('Saving donor {} info to {}'.format(donor, exp_fname))
+            LGR.info('Saving donor {} info to {}'.format(donor, exp_fname))
             exp.to_csv(exp_fname)
     else:
         expression.to_csv(opts.output_file)
diff --git a/abagen/data/donor_info.csv b/abagen/data/donor_info.csv
index 5d800d6..3c7c990 100644
--- a/abagen/data/donor_info.csv
+++ b/abagen/data/donor_info.csv
@@ -1,7 +1,7 @@
-donor,age,sex,ethnicity,medical_conditions,post_mortem_interval_hours
-H0351 2001,24,M,African American,History of asthma,23
-H0351 2002,39,M,African American,Possible small pituitary adenoma; single neurofibrillary tangle in entorhinal cortex,10
-H0351 1009,57,M,Caucasian,History of atherosclerotic cardiovascular disease,25.5
-H0351 1012,31,M,Caucasian,Sudden cardiac arrest; benign spindle cell proliferation and dystrophic calcification in temporal horn of lateral ventricle,17.5
-H0351 1015,49,F,Hispanic,"Splenectomy, hypothyroidism treated with Levothroid; modest numbers of hemosiderin laden macrophages noted in Virchow-Robin spaces in parietal and occipital lobes, mild arteriosclerosis",30
-H0351 1016,55,M,Caucasian,"Coronary artery atherosclerosis, prescriptions for clotting and high cholesterol",18
\ No newline at end of file
+donor,uid,age,sex,ethnicity,medical_conditions,post_mortem_interval_hours
+H0351.2001,9861,24,M,African American,History of asthma,23
+H0351.2002,10021,39,M,African American,Possible small pituitary adenoma; single neurofibrillary tangle in entorhinal cortex,10
+H0351.1009,12876,57,M,Caucasian,History of atherosclerotic cardiovascular disease,25.5
+H0351.1012,14380,31,M,Caucasian,Sudden cardiac arrest; benign spindle cell proliferation and dystrophic calcification in temporal horn of lateral ventricle,17.5
+H0351.1015,15496,49,F,Hispanic,"Splenectomy, hypothyroidism treated with Levothroid; modest numbers of hemosiderin laden macrophages noted in Virchow-Robin spaces in parietal and occipital lobes, mild arteriosclerosis",30
+H0351.1016,15697,55,M,Caucasian,"Coronary artery atherosclerosis, prescriptions for clotting and high cholesterol",18
diff --git a/abagen/datasets/__init__.py b/abagen/datasets/__init__.py
index e1e4ee6..e0140c5 100644
--- a/abagen/datasets/__init__.py
+++ b/abagen/datasets/__init__.py
@@ -6,11 +6,12 @@ dataset
 __all__ = [
     'fetch_microarray', 'fetch_raw_mri', 'fetch_desikan_killiany',
     'fetch_gene_group', 'fetch_rnaseq', 'fetch_donor_info', 'fetch_freesurfer',
-    'fetch_fsaverage5', 'fetch_fsnative', 'WELL_KNOWN_IDS', '_get_dataset_dir'
+    'fetch_fsaverage5', 'fetch_fsnative', 'check_donors', 'WELL_KNOWN_IDS',
+    '_get_dataset_dir'
 ]
 
 from .fetchers import (fetch_microarray, fetch_raw_mri, fetch_desikan_killiany,
                        fetch_gene_group, fetch_rnaseq, fetch_donor_info,
                        fetch_freesurfer, fetch_fsaverage5, fetch_fsnative,
-                       WELL_KNOWN_IDS)
+                       check_donors, WELL_KNOWN_IDS)
 from .utils import _get_dataset_dir
diff --git a/abagen/images.py b/abagen/images.py
index 9d9d3fa..7f47c53 100644
--- a/abagen/images.py
+++ b/abagen/images.py
@@ -1,16 +1,19 @@
 # -*- coding: utf-8 -*-
 
+import logging
 import os
 
 import nibabel as nib
 import numpy as np
 import pandas as pd
 
-from .datasets import fetch_fsaverage5, fetch_fsnative
+from .datasets import (check_donors, fetch_fsaverage5, fetch_fsnative,
+                       WELL_KNOWN_IDS)
 from .samples_ import ONTOLOGY
-from .utils import labeltable_to_df, load_gifti
+from .utils import labeltable_to_df, load_gifti, first_entry
 from . import matching, transforms
 
+LGR = logging.getLogger('abagen')
 BACKGROUND = [
     'unknown', 'corpuscallosum',
     'Background+FreeSurfer_Defined_Medial_Wall', '???'
@@ -411,3 +414,67 @@ def check_atlas_info(atlas_info, labels):
                          .format(valid_structures, struct_diff))
 
     return atlas_info
+
+
+def coerce_atlas_to_dict(atlas, donors, atlas_info=None, data_dir=None):
+    """
+    Coerces `atlas` to dict with keys `donors`
+
+    If already a dictionary, confirms that `atlas` has entries for all values
+    in `donors`
+
+    Parameters
+    ----------
+    atlas : niimg-like object
+        A parcellation image in MNI space, where each parcel is identified by a
+        unique integer ID
+    donors : array_like
+        Donors that should have entries in returned `atlas` dictionary
+    atlas_info : os.PathLike or pandas.DataFrame, optional
+        Filepath to or pre-loaded dataframe containing information about
+        `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/brainstem",
+        "cerebellum"). If provided, this will constrain matching of tissue
+        samples to regions in `atlas`. Default: None
+    data_dir : str, optional
+        Directory where data should be downloaded and unpacked. Only used if
+        provided `atlas` is a dictionary of surface files. Default: $HOME/
+        abagen-data
+
+    Returns
+    -------
+    atlas : dict
+        Dict where keys are `donors` and values are `atlas`. If a dict was
+        provided it is checked to ensure
+    group_atlas : bool
+        Whether one atlas was provided for all donors (True) instead of
+        donor-specific atlases (False)
+    """
+
+    donors = check_donors(donors)
+    group_atlas = True
+
+    # FIXME: so that we're not depending on type checks so much :grimacing:
+    if isinstance(atlas, dict):
+        atlas = {
+            WELL_KNOWN_IDS.subj[donor]: check_atlas(atl, atlas_info,
+                                                    donor, data_dir)
+            for donor, atl in atlas.items()
+        }
+        # if it's a group atlas they should all be the same object
+        base = first_entry(atlas)
+        group_atlas = all(base is atl for atl in atlas.values())
+        missing = set(donors) - set(atlas)
+        if len(missing) > 0:
+            raise ValueError('Provided `atlas` does not have entry for all '
+                             f'requested donors. Missing donors: {donors}.')
+        LGR.info('Donor-specific atlases provided; using native coords for '
+                 'tissue samples')
+    else:
+        atlas = check_atlas(atlas, atlas_info)
+        atlas = {donor: atlas for donor in donors}
+        LGR.info('Group-level atlas provided; using MNI coords for '
+                 'tissue samples')
+
+    return atlas, group_atlas
diff --git a/abagen/probes_.py b/abagen/probes_.py
index db05fc2..c6bd974 100644
--- a/abagen/probes_.py
+++ b/abagen/probes_.py
@@ -16,7 +16,7 @@ from scipy import stats as sstats
 
 from . import datasets, io, utils
 
-lgr = logging.getLogger('abagen')
+LGR = logging.getLogger('abagen')
 
 
 def reannotate_probes(probes):
@@ -45,7 +45,7 @@ def reannotate_probes(probes):
        data. NeuroImage, 189, 353-367.
     """
 
-    lgr.info('Reannotating probes with information from Arnatkevic̆iūtė '
+    LGR.info('Reannotating probes with information from Arnatkevic̆iūtė '
              'et al., 2019, NeuroImage')
 
     # load in reannotated probes
@@ -102,7 +102,7 @@ def filter_probes(pacall, annotation, probes, threshold=0.5):
 
     threshold = np.clip(threshold, 0.0, 1.0)
 
-    lgr.info(f'Filtering probes with intensity-based threshold of {threshold}')
+    LGR.info(f'Filtering probes with intensity-based threshold of {threshold}')
 
     probes = io.read_probes(probes)
     signal, n_samp = np.zeros(len(probes), dtype=int), 0
@@ -116,7 +116,7 @@ def filter_probes(pacall, annotation, probes, threshold=0.5):
     # calculate proportion of signal to noise for given probe across samples
     keep = (signal / n_samp) >= threshold
 
-    lgr.info(f'{keep.sum()} probes survive intensity-based filtering')
+    LGR.info(f'{keep.sum()} probes survive intensity-based filtering')
 
     return probes[keep]
 
@@ -258,6 +258,14 @@ def _diff_stability(expression, probes, annotation, *args, **kwargs):
 
 def _rnaseq(expression, probes, annotation, *args, **kwargs):
     """
+    Picks one probe to represent `expression` data for each gene in `probes`
+
+    If there are multiple probes with expression data for the same gene, this
+    function will calculate the similarity between each probes' microarray
+    expression and RNAseq expression data of the relevant gene, selecting the
+    probe with the greatest similarity to the RNAseq data. Regions are defined
+    by the "structure_id" column in `annotation`; similarity is calculated by
+    the Spearman correlation coefficient.
 
     Parameters
     ----------
@@ -707,7 +715,7 @@ def collapse_probes(microarray, annotation, probes, method='diff_stability',
         raise ValueError(f'Provided `donor_probes` "{donor_probes}" is '
                          f'invalid; must be one of {valid_probes}')
 
-    lgr.info(f'Reducing probes indexing same gene with method: {method}')
+    LGR.info(f'Reducing probes indexing same gene with method: {method}')
     # subset microarray data for pre-selected probes + samples
     # this will also left/right mirror samples, if previously requested
     probes = io.read_probes(probes)
@@ -749,6 +757,6 @@ def collapse_probes(microarray, annotation, probes, method='diff_stability',
             microarray[donor] = micro.sort_index(axis=1)
 
     n_genes = utils.first_entry(microarray).shape[-1]
-    lgr.info(f'{n_genes} genes remain after probe filtering + selection')
+    LGR.info(f'{n_genes} genes remain after probe filtering + selection')
 
     return microarray
diff --git a/abagen/reporting.py b/abagen/reporting.py
new file mode 100644
index 0000000..02dace7
--- /dev/null
+++ b/abagen/reporting.py
@@ -0,0 +1,578 @@
+# -*- coding: utf-8 -*-
+"""
+Functions for generating workflow methods reports
+
+Note: all text contained within this module is released under a CC-0 license
+(https://creativecommons.org/publicdomain/zero/1.0/).
+"""
+
+import numpy as np
+
+from . import __version__
+from .datasets import check_donors, fetch_donor_info
+from .images import coerce_atlas_to_dict
+from .utils import first_entry
+
+
+METRICS = {
+    np.mean: 'mean',
+    np.median: 'median'
+}
+REFERENCES = dict(
+    A2019N=('Arnatkevic̆iūtė, A., Fulcher, B. D., & Fornito, A. (2019). '
+            'A practical guide to linking brain-wide gene expression and '
+            'neuroimaging data. Neuroimage, 189, 353-367.'),
+    F2016P=('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.'),
+    F2013J=('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.'),
+    H2012N=('Hawrylycz, M. J., Lein, E. S., Guillozet-Bongaarts, A. L., '
+            'Shen, E. H., Ng, L., Miller, J. A., ... & Jones, A. R. '
+            '(2012). An anatomically comprehensive atlas of the adult '
+            'human brain transcriptome. Nature, 489(7416), 391-399.'),
+    H2015N=('Hawrylycz, M., Miller, J. A., Menon, V., Feng, D., '
+            'Dolbeare, T., Guillozet-Bongaarts, A. L., ... & Lein, E. '
+            '(2015). Canonical genetic signatures of the adult human '
+            'brain. Nature Neuroscience, 18(12), 1832.'),
+    P2017G=('Parkes, L., Fulcher, B. D., Yücel, M., & Fornito, A. '
+            '(2017). Transcriptional signatures of connectomic '
+            'subregions of the human striatum. Genes, Brain and '
+            'Behavior, 16(7), 647-663.'),
+    Q2002N=('Quackenbush, J. (2002). Microarray data normalization and '
+            'transformation. Nature Genetics, 32(4), 496-501.'),
+    R2018N=('Romero-Garcia, R., Whitaker, K. J., Váša, F., Seidlitz, '
+            'J., Shinn, M., Fonagy, P., ... & NSPN Consortium. (2018). '
+            'Structural covariance networks are coupled to expression '
+            'of genes enriched in supragranular layers of the human '
+            'cortex. NeuroImage, 171, 256-267.')
+)
+
+
+class Report:
+    """ Generates report of methods for :func:`abagen.get_expression_data()`
+
+    Parameters
+    ----------
+
+    """
+
+    def __init__(self, atlas, atlas_info=None, *, ibf_threshold=0.5,
+                 probe_selection='diff_stability', donor_probes='aggregate',
+                 lr_mirror=None, exact=True, tolerance=2, sample_norm='srs',
+                 gene_norm='srs', norm_matched=True, norm_structures=False,
+                 region_agg='donors', agg_metric='mean', corrected_mni=True,
+                 reannotated=True, donors='all', return_donors=False,
+                 data_dir=None):
+        atlas, self.group_atlas = \
+            coerce_atlas_to_dict(atlas, donors, atlas_info=atlas_info,
+                                 data_dir=data_dir)
+        self.atlas = first_entry(atlas)
+        self.atlas_info = self.atlas.atlas_info
+        self.ibf_threshold = ibf_threshold
+        self.probe_selection = probe_selection
+        self.donor_probes = donor_probes
+        self.lr_mirror = lr_mirror
+        self.exact = exact
+        self.tolerance = tolerance
+        self.sample_norm = sample_norm
+        self.gene_norm = gene_norm
+        self.norm_matched = norm_matched
+        self.norm_structures = norm_structures
+        self.region_agg = region_agg
+        self.agg_metric = METRICS.get(agg_metric, agg_metric)
+        self.corrected_mni = corrected_mni
+        self.reannotated = reannotated
+        self.donors = donors
+        self.return_donors = return_donors
+        self.body = self.gen_report()
+
+    def gen_report(self):
+        """ Generates body of report
+        """
+
+        report = ''
+        report += """
+        Regional microarry expression data were obtained from {n_donors}
+        post-mortem brains ({n_female} female, ages {min}--{max}, {mean:.2f}
+        +/- {std:.2f}) provided by the Allen Human Brain Atlas (AHBA,
+        https://human.brain-map.org; [H2012N]). Data were processed with the
+        abagen toolbox (version {vers}; https://github.com/rmarkello/abagen)
+        """.format(**_get_donor_demographics(self.donors), vers=__version__)
+
+        if self.atlas.volumetric and self.group_atlas:
+            report += """
+            using a {n_region}-region volumetric atlas in MNI space.<br>
+            """.format(n_region=len(self.atlas.labels))
+        elif not self.atlas.volumetric and self.group_atlas:
+            report += """
+            using a {n_region}-region surface-based atlas in MNI space.<br>
+            """.format(n_region=len(self.atlas.labels))
+        elif self.atlas.volumetric and not self.group_atlas:
+            report += """
+            using a {n_region}-region volumetric atlas, independently aligned
+            to each donor's native MRI space.<br>
+            """.format(n_region=len(self.atlas.labels))
+        else:
+            report += ".<br>"
+
+        if self.reannotated:
+            report += """
+            First, microarray probes were reannotated using data provided by
+            [A2019N]; probes not matched to a valid Entrez ID were discarded.
+            """
+        else:
+            report += """
+            First, microarray probes not matched to a valid Entrez ID were
+            discarded.
+            """
+
+        if self.ibf_threshold > 0:
+            report += """
+            Next, probes were filtered based on their expression intensity
+            relative to background noise [Q2002N], such that probes with
+            intensity less than the background in >={threshold:.2f}% of samples
+            across donors were discarded, yielding {{n_probes:,}} probes.
+            """.format(threshold=self.ibf_threshold * 100)
+
+        if self.probe_selection == 'average':
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            calculated the mean expression across probes.
+            """
+        elif self.probe_selection == 'diff_stability':
+            report += r"""
+            When multiple probes indexed the expression of the same gene, we
+            selected and used the probe with the most consistent pattern of
+            regional variation across donors (i.e., differential stability;
+            [H2015N]), calculated with:<br>
+
+            $$ \Delta_{{S}}(p) = \frac{{1}}{{\binom{{N}}{{2}}}} \,
+            \sum_{{i=1}}^{{N-1}} \sum_{{j=i+1}}^{{N}}
+            \rho[B_{{i}}(p), B_{{j}}(p)] $$<br>
+
+            where $ \rho $ is Spearman's rank correlation of the expression of
+            a single probe, p, across regions in two donors $B_{{i}}$ and
+            $B_{{j}}$, and N is the total number of donors. Here, regions
+            correspond to the structural designations provided in the ontology
+            from the AHBA.
+            """
+        elif self.probe_selection == "pc_loading":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the probe with the highest loading on the first
+            principal component taken from a decomposition of probe expression
+            across samples from all donors [P2017G].
+            """
+        elif self.probe_selection == "max_intensity":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the probe with the highest mean intensity across
+            samples.
+            """
+        elif self.probe_selection == "max_variance":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the probe with the highest variance across
+            samples.
+            """
+        elif self.probe_selection == "corr_intensity":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the expression profile of a single representative
+            probe. For genes where only two probes were available we selected
+            the probe with the highest mean intensity across samples; where
+            three or more probes were available, we calculated the correlation
+            of probe expression across samples and selected the probe with the
+            highest average correlation.
+            """
+        elif self.probe_selection == "corr_variance":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the expression profile of a single representative
+            probe. For genes where only two probes were available we selected
+            the probe with the highest variance across samples; where three or
+            more probes were available, we calculated the correlation of probe
+            expression across samples and selected the probe with the highest
+            average correlation.
+            """
+        elif self.probe_selection == "rnaseq":
+            report += """
+            When multiple probes indexed the expression of the same gene we
+            selected and used the probe with the most consistent pattern of
+            regional expression to RNA-seq data (available for two donors in
+            the AHBA). That is, we calculated the Spearman’s rank correlation
+            between each probes' microarray expression and RNA-seq expression
+            data of the corresponding gene, and selected the probe with the
+            highest correspondence. Here, regions correspond to the structural
+            designations provided in the ontology from the AHBA.
+            """
+
+        if (self.donor_probes == "aggregate" and self.probe_selection not in
+                ['average', 'diff_stability', 'rnaseq']):
+            report += """
+            The selection of probes was performed using sample expression data
+            aggregated across all donors.<br>
+            """
+        elif (self.donor_probes == "aggregate" and self.probe_selection in
+                ['average', 'diff_stability', 'rnaseq']):
+            report += "<br>"
+        elif self.donor_probes == "independent":
+            report += """
+            The selection of probes was performed independently for each donor,
+            such that the probes chosen to represent each gene could differ
+            across donors.<br>
+            """
+        elif self.donor_probes == "common":
+            report += """
+            The selection of probes was performed independently for each donor,
+            and the probe most commonly selected across all donors was chosen
+            to represent the expression of each gene for all donors.<br>
+            """
+
+        if self.corrected_mni and self.group_atlas:
+            report += """
+            The MNI coordinates of tissue samples were updated to those
+            generated via non-linear registration using the Advanced
+            Normalization Tools (ANTs; https://github.com/chrisfilo/alleninf).
+            """
+
+        if self.lr_mirror == 'bidirectional':
+            report += """
+            To increase spatial coverage, tissue samples were mirrored
+            bilaterally across the left and right hemispheres [R2018N].
+            """
+        elif self.lr_mirror == 'leftright':
+            report += """
+            To increase spatial coverage, tissue samples in the left hemisphere
+            were mirrored into the right hemisphere [R2018N].
+            """
+        elif self.lr_mirror == 'rightleft':
+            report += """
+            To increase spatial coverage, tissue samples in the right
+            hemisphere were mirrored into the left hemisphere [R2018N].
+            """
+
+        if self.tolerance == 0 and not self.atlas.volumetric:
+            report += """
+            Samples were assigned to brain regions by minimizing the Euclidean
+            distance between the {space} coordinates of each sample and the
+            nearest surface vertex. Samples where the Euclidean distance to the
+            nearest vertex was greater than the mean distance for all samples
+            belonging to that donor were excluded.
+            """.format(space='MNI' if self.group_atlas else 'native voxel')
+        elif self.tolerance == 0 and self.atlas.volumetric:
+            report += """
+            Samples were assigned to brain regions in the provided atlas only
+            if their {space} coordinates were directly within a voxel belonging
+            to a parcel.
+            """.format(space='MNI' if self.group_atlas else 'native voxel')
+        elif self.tolerance > 0 and not self.atlas.volumetric:
+            report += """
+            Samples were assigned to brain regions by minimizing the Euclidean
+            distance between the {space} coordinates of each sample and the
+            nearest surface vertex. Samples where the Euclidean distance to the
+            nearest vertex was more than {tolerance} standard deviations above
+            the mean distance for all samples belonging to that donor were
+            excluded.
+            """.format(space='MNI' if self.group_atlas else 'native voxel',
+                       tolerance=self.tolerance)
+        elif self.tolerance > 0 and self.atlas.volumetric:
+            report += """
+            Samples were assigned to brain regions in the provided atlas if
+            their {space} coordinates were within {tolerance} mm of a given
+            parcel.
+            """.format(space='MNI' if self.group_atlas else 'native voxel',
+                       tolerance=self.tolerance)
+
+        if self.atlas_info is not None:
+            report += """
+            To reduce the potential for misassignment, sample-to-region
+            matching was constrained by hemisphere and gross structural
+            divisions (i.e., cortex, subcortex/brainstem, and cerebellum, such
+            that e.g., a sample in the left cortex could only be assigned to an
+            atlas parcel in the left cortex; [A2019N]).
+            """
+
+        if not self.exact:
+            report += """
+            If a brain region was not assigned any sample based on the above
+            procedure, the sample closest to the centroid of that region was
+            selected in order to ensure that all brain regions were assigned a
+            value.
+            """
+
+        if self.norm_matched:
+            report += """
+            All tissue samples not assigned to a brain region in the provided
+            atlas were discarded.<br>
+            """
+        else:
+            report += "<br>"
+
+        if self.sample_norm is not None:
+            report += """
+            Inter-subject variation was addressed {}
+            """.format(_get_norm_procedure(self.sample_norm, 'sample_norm'))
+
+        if self.gene_norm is None:
+            pass
+        elif self.gene_norm == self.sample_norm:
+            report += """
+            Gene expression values were then normalized across tissue samples
+            using an identical procedure.
+            """
+        elif (self.gene_norm != self.sample_norm
+                and self.sample_norm is not None):
+            report += """
+            Inter-subject variation in gene expression was then addressed {}
+            """.format(_get_norm_procedure(self.gene_norm, 'gene_norm'))
+        elif self.gene_norm != self.sample_norm and self.sample_norm is None:
+            report += """
+            Inter-subject variation was addressed {}
+            """.format(_get_norm_procedure(self.gene_norm, 'gene_norm'))
+
+        if not self.norm_matched and not self.norm_structures:
+            report += """
+            All available tissue samples were used in the normalization process
+            regardless of whether they were assigned to a brain region.
+            """
+        elif not self.norm_matched and self.norm_structures:
+            report += """
+            All available tissue samples were used in the normalization process
+            regardless of whether they were matched to a brain region; however,
+            normalization was performed separately for samples in distinct
+            structural classes (i.e., cortex, subcortex/brainstem, cerebellum).
+            """
+        elif self.norm_matched and self.norm_structures:
+            report += """
+            Normalization was performed separately for samples in distinct
+            structural classes (i.e., cortex, subcortex/brainstem, cerebellum).
+            """
+
+        if not self.norm_matched and self.gene_norm is not None:
+            report += """
+            Tissue samples not matched to a brain region were discarded after
+            normalization.
+            """
+
+        if self.region_agg == 'donors' and self.agg_metric == 'mean':
+            report += """
+            Samples assigned to the same brain region were averaged separately
+            for each donor{donors}, yielding a regional expression matrix
+            {n_donors}with {n_region} rows, corresponding to brain regions,
+            and {{n_genes:,}} columns, corresponding to the retained genes.
+            """.format(donors=' and then across donors' if
+                              not self.return_donors else '',
+                       n_donors=' for each donor ' if
+                                self.return_donors else '',
+                       n_region=len(self.atlas.labels))
+        elif self.region_agg == 'samples' and self.agg_metric == 'mean':
+            report += """
+            Samples assigned to the same brain region were averaged across all
+            donors, yielding a regional expression matrix with {n_region}
+            rows, corresponding to brain regions, and {{n_genes:,}} columns,
+            corresponding to the retained genes.
+            """.format(n_region=len(self.atlas.labels))
+        elif self.region_agg == 'donors' and self.agg_metric == 'median':
+            report += """
+            The median value of samples assigned to the same brain region was
+            computed separately for each donor{donors}, yielding a regional
+            expression matrix{n_donors}with {n_region} rows, corresponding to
+            brain regions, and {{n_genes:,}} columns, corresponding to the
+            retained genes.
+            """.format(donors=' and then across donors' if
+                              not self.return_donors else '',
+                       n_donors=' for each donor ' if
+                                self.return_donors else '',
+                       n_region=len(self.atlas.labels))
+        elif self.region_agg == 'samples' and self.agg_metric == 'median':
+            report += """
+            The median value of samples assigned to the same brain region was
+            computed across donors, yielding a regional expression matrix with
+            {n_region} rows, corresponding to brain regions, and {{n_genes:,}}
+            columns, corresponding to the retained genes.
+            """.format(n_region=len(self.atlas.labels))
+
+        report += str(_add_references(report))
+        return _sanitize_text(report)
+
+
+def _sanitize_text(text):
+    """ Cleans up line and paragraph breaks in `text`
+    """
+
+    text = ' '.join(text.replace('\n', ' ').split())
+    return text.replace('<br> ', '\n\n').replace('<p>', '\n')
+
+
+def _get_donor_demographics(donors):
+    """
+    Gets demographic info about the requested `donors`
+
+    Parameters
+    ----------
+    donors : list, optional
+        List of donors. Can be either donor numbers or UID. If not specified
+        will use all available donors. Note that donors '9861' and '10021'
+        have samples from both left + right hemispheres; all other donors have
+        samples from the left hemisphere only. Default: 'all'
+
+    Returns
+    -------
+    info : dict
+        With keys ['n_donor', 'n_female','min', 'max', 'mean', 'std']. The
+        values represent features of the age distribution of requested donors.
+    """
+
+    donors = [int(i) for i in check_donors(donors)]
+    info = fetch_donor_info().set_index('uid').loc[donors]
+    age_info = info['age'].describe()
+
+    dinfo = dict(n_donors=len(info), n_female=sum(info['sex'] == 'F'))
+    dinfo.update(age_info.loc[['min', 'max', 'mean', 'std']].to_dict())
+
+    return dinfo
+
+
+def _get_norm_procedure(norm, parameter):
+    """
+    Returns methods description of the provided `norm`
+
+    Parameters
+    ----------
+    norm : str
+        Normalization procedure
+
+    Returns
+    -------
+    procedure : str
+        Reporting procedures
+    """
+
+    if parameter not in ('sample_norm', 'gene_norm'):
+        raise ValueError(f'Invalid norm parameter {parameter}')
+
+    if parameter == 'sample_norm':
+        mod = ('tissue sample expression values across genes')
+        suff = ('tissue sample across genes')
+    else:
+        mod = ('gene expression values across tissue samples')
+        suff = ('gene across tissue samples')
+
+    sigmoid = r"""
+    independently for each donor using a sigmoid function [F2016P]:<br>
+
+    $$ x_{{norm}} = \frac{{1}}{{1 + \exp(-\frac{{(x - \overline{{x}})}}
+    {{\sigma_{{x}}}})}} $$<br>
+
+    where $\bar{{x}}$ is the arithmetic mean and $\sigma_{{x}}$ is the
+    sample standard deviation of the expression of a single
+    """
+    robust_sigmoid = r"""
+    using a robust sigmoid function [F2013J]:<br>
+
+    $$ x_{{norm}} = \frac{{1}}{{1 + \exp(-\frac{{(x-\langle x \rangle)}}
+    {{\text{{IQR}}_{{x}}}})}} $$<br>
+
+    where $\langle x \rangle$ is the median and $\text{{IQR}}_{{x}}$
+    is the normalized interquartile range of the expression of a single
+    """
+    rescale = r"""
+    Normalized expression values were then rescaled to the unit interval: <br>
+
+    $$ x_{{scaled}} = \frac{{x_{{norm}} - \min(x_{{norm}})}}
+    {{\max(x_{{norm}}) - \min(x_{{norm}})}} $$<br>
+    """
+
+    procedure = ''
+    if norm in ['center', 'demean']:
+        procedure += """
+        by demeaning {mod} independently for each donor.
+        """.format(mod=mod)
+    elif norm == 'zscore':
+        procedure += """
+        by mean- and variance-normalizing (i.e., z-scoring) {mod} independently
+        for each donor.
+        """.format(mod=mod)
+    elif norm == 'minmax':
+        procedure += r"""
+        by rescaling {mod} to the unit interval independently for each donor:
+        <br>
+
+        $$ x_{{{{scaled}}}} = \frac{{{{x - \min(x)}}}}{{{{\max(x) - \min(x)}}}}
+        $$<br>
+        """.format(mod=mod)
+    elif norm in ['sigmoid', 'sig']:
+        procedure += r"""
+        by normalizing {mod} {sigmoid} {suff}.
+        """.format(mod=mod, sigmoid=sigmoid, suff=suff)
+    elif norm in ['scaled_sigmoid', 'scaled_sig']:
+        procedure += r"""
+        by normalizing {mod} {sigmoid} {suff}. {rescale}
+        """.format(mod=mod, sigmoid=sigmoid, suff=suff, rescale=rescale)
+    elif norm in ['scaled_sigmoid_quantiles', 'scaled_sig_qnt']:
+        procedure += r"""
+        by normalizing {mod} {sigmoid} {suff}, calculated using only data in
+        the 5–95th percentile range to downweight the impact of outliers.
+        {rescale}
+        """.format(mod=mod, sigmoid=sigmoid, suff=suff, rescale=rescale)
+    elif norm in ['robust_sigmoid', 'rsig', 'rs']:
+        procedure += r"""
+        by normalizing {mod} {robust_sigmoid} {suff}.
+        """.format(mod=mod, robust_sigmoid=robust_sigmoid, suff=suff)
+    elif norm in ['scaled_robust_sigmoid', 'scaled_rsig', 'srs']:
+        procedure += r"""
+        by normalizing {mod} {robust_sigmoid} {suff}. {rescale}
+        """.format(mod=mod, robust_sigmoid=robust_sigmoid, suff=suff,
+                   rescale=rescale)
+    elif norm in ['mixed_sigmoid', 'mixed_sig']:
+        procedure += r"""
+        by normalizing {mod} using a mixed sigmoid function [F2013J]:<br>
+
+        $$ x_{{{{norm}}}} = \left\{{{{\begin{{{{array}}}}{{{{r r}}}}
+        \frac{{{{1}}}}{{{{1 + \exp(-\frac{{{{(x-\overline{{{{x}}}})}}}}
+        {{{{\sigma_{{{{x}}}}}}}})}}}} ,& \text{{{{IQR}}}}_{{{{x}}}} = 0
+        \frac{{{{1}}}}{{{{1 + \exp(-\frac{{{{(x-\langle x \rangle)}}}}
+        {{{{\text{{{{IQR}}}}_{{{{x}}}}}}}})}}}} ,& \text{{{{IQR}}}}_{{{{x}}}}
+        \neq 0 \end{{{{array}}}}\right. $$<br>
+
+        where $\bar{{{{x}}}}$ is the arithmetic mean, $\sigma_{{{{x}}}}$ is the
+        sample standard deviation, $\langle x \rangle$ is the median, and
+        $\text{{{{IQR}}}}_{{{{x}}}}$ is the normalized interquartile range of
+        the expression value of a single {suff}. {rescale}
+        """.format(mod=mod, suff=suff, rescale=rescale)
+
+    return procedure
+
+
+def _add_references(report):
+    """
+    Detects references in `report` and generates list
+
+    Parameters
+    ----------
+    report : str
+        Report body
+
+    Returns
+    -------
+    references : str
+        List of references to be appended to `report`
+    """
+
+    refreport = ''
+    for ref, cite in REFERENCES.items():
+        if ref in report:
+            refreport += f'[{ref}]: {cite}<p>'
+
+    if len(refreport) > 0:
+        refreport = '<br> REFERENCES<p>----------<p>' + refreport
+
+    if refreport.endswith('<p> '):
+        refreport = refreport[:-4]
+
+    return refreport
diff --git a/abagen/samples_.py b/abagen/samples_.py
index 1dbb558..95ec777 100644
--- a/abagen/samples_.py
+++ b/abagen/samples_.py
@@ -13,7 +13,7 @@ import pandas as pd
 
 from . import io, transforms, utils
 
-lgr = logging.getLogger('abagen')
+LGR = logging.getLogger('abagen')
 
 # AHBA structure IDs corresponding to different brain parts
 ONTOLOGY = nib.volumeutils.Recoder(
@@ -429,7 +429,7 @@ def aggregate_samples(microarray, labels=None, region_agg='donors',
                          '`return_donors=True`.')
     metric = utils.check_metric(agg_metric)
 
-    lgr.info('Aggregating samples to regions with provided region_agg: {}'
+    LGR.info('Aggregating samples to regions with provided region_agg: {}'
              .format(region_agg))
 
     if region_agg == 'donors':
@@ -445,7 +445,7 @@ def aggregate_samples(microarray, labels=None, region_agg='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 '
+        LGR.info('Dropping {} gene from concatenated expression data due to '
                  'poor normalization'.format(drop.sum()))
         microarray = microarray.drop(drop[drop].index, axis=1)
 
diff --git a/docs/api.rst b/docs/api.rst
index d2320b9..cb79805 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -100,6 +100,22 @@ Reference API
 
    abagen.AtlasTree
 
+.. _ref_reporting:
+
+:mod:`abagen.reporting` - Functions for generating reports
+----------------------------------------------------------
+.. automodule:: abagen.reporting
+   :no-members:
+   :no-inherited-members:
+
+.. currentmodule:: abagen.reporting
+
+.. autosummary::
+   :template: class.rst
+   :toctree: generated/
+
+   abagen.reporting.Report
+
 .. _ref_io:
 
 :mod:`abagen.io` - Loading AHBA data files
diff --git a/docs/usage.rst b/docs/usage.rst
index 423a20d..ff47938 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -27,3 +27,4 @@ abagen/issues>`_.
    user_guide/normalization.rst
    user_guide/aggregation.rst
    user_guide/mask.rst
+   user_guide/reporting.rst
diff --git a/docs/user_guide/reporting.rst b/docs/user_guide/reporting.rst
new file mode 100644
index 0000000..0bfd748
--- /dev/null
+++ b/docs/user_guide/reporting.rst
@@ -0,0 +1,132 @@
+.. _usage_reporting:
+
+Generating reporting methods
+============================
+
+There are a lot of options and parameters to choose from when processing data
+with :func:`abagen.get_expression_data` and while we've attempted to select
+reasonable defaults, we don't want to limit your options—you are free to pick
+and choose any combination of inputs to process the AHBA data! That said, we
+also wanted to make it easy for you to report *exactly* what was done to the
+data based on the parameters you choose, so to that end we have added the
+``return_report`` parameter to :func:`abagen.get_expression_data`.
+
+When ``return_report=True``, the workflow will return an extra output. That is,
+in addition to the default regional microarray expression dataframe, a string
+will be returned that describes, in detail, all the processing that was done to
+the AHBA in the process of generating the expression matrix. We have tried to
+write this in such a way that you can simply copy-and-paste the provided text
+into the methods section of a paper, though you are of course free to edit it
+as you see fit (though if you feel edits are necessary please let us know and
+we can modify the generation more permanently!).
+
+Example report
+--------------
+
+A report can be generated with:
+
+.. code-block::
+
+    >>> expression, report = abagen.get_expression_data(atlas['image'], atlas['info'],
+    ...                                                 return_report=True)
+
+Alternatively, you can use the ``abagen.reporting`` module to generate a report
+directly without having to re-run the entire pipeline. (Note that the ``Report``
+class accepts (nearly) all the same parameters as the ``get_expression_data()``
+workflow.)
+
+.. doctest::
+
+    >>> from abagen import reporting
+    >>> generator = reporting.Report(atlas['image'], atlas['info'])
+    >>> report = generator.gen_report()
+
+The returned ``report`` (with default parameters) will look something like this
+example:
+
+    *Regional microarry expression data were obtained from 6 post-mortem brains
+    (1 female, ages 24.0--57.0, 42.50 +/- 13.38) provided by the Allen Human
+    Brain Atlas (AHBA, https://human.brain-map.org; [H2012N]). Data were
+    processed with the abagen toolbox (version 1.0;
+    https://github.com/rmarkello/abagen) using a 83-region volumetric atlas
+    in MNI space.*
+
+    *First, microarray probes were reannotated using data provided by [A2019N];
+    probes not matched to a valid Entrez ID were discarded. Next, probes were
+    filtered based on their expression intensity relative to background noise
+    [Q2002N], such that probes with intensity less than the background in
+    >=50.00% of samples across donors were discarded, yielding 31,569 probes.
+    When multiple probes indexed the expression of the same gene, we selected
+    and used the probe with the most consistent pattern of regional variation
+    across donors (i.e., differential stability; [H2015N]), calculated with:*
+
+    :math:`\Delta_{S}(p) = \frac{1}{\binom{N}{2}} \, \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} \rho[B_{i}(p), B_{j}(p)]`
+
+    *where* :math:`\rho` *is Spearman's rank correlation of the expression of a
+    single probe, p, across regions in two donor brains* :math:`B_{i}` *and*
+    :math:`B_{j}`, *and N is the total number of donors. Here, regions
+    correspond to the structural designations provided in the ontology from the
+    AHBA.*
+
+    *The MNI coordinates of tissue samples were updated to those generated via
+    non-linear registration using the Advanced Normalization Tools (ANTs;
+    https://github.com/chrisfilo/alleninf). Samples were assigned to brain
+    regions in the provided atlas if their MNI coordinates were within 2 mm of
+    a given parcel. To reduce the potential for misassignment, sample-to-region
+    matching was constrained by hemisphere and gross structural divisions
+    (i.e., cortex, subcortex/brainstem, and cerebellum, such that e.g., a
+    sample in the left cortex could only be assigned to an atlas parcel in the
+    left cortex; [A2019N]). All tissue samples not assigned to a brain region
+    in the provided atlas were discarded.*
+
+    *Inter-subject variation was addressed by normalizing tissue sample
+    expression values across genes using a robust sigmoid function [F2013J]:*
+
+    :math:`x_{norm} = \frac{1}{1 + \exp(-\frac{(x-\langle x \rangle)} {\text{IQR}_{x}})}`
+
+    *where* :math:`\langle x \rangle` *is the median and* :math:`\text{IQR}_{x}`
+    *is the normalized interquartile range of the expression of a single tissue
+    sample across genes. Normalized expression values were then rescaled to the
+    unit interval:*
+
+    :math:`x_{scaled} = \frac{x_{norm} - \min(x_{norm})} {\max(x_{norm}) - \min(x_{norm})}`
+
+    *Gene expression values were then normalized across tissue samples using an
+    identical procedure. Samples assigned to the same brain region were
+    averaged separately for each donor and then across donors, yielding a
+    regional expression matrix with 83 rows, corresponding to brain regions,
+    and 15,633 columns, corresponding to the retained genes.*
+
+    *REFERENCES* |br|
+    -\-\-\-\-\-\-\-\-\- |br|
+    *[A2019N]: Arnatkevic̆iūtė, A., Fulcher, B. D., & Fornito, A. (2019). A
+    practical guide to linking brain-wide gene expression and neuroimaging
+    data. Neuroimage, 189, 353-367.* |br|
+    *[F2013J]: 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.* |br|
+    *[H2012N]: Hawrylycz, M. J., Lein, E. S., Guillozet-Bongaarts, A. L., Shen,
+    E. H., Ng, L., Miller, J. A., ... & Jones, A. R. (2012). An anatomically
+    comprehensive atlas of the adult human brain transcriptome. Nature,
+    489(7416), 391-399.* |br|
+    *[H2015N]: Hawrylycz, M., Miller, J. A., Menon, V., Feng, D., Dolbeare,
+    T., Guillozet-Bongaarts, A. L., ... & Lein, E. (2015). Canonical genetic
+    signatures of the adult human brain. Nature Neuroscience, 18(12), 1832.*
+    |br|
+    *[Q2002N]: Quackenbush, J. (2002). Microarray data normalization and
+    transformation. Nature Genetics, 32(4), 496-501.*
+
+Note that due to text formatting limitations in Python, relevant equations used
+for e.g., normalizing the expression data will be provided in LaTeX format
+(i.e., surrounded by `$$` characters and with TeX math commands).
+
+.. important::
+    Please note that we explicitly release all text in the ``abagen.reporting``
+    `module <ref_reporting>`_ (used to generate the above-referenced reports)
+    under a `CC0 license <https://creativecommons.org/publicdomain/zero/1.0/>`_
+    such that it can be used in manuscripts without modification.
+
+.. |br| raw:: html
+
+   <br>


Test patch:
diff --git a/abagen/tests/datasets/test_fetchers.py b/abagen/tests/datasets/test_fetchers.py
index a6e08bf..7bc3352 100644
--- a/abagen/tests/datasets/test_fetchers.py
+++ b/abagen/tests/datasets/test_fetchers.py
@@ -128,12 +128,12 @@ def test_fetch_desikan_killiany():
 
 def test_fetch_donor_info():
     cols = [
-        'donor', 'age', 'sex', 'ethnicity', 'medical_conditions',
+        'donor', 'uid', 'age', 'sex', 'ethnicity', 'medical_conditions',
         'post_mortem_interval_hours'
     ]
     donor_info = fetchers.fetch_donor_info()
     assert isinstance(donor_info, pd.DataFrame)
-    assert donor_info.shape == (6, 6)
+    assert donor_info.shape == (6, 7)
     assert all(donor_info.columns == cols)
 
 
diff --git a/abagen/tests/test_allen.py b/abagen/tests/test_allen.py
index 2a224c7..f84be90 100644
--- a/abagen/tests/test_allen.py
+++ b/abagen/tests/test_allen.py
@@ -3,7 +3,6 @@
 Tests for abagen.allen module
 """
 
-import nibabel as nib
 import numpy as np
 import pandas as pd
 import pytest
@@ -109,39 +108,3 @@ def test_get_samples_in_mask(testfiles, atlas):
 
     # providing the cortical mask (atlas) should reduce # of samples returned
     assert len(allexp) > len(cortexp)
-
-
-def test_coerce_atlas_to_dict(testfiles, atlas):
-    img, info = atlas['image'], atlas['info']
-    donors = ['12876', '15496']
-
-    # test providing single atlas file
-    atl, same = allen.coerce_atlas_to_dict(img, donors, info)
-    assert same
-    assert sorted(atl.keys()) == sorted(donors)
-    imgs = list(atl.values())
-    assert all(imgs[0] is a for a in imgs[1:])
-    assert isinstance(imgs[0].atlas_info, pd.DataFrame)
-
-    # test providing pre-loaded atlas file
-    atl, same = allen.coerce_atlas_to_dict(nib.load(img), donors, info)
-    assert same
-    assert sorted(atl.keys()) == sorted(donors)
-    imgs = list(atl.values())
-    assert all(imgs[0] is a for a in imgs[1:])
-    assert isinstance(imgs[0].atlas_info, pd.DataFrame)
-
-    # test providing dictionary for atlas
-    atlas_dict = {d: img for d in donors}
-    atl, same = allen.coerce_atlas_to_dict(atlas_dict, donors)
-    assert not same
-    assert sorted(atl.keys()) == sorted(donors)
-    imgs = list(atl.values())
-    assert not any(imgs[0] is a for a in imgs[1:])
-    assert imgs[0].atlas_info is None
-
-    with pytest.raises(ValueError):
-        allen.coerce_atlas_to_dict(atlas_dict, donors + ['9861'])
-
-    with pytest.raises(TypeError):
-        allen.coerce_atlas_to_dict('notanatlas', donors)
diff --git a/abagen/tests/test_images.py b/abagen/tests/test_images.py
index 59b4f6b..07c4281 100644
--- a/abagen/tests/test_images.py
+++ b/abagen/tests/test_images.py
@@ -197,3 +197,39 @@ def test_check_atlas_info(atlas):
     bad_struct_df.loc[1, 'structure'] = 'notastructure'
     with pytest.raises(ValueError):
         images.check_atlas_info(bad_struct_df, labels)
+
+
+def test_coerce_atlas_to_dict(testfiles, atlas):
+    img, info = atlas['image'], atlas['info']
+    donors = ['12876', '15496']
+
+    # test providing single atlas file
+    atl, same = images.coerce_atlas_to_dict(img, donors, info)
+    assert same
+    assert sorted(atl.keys()) == sorted(donors)
+    imgs = list(atl.values())
+    assert all(imgs[0] is a for a in imgs[1:])
+    assert isinstance(imgs[0].atlas_info, pd.DataFrame)
+
+    # test providing pre-loaded atlas file
+    atl, same = images.coerce_atlas_to_dict(nib.load(img), donors, info)
+    assert same
+    assert sorted(atl.keys()) == sorted(donors)
+    imgs = list(atl.values())
+    assert all(imgs[0] is a for a in imgs[1:])
+    assert isinstance(imgs[0].atlas_info, pd.DataFrame)
+
+    # test providing dictionary for atlas
+    atlas_dict = {d: img for d in donors}
+    atl, same = images.coerce_atlas_to_dict(atlas_dict, donors)
+    assert not same
+    assert sorted(atl.keys()) == sorted(donors)
+    imgs = list(atl.values())
+    assert not any(imgs[0] is a for a in imgs[1:])
+    assert imgs[0].atlas_info is None
+
+    with pytest.raises(ValueError):
+        images.coerce_atlas_to_dict(atlas_dict, donors + ['9861'])
+
+    with pytest.raises(TypeError):
+        images.coerce_atlas_to_dict('notanatlas', donors)
diff --git a/abagen/tests/test_reporting.py b/abagen/tests/test_reporting.py
new file mode 100644
index 0000000..21e8ffe
--- /dev/null
+++ b/abagen/tests/test_reporting.py
@@ -0,0 +1,58 @@
+# -*- coding: utf-8 -*-
+"""
+Tests for abagen.reporting module
+"""
+
+from abagen import reporting
+
+
+def test_sanitize_text():
+    assert reporting._sanitize_text('normal text') == 'normal text'
+    assert reporting._sanitize_text('header\n\n<br> this is\n a   test ') == \
+        'header \n\nthis is a test'
+
+
+def test_get_donor_demographics():
+    out = reporting._get_donor_demographics('all')
+    keys = ('n_donors', 'n_female', 'min', 'max', 'mean', 'std')
+    assert all(key in out for key in keys)
+    assert out['n_donors'] == 6 and out['n_female'] == 1
+
+    out = reporting._get_donor_demographics(['9861', '10021'])
+    assert all(key in out for key in keys)
+    assert out['n_donors'] == 2 and out['n_female'] == 0
+
+
+def test_get_norm_procedure():
+    pass
+
+
+def test_add_references():
+    assert reporting._add_references('') == ''
+    assert reporting._add_references('no references') == ''
+    assert reporting._add_references('[MISSNO]') == ''
+    assert reporting._add_references('[A2019N]') != ''
+
+
+def test_reports(atlas, surface):
+    # not sure how else to check this beyond a simple smoke test
+    # i'm certainly not gonna recode all the logic here...
+    volreport = reporting.Report(atlas['image'], atlas['info'])
+    assert isinstance(volreport, reporting.Report)
+    assert volreport.body != ''
+    assert "83-region volumetric atlas in MNI space" in volreport.body
+
+    # check that surface returns different output than volumetric atlas
+    surfreport = reporting.Report(surface['image'], surface['info'])
+    assert isinstance(surfreport, reporting.Report)
+    assert surfreport.body != ''
+    assert "68-region surface-based atlas in MNI space" in surfreport.body
+    assert surfreport.body != volreport.body
+
+    # check a couple of specific parameter choices :man_shrugging:
+    report = reporting.Report(atlas['image'], atlas['info'],
+                              lr_mirror='bidirectional', exact=False,
+                              sample_norm='srs', gene_norm='zscore').body
+    assert 'tissue samples were mirrored' in report
+    assert 'sample closest to the centroid of that region' in report
+    assert 'variation in gene expression was then addressed' in report