Repository: apriha/lineage

Issue:
Add summary info
Add capability to get information that summarizes the data (e.g., SNP data source(s), SNP count, assembly name, chromosomes, etc.). Additionally, add summary info to the output file generated by lineage.

Discussion:

Patch:
diff --git a/README.rst b/README.rst
index a677434..fe0f127 100644
--- a/README.rst
+++ b/README.rst
@@ -179,8 +179,7 @@ mother and father for that segment.
 With that background, let's find the shared DNA between the ``User662`` and ``User663`` datasets,
 calculating the centiMorgans of shared DNA and plotting the results:
 
->>> one_chrom_shared_dna, two_chrom_shared_dna, one_chrom_shared_genes, two_chrom_shared_genes = \
-l.find_shared_dna(user662, user663, cM_threshold=0.75, snp_threshold=1100)
+>>> one_chrom_shared_dna, two_chrom_shared_dna, one_chrom_shared_genes, two_chrom_shared_genes = l.find_shared_dna(user662, user663, cM_threshold=0.75, snp_threshold=1100)
 Downloading resources/genetic_map_HapMapII_GRCh37.tar.gz
 Downloading resources/cytoBand_hg19.txt.gz
 Saving output/shared_dna_User662_User663.png
@@ -222,8 +221,7 @@ Loading resources/4584.ftdna-illumina.3483.csv.gz
 
 Now let's find the shared genes:
 
->>> one_chrom_shared_dna, two_chrom_shared_dna, one_chrom_shared_genes, two_chrom_shared_genes = \
-l.find_shared_dna(user4583, user4584, shared_genes=True)
+>>> one_chrom_shared_dna, two_chrom_shared_dna, one_chrom_shared_genes, two_chrom_shared_genes = l.find_shared_dna(user4583, user4584, shared_genes=True)
 Saving output/shared_dna_User4583_User4584.png
 Saving output/shared_dna_one_chrom_User4583_User4584.csv
 Downloading resources/knownGene_hg19.txt.gz
diff --git a/lineage/individual.py b/lineage/individual.py
index 8f6c1a3..495a193 100644
--- a/lineage/individual.py
+++ b/lineage/individual.py
@@ -18,16 +18,17 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 """
 
-import gzip
+import datetime
 import os
 import re
-import zipfile
 
 import numpy as np
 import pandas as pd
 from pandas.api.types import CategoricalDtype
 
 import lineage
+from lineage.snps import (SNPs, get_assembly_name, get_chromosomes, get_chromosomes_summary,
+                          get_snp_count)
 
 class Individual(object):
     """ Object used to represent and interact with an individual.
@@ -57,6 +58,7 @@ class Individual(object):
         self._ensembl_rest_client = ensembl_rest_client
         self._snps = None
         self._assembly = None
+        self._source = []
         self._discrepant_positions_file_count = 0
         self._discrepant_genotypes_file_count = 0
 
@@ -89,6 +91,38 @@ class Individual(object):
         else:
             return None
 
+    @property
+    def snp_count(self):
+        """ Count of SNPs loaded for this ``Individual``.
+
+        Returns
+        -------
+        int
+        """
+        return get_snp_count(self._snps)
+
+    @property
+    def chromosomes(self):
+        """ Get the chromosomes of this ``Individual``'s SNPs.
+
+        Returns
+        -------
+        list
+            list of str chromosomes (e.g., ['1', '2', '3', 'MT'], empty list if no chromosomes
+        """
+        return get_chromosomes(self._snps)
+
+    @property
+    def chromosomes_summary(self):
+        """ Summary of the chromosomes represented by this ``Individual``'s SNPs.
+
+        Returns
+        -------
+        str
+            human-readable listing of chromosomes (e.g., '1-3, MT'), empty str if no chromosomes
+        """
+        return get_chromosomes_summary(self._snps)
+
     @property
     def assembly(self):
         """ Get the assembly of this ``Individual``'s SNPs.
@@ -99,6 +133,26 @@ class Individual(object):
         """
         return self._assembly
 
+    @property
+    def assembly_name(self):
+        """ Get the name of the assembly of this ``Individual``'s SNPs.
+
+        Returns
+        -------
+        str
+        """
+        return get_assembly_name(self._assembly)
+
+    @property
+    def source(self):
+        """ Summary of the SNP data source(s) for this ``Individual``'s SNPs.
+
+        Returns
+        -------
+        str
+        """
+        return ', '.join(self._source)
+
     def load_snps(self, raw_data, discrepant_snp_positions_threshold=100,
                   discrepant_genotypes_threshold=10000):
         """ Load raw genotype data.
@@ -117,11 +171,11 @@ class Individual(object):
         if type(raw_data) is list:
             for file in raw_data:
                 print('Loading ' + os.path.relpath(file))
-                self._add_snps(self._read_raw_data(file), discrepant_snp_positions_threshold,
+                self._add_snps(SNPs(file), discrepant_snp_positions_threshold,
                                discrepant_genotypes_threshold)
         elif type(raw_data) is str:
             print('Loading ' + os.path.relpath(raw_data))
-            self._add_snps(self._read_raw_data(raw_data), discrepant_snp_positions_threshold,
+            self._add_snps(SNPs(raw_data), discrepant_snp_positions_threshold,
                            discrepant_genotypes_threshold)
         else:
             raise TypeError('invalid filetype')
@@ -131,27 +185,48 @@ class Individual(object):
 
         Returns
         -------
-        bool
-            true if SNPs saved to file in output directory
+        str
+            path to file in output directory if SNPs were saved, else empty str
         """
         if self._snps is not None:
             try:
                 if lineage.create_dir(self._output_dir):
                     output_dir = self._output_dir
                 else:
-                    return False
+                    return ''
 
                 file = os.path.join(output_dir, self.get_var_name() + '.csv')
+
+                comment = '# Generated by lineage v{}, https://github.com/apriha/lineage\n' \
+                          '# Generated at {} UTC\n' \
+                          '# Assembly: {}\n' \
+                          '# SNPs: {}\n' \
+                          '# Chromosomes: {}\n' \
+                          '# Source(s): {}\n'
+
+                comment = comment.format(lineage.__version__,
+                                         datetime.datetime.utcnow().strftime('%Y-%m-%d %H:%M:%S'),
+                                         self.assembly_name,
+                                         '{:,}'.format(self.snp_count),
+                                         self.chromosomes_summary,
+                                         self.source)
+
                 print('Saving ' + os.path.relpath(file))
-                self._snps.to_csv(file, na_rep='--', header=['chromosome', 'position', 'genotype'])
+
+                with open(file, 'w') as f:
+                    f.write(comment)
+
+                # https://stackoverflow.com/a/29233924/4727627
+                with open(file, 'a') as f:
+                    self._snps.to_csv(f, na_rep='--', header=['chromosome', 'position', 'genotype'])
+
+                return file
             except Exception as err:
                 print(err)
-                return False
+                return ''
         else:
             print('no SNPs to save...')
-            return False
-
-        return True
+            return ''
 
     def get_var_name(self):
         """ Clean a string so that it can be a valid Python variable
@@ -329,197 +404,37 @@ class Individual(object):
 
         return complement
 
-    def _read_raw_data(self, file):
-        if not os.path.exists(file):
-            print(file + ' does not exist; skipping')
-            return None
-
-        # peek into files to determine the data format
-        if '.zip' in file:
-            with zipfile.ZipFile(file) as z:
-                with z.open(z.namelist()[0], 'r') as f:
-                    # https://stackoverflow.com/a/606199
-                    line = f.readline().decode('utf-8')
-        elif '.gz' in file:
-            with gzip.open(file, 'rt') as f:
-                line = f.readline()
-        else:
-            with open(file, 'r') as f:
-                line = f.readline()
-
-        if '23andMe' in line:
-            return self._read_23andme(file)
-        elif 'Ancestry' in line:
-            return self._read_ancestry(file)
-        elif line[:4] == 'RSID':
-            return self._read_ftdna(file)
-        elif line[:4] == 'rsid':
-            return self._read_generic_csv(file)
-        else:
-            return None
-
-    @staticmethod
-    def _read_23andme(file):
-        """ Read and parse 23andMe file.
-
-        https://www.23andme.com
-
-        Parameters
-        ----------
-        file : str
-            path to file
-
-        Returns
-        -------
-        pandas.DataFrame
-            individual's genetic data normalized for use with `lineage`
-        """
-        try:
-            return pd.read_csv(file, comment='#', sep='\t', na_values='--',
-                               names=['rsid', 'chrom', 'pos', 'genotype'],
-                               index_col=0, dtype={'chrom': object})
-        except Exception as err:
-            print(err)
-            return None
-
-    @staticmethod
-    def _read_ftdna(file):
-        """ Read and parse Family Tree DNA (FTDNA) file.
-
-        https://www.familytreedna.com
-
-        Parameters
-        ----------
-        file : str
-            path to file
-
-        Returns
-        -------
-        pandas.DataFrame
-            individual's genetic data normalized for use with `lineage`
-        """
-        try:
-            df = pd.read_csv(file, skiprows=1, na_values='--',
-                             names=['rsid', 'chrom', 'pos', 'genotype'],
-                             index_col=0, dtype={'chrom': object})
-
-            # remove incongruous data
-            df = df.drop(df.loc[df['chrom'] == '0'].index)
-            df = df.drop(df.loc[df.index == 'RSID'].index)  # second header for concatenated data
-
-            # if second header existed, pos dtype will be object (should be np.int64)
-            df['pos'] = df['pos'].astype(np.int64)
-
-            return df
-        except Exception as err:
-            print(err)
-            return None
-
-    @staticmethod
-    def _read_ancestry(file):
-        """ Read and parse Ancestry.com file.
-
-        http://www.ancestry.com
-
-        Parameters
-        ----------
-        file : str
-            path to file
-
-        Returns
-        -------
-        pandas.DataFrame
-            individual's genetic data normalized for use with `lineage`
-        """
-        try:
-            df = pd.read_csv(file, comment='#', header=0, sep='\t', na_values=0,
-                             names=['rsid', 'chrom', 'pos', 'allele1', 'allele2'],
-                             index_col=0, dtype={'chrom': object})
-
-            # create genotype column from allele columns
-            df['genotype'] = df['allele1'] + df['allele2']
-
-            # delete allele columns
-            # http://stackoverflow.com/a/13485766
-            del df['allele1']
-            del df['allele2']
-
-            # https://redd.it/5y90un
-            df.ix[np.where(df['chrom'] == '23')[0], 'chrom'] = 'X'
-            df.ix[np.where(df['chrom'] == '24')[0], 'chrom'] = 'Y'
-            df.ix[np.where(df['chrom'] == '25')[0], 'chrom'] = 'PAR'
-            df.ix[np.where(df['chrom'] == '26')[0], 'chrom'] = 'MT'
-
-            return df
-        except Exception as err:
-            print(err)
-            return None
-
-    @staticmethod
-    def _read_generic_csv(file):
-        """ Read and parse generic CSV file.
-
-        Notes
-        -----
-        Assumes columns are 'rsid', 'chrom' / 'chromosome', 'pos' / 'position', and 'genotype';
-        values are comma separated; unreported genotypes are indicated by '--'; and one header row
-        precedes data. For example:
-
-            rsid,chromosome,position,genotype
-            rs1,1,1,AA
-            rs2,1,2,CC
-            rs3,1,3,--
-
-        Parameters
-        ----------
-        file : str
-            path to file
-
-        Returns
-        -------
-        pandas.DataFrame
-            individual's genetic data normalized for use with `lineage`
-        """
-        try:
-            df = pd.read_csv(file, skiprows=1, na_values='--',
-                             names=['rsid', 'chrom', 'pos', 'genotype'],
-                             index_col=0, dtype={'chrom': object, 'pos': np.int64})
-
-            return df
-        except Exception as err:
-            print(err)
-            return None
-
     def _add_snps(self, snps, discrepant_snp_positions_threshold, discrepant_genotypes_threshold):
         """ Add SNPs to this Individual.
 
         Parameters
         ----------
-        snps : pandas.DataFrame
+        snps : SNPs
             SNPs to add
         discrepant_snp_positions_threshold : int
             see above
         discrepant_genotypes_threshold : int
             see above
         """
-        if snps is None:
+        if snps.snps is None:
             return
 
-        # ensure there area always two X alleles
-        snps = self._double_single_alleles(snps, 'X')
-
-        assembly = self._detect_assembly(snps)
+        assembly = snps.assembly
+        source = snps.source
 
-        if assembly is None:
-            print('assembly not detected, assuming build 37')
-            assembly = 37
+        if not snps.assembly_detected:
+            print('assembly not detected, assuming build {}'.format(snps.assembly))
 
         if self._assembly is None:
             self._assembly = assembly
         elif self._assembly != assembly:
             print('assembly / build mismatch between current assembly of SNPs and SNPs being loaded')
 
+        # ensure there area always two X alleles
+        snps = self._double_single_alleles(snps.snps, 'X')
+
         if self._snps is None:
+            self._source.append(source)
             self._snps = snps
         else:
             common_snps = self._snps.join(snps, how='inner', rsuffix='_added')
@@ -571,6 +486,7 @@ class Individual(object):
                 return
 
             # add new SNPs
+            self._source.append(source)
             self._snps = self._snps.combine_first(snps)
             self._snps.loc[discrepant_genotypes.index, 'genotype'] = np.nan
 
@@ -634,61 +550,3 @@ class Individual(object):
     def _natural_sort_key(s, natural_sort_re=re.compile('([0-9]+)')):
         return [int(text) if text.isdigit() else text.lower()
                 for text in re.split(natural_sort_re, s)]
-
-    def _detect_assembly(self, snps):
-        """ Detect assembly of SNPs being loaded.
-
-        Use the coordinates of common SNPs to identify the assembly / build of a genotype file
-        that is being loaded.
-
-        Notes
-        -----
-        rs3094315 : plus strand in 36, 37, and 38
-        rs11928389 : plus strand in 36, minus strand in 37 and 38
-        rs2500347 : plus strand in 36 and 37, minus strand in 38
-        rs964481 : plus strand in 36, 37, and 38
-
-        Parameters
-        ----------
-        snps : pandas.DataFrame
-            SNPs to add
-
-        Returns
-        -------
-        int
-            detected assembly of SNPs, else None
-
-        References
-        ----------
-        ..[1] Yates et. al. (doi:10.1093/bioinformatics/btu613),
-          http://europepmc.org/search/?query=DOI:10.1093/bioinformatics/btu613
-        ..[2] Zerbino et. al. (doi.org/10.1093/nar/gkx1098), https://doi.org/10.1093/nar/gkx1098
-        ..[3] Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K.
-          dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001 Jan 1;29(1):308-11.
-        ..[4] Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center
-          for Biotechnology Information, National Library of Medicine. dbSNP accession: rs3094315,
-          rs11928389, rs2500347, and rs964481 (dbSNP Build ID: 151). Available from:
-          http://www.ncbi.nlm.nih.gov/SNP/
-        """
-        assembly = None
-
-        rsids = ['rs3094315', 'rs11928389', 'rs2500347', 'rs964481']
-        df = pd.DataFrame({36: [742429, 50908372, 143649677, 27566744],
-                           37: [752566, 50927009, 144938320, 27656823],
-                           38: [817186, 50889578, 148946169, 27638706]}, index=rsids)
-
-        for rsid in rsids:
-            if rsid in snps.index:
-                assembly = self._lookup_assembly_with_snp_pos(snps.loc[rsid].pos, df.loc[rsid])
-
-            if assembly is not None:
-                break
-
-        return assembly
-
-    @staticmethod
-    def _lookup_assembly_with_snp_pos(pos, s):
-        try:
-            return s.loc[s == pos].index[0]
-        except:
-            return None
diff --git a/lineage/snps.py b/lineage/snps.py
new file mode 100644
index 0000000..8f5884b
--- /dev/null
+++ b/lineage/snps.py
@@ -0,0 +1,426 @@
+""" Utilities for reading and analyzing SNPs within the `lineage` framework. """
+
+"""
+Copyright (C) 2018 Andrew Riha
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+from itertools import groupby, count
+import gzip
+import os
+import zipfile
+
+import numpy as np
+import pandas as pd
+
+
+class SNPs(object):
+
+    def __init__(self, file):
+        self.snps, self.source = self._read_raw_data(file)
+        self.assembly = None
+        self.assembly_detected = False
+
+        if self.snps is not None:
+            self.assembly = detect_assembly(self.snps)
+
+            if self.assembly is None:
+                self.assembly = 37  # assume GRCh37 if not detected
+            else:
+                self.assembly_detected = True
+
+    @property
+    def assembly_name(self):
+        """ Get the name of the assembly of ``SNPs``.
+
+        Returns
+        -------
+        str
+        """
+        return get_assembly_name(self.assembly)
+
+    @property
+    def snp_count(self):
+        """ Count of SNPs.
+
+        Returns
+        -------
+        int
+        """
+        return get_snp_count(self.snps)
+
+    @property
+    def chromosomes(self):
+        """ Chromosomes of ``SNPs``.
+
+        Returns
+        -------
+        list
+            list of str chromosomes (e.g., ['1', '2', '3', 'MT'], empty list if no chromosomes
+        """
+        return get_chromosomes(self.snps)
+
+    @property
+    def chromosomes_summary(self):
+        """ Summary of the chromosomes of ``SNPs``.
+
+        Returns
+        -------
+        str
+            human-readable listing of chromosomes (e.g., '1-3, MT'), empty str if no chromosomes
+        """
+        return get_chromosomes_summary(self.snps)
+
+    def _read_raw_data(self, file):
+        try:
+            if not os.path.exists(file):
+                print(file + ' does not exist; skipping')
+                return None, ''
+
+            # peek into files to determine the data format
+            if '.zip' in file:
+                with zipfile.ZipFile(file) as z:
+                    with z.open(z.namelist()[0], 'r') as f:
+                        # https://stackoverflow.com/a/606199
+                        line = f.readline().decode('utf-8')
+            elif '.gz' in file:
+                with gzip.open(file, 'rt') as f:
+                    line = f.readline()
+            else:
+                with open(file, 'r') as f:
+                    line = f.readline()
+
+            if '23andMe' in line:
+                return self._read_23andme(file)
+            elif 'Ancestry' in line:
+                return self._read_ancestry(file)
+            elif line[:4] == 'RSID':
+                return self._read_ftdna(file)
+            elif 'lineage' in line:
+                return self._read_lineage_csv(file)
+            elif line[:4] == 'rsid':
+                return self._read_generic_csv(file)
+            else:
+                return None, ''
+        except Exception as err:
+            print(err)
+            return None, ''
+
+    @staticmethod
+    def _read_23andme(file):
+        """ Read and parse 23andMe file.
+
+        https://www.23andme.com
+
+        Parameters
+        ----------
+        file : str
+            path to file
+
+        Returns
+        -------
+        pandas.DataFrame
+            individual's genetic data normalized for use with `lineage`
+        str
+            name of data source
+        """
+        df = pd.read_csv(file, comment='#', sep='\t', na_values='--',
+                               names=['rsid', 'chrom', 'pos', 'genotype'],
+                               index_col=0, dtype={'chrom': object})
+
+        return df, '23andMe'
+
+    @staticmethod
+    def _read_ftdna(file):
+        """ Read and parse Family Tree DNA (FTDNA) file.
+
+        https://www.familytreedna.com
+
+        Parameters
+        ----------
+        file : str
+            path to file
+
+        Returns
+        -------
+        pandas.DataFrame
+            individual's genetic data normalized for use with `lineage`
+        str
+            name of data source
+        """
+        df = pd.read_csv(file, skiprows=1, na_values='--',
+                         names=['rsid', 'chrom', 'pos', 'genotype'],
+                         index_col=0, dtype={'chrom': object})
+
+        # remove incongruous data
+        df = df.drop(df.loc[df['chrom'] == '0'].index)
+        df = df.drop(df.loc[df.index == 'RSID'].index)  # second header for concatenated data
+
+        # if second header existed, pos dtype will be object (should be np.int64)
+        df['pos'] = df['pos'].astype(np.int64)
+
+        return df, 'FTDNA'
+
+    @staticmethod
+    def _read_ancestry(file):
+        """ Read and parse Ancestry.com file.
+
+        http://www.ancestry.com
+
+        Parameters
+        ----------
+        file : str
+            path to file
+
+        Returns
+        -------
+        pandas.DataFrame
+            individual's genetic data normalized for use with `lineage`
+        str
+            name of data source
+        """
+        df = pd.read_csv(file, comment='#', header=0, sep='\t', na_values=0,
+                         names=['rsid', 'chrom', 'pos', 'allele1', 'allele2'],
+                         index_col=0, dtype={'chrom': object})
+
+        # create genotype column from allele columns
+        df['genotype'] = df['allele1'] + df['allele2']
+
+        # delete allele columns
+        # http://stackoverflow.com/a/13485766
+        del df['allele1']
+        del df['allele2']
+
+        # https://redd.it/5y90un
+        df.ix[np.where(df['chrom'] == '23')[0], 'chrom'] = 'X'
+        df.ix[np.where(df['chrom'] == '24')[0], 'chrom'] = 'Y'
+        df.ix[np.where(df['chrom'] == '25')[0], 'chrom'] = 'PAR'
+        df.ix[np.where(df['chrom'] == '26')[0], 'chrom'] = 'MT'
+
+        return df, 'AncestryDNA'
+
+    @staticmethod
+    def _read_lineage_csv(file):
+        """ Read and parse CSV file generated by lineage.
+
+        Parameters
+        ----------
+        file : str
+            path to file
+
+        Returns
+        -------
+        pandas.DataFrame
+            individual's genetic data normalized for use with `lineage`
+        str
+            name of data source
+        """
+        df = pd.read_csv(file, comment='#', header=0, na_values='--',
+                         names=['rsid', 'chrom', 'pos', 'genotype'],
+                         index_col=0, dtype={'chrom': object, 'pos': np.int64})
+
+        return df, 'lineage'
+
+    @staticmethod
+    def _read_generic_csv(file):
+        """ Read and parse generic CSV file.
+
+        Notes
+        -----
+        Assumes columns are 'rsid', 'chrom' / 'chromosome', 'pos' / 'position', and 'genotype';
+        values are comma separated; unreported genotypes are indicated by '--'; and one header row
+        precedes data. For example:
+
+            rsid,chromosome,position,genotype
+            rs1,1,1,AA
+            rs2,1,2,CC
+            rs3,1,3,--
+
+        Parameters
+        ----------
+        file : str
+            path to file
+
+        Returns
+        -------
+        pandas.DataFrame
+            individual's genetic data normalized for use with `lineage`
+        str
+            name of data source
+        """
+        df = pd.read_csv(file, skiprows=1, na_values='--',
+                         names=['rsid', 'chrom', 'pos', 'genotype'],
+                         index_col=0, dtype={'chrom': object, 'pos': np.int64})
+
+        return df, 'generic'
+
+
+def detect_assembly(snps):
+    """ Detect assembly of SNPs.
+
+    Use the coordinates of common SNPs to identify the assembly / build of a genotype file
+    that is being loaded.
+
+    Notes
+    -----
+    rs3094315 : plus strand in 36, 37, and 38
+    rs11928389 : plus strand in 36, minus strand in 37 and 38
+    rs2500347 : plus strand in 36 and 37, minus strand in 38
+    rs964481 : plus strand in 36, 37, and 38
+
+    Parameters
+    ----------
+    snps : pandas.DataFrame
+        SNPs to add
+
+    Returns
+    -------
+    int
+        detected assembly of SNPs, else None
+
+    References
+    ----------
+    ..[1] Yates et. al. (doi:10.1093/bioinformatics/btu613),
+      http://europepmc.org/search/?query=DOI:10.1093/bioinformatics/btu613
+    ..[2] Zerbino et. al. (doi.org/10.1093/nar/gkx1098), https://doi.org/10.1093/nar/gkx1098
+    ..[3] Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K.
+      dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001 Jan 1;29(1):308-11.
+    ..[4] Database of Single Nucleotide Polymorphisms (dbSNP). Bethesda (MD): National Center
+      for Biotechnology Information, National Library of Medicine. dbSNP accession: rs3094315,
+      rs11928389, rs2500347, and rs964481 (dbSNP Build ID: 151). Available from:
+      http://www.ncbi.nlm.nih.gov/SNP/
+    """
+    def lookup_assembly_with_snp_pos(pos, s):
+        try:
+            return s.loc[s == pos].index[0]
+        except:
+            return None
+
+    assembly = None
+
+    rsids = ['rs3094315', 'rs11928389', 'rs2500347', 'rs964481']
+    df = pd.DataFrame({36: [742429, 50908372, 143649677, 27566744],
+                       37: [752566, 50927009, 144938320, 27656823],
+                       38: [817186, 50889578, 148946169, 27638706]}, index=rsids)
+
+    for rsid in rsids:
+        if rsid in snps.index:
+            assembly = lookup_assembly_with_snp_pos(snps.loc[rsid].pos, df.loc[rsid])
+
+        if assembly is not None:
+            break
+
+    return assembly
+
+
+def get_assembly_name(assembly):
+    """ Get the name of an assembly.
+
+    Parameters
+    ----------
+    assembly : int {36, 37, 38}
+
+    Returns
+    -------
+    str
+        empty str if `assembly` is None
+    """
+
+    if assembly is None:
+        return ''
+    elif assembly == 36:
+        return 'NCBI36'
+    elif assembly == 38:
+        return 'GRCh38'
+    else:
+        return 'GRCh37'
+
+
+def get_snp_count(snps):
+    """ Count of SNPs.
+
+    Parameters
+    ----------
+    snps : pandas.DataFrame
+
+    Returns
+    -------
+    int
+    """
+
+    if snps is not None:
+        return len(snps)
+    else:
+        return 0
+
+
+def get_chromosomes(snps):
+    """ Get the chromosomes of SNPs.
+
+    Parameters
+    ----------
+    snps : pandas.DataFrame
+
+    Returns
+    -------
+    list
+        list of str chromosomes (e.g., ['1', '2', '3', 'MT'], empty list if no chromosomes
+    """
+
+    if snps is not None:
+        return list(pd.unique(snps['chrom']))
+    else:
+        return []
+
+
+def get_chromosomes_summary(snps):
+    """ Summary of the chromosomes of SNPs.
+
+    Parameters
+    ----------
+    snps : pandas.DataFrame
+
+    Returns
+    -------
+    str
+        human-readable listing of chromosomes (e.g., '1-3, MT'), empty str if no chromosomes
+    """
+
+    if snps is not None:
+        chroms = list(pd.unique(snps['chrom']))
+
+        int_chroms = [int(chrom) for chrom in chroms if chrom.isdigit()]
+        str_chroms = [chrom for chrom in chroms if not chrom.isdigit()]
+
+        # https://codereview.stackexchange.com/a/5202
+        def as_range(iterable):
+            l = list(iterable)
+            if len(l) > 1:
+                return '{0}-{1}'.format(l[0], l[-1])
+            else:
+                return '{0}'.format(l[0])
+
+        # create str representations
+        int_chroms = ', '.join(as_range(g) for _, g in
+                               groupby(int_chroms, key=lambda n, c=count(): n - next(c)))
+        str_chroms = ', '.join(str_chroms)
+
+        if str_chroms != '':
+            int_chroms += ', '
+
+        return int_chroms + str_chroms
+    else:
+        return ''


Test patch:
diff --git a/tests/input/chromosomes.csv b/tests/input/chromosomes.csv
new file mode 100644
index 0000000..957d35f
--- /dev/null
+++ b/tests/input/chromosomes.csv
@@ -0,0 +1,6 @@
+rsid,chromosome,position,genotype
+rs1,1,1,AA
+rs2,2,1,CC
+rs3,3,1,GG
+rs4,5,1,TT
+rs5,MT,1,A
diff --git a/tests/test_individual.py b/tests/test_individual.py
index 1be1772..872e814 100644
--- a/tests/test_individual.py
+++ b/tests/test_individual.py
@@ -17,6 +17,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 import gzip
+import os
 import shutil
 import zipfile
 
@@ -67,6 +68,7 @@ def test_name(l):
 def test_snps_23andme(l, generic_snps):
     # https://www.23andme.com
     ind = l.create_individual('', 'tests/input/23andme.txt')
+    assert ind.source == '23andMe'
     pd.testing.assert_frame_equal(ind.snps, generic_snps)
 
 
@@ -75,12 +77,14 @@ def test_snps_23andme_zip(l, generic_snps):
         # https://stackoverflow.com/a/16104667
         f.write('tests/input/23andme.txt', arcname='23andme.txt')
     ind = l.create_individual('', 'tests/input/23andme.txt.zip')
+    assert ind.source == '23andMe'
     pd.testing.assert_frame_equal(ind.snps, generic_snps)
 
 
 def test_snps_ftdna(l, generic_snps):
     # https://www.familytreedna.com
     ind = l.create_individual('', 'tests/input/ftdna.csv')
+    assert ind.source == 'FTDNA'
     pd.testing.assert_frame_equal(ind.snps, generic_snps)
 
 
@@ -89,29 +93,77 @@ def test_snps_ftdna_gzip(l, generic_snps):
         with gzip.open('tests/input/ftdna.csv.gz', 'wb') as f_out:
             shutil.copyfileobj(f_in, f_out)
     ind = l.create_individual('', 'tests/input/ftdna.csv.gz')
+    assert ind.source == 'FTDNA'
     pd.testing.assert_frame_equal(ind.snps, generic_snps)
 
 
 def test_snps_ancestry(l, generic_snps):
     # https://www.ancestry.com
     ind = l.create_individual('', 'tests/input/ancestry.txt')
+    assert ind.source == 'AncestryDNA'
     pd.testing.assert_frame_equal(ind.snps, generic_snps)
 
 
+def test_source_lineage(l):
+    ind = l.create_individual('', 'tests/input/chromosomes.csv')
+    assert ind.source == 'generic'
+    file = ind.save_snps()
+    ind_saved_snps = l.create_individual('', file)
+    assert ind_saved_snps.source == 'lineage'
+    pd.testing.assert_frame_equal(ind.snps, ind_saved_snps.snps)
+
+
+def test_source_generic(l):
+    ind = l.create_individual('', 'tests/input/NCBI36.csv')
+    assert ind.source == 'generic'
+
+
 def test_snps_None(l):
     ind = l.create_individual('')
     assert ind.snps is None
 
 
+def test_snp_count(l):
+    ind = l.create_individual('', 'tests/input/NCBI36.csv')
+    assert ind.snp_count == 4
+
+
+def test_snp_count_None(l):
+    ind = l.create_individual('')
+    assert ind.snp_count == 0
+
+
+def test_chromosomes(l):
+    ind = l.create_individual('', 'tests/input/chromosomes.csv')
+    assert ind.chromosomes == ['1', '2', '3', '5', 'MT']
+
+
+def test_chromosomes_None(l):
+    ind = l.create_individual('')
+    assert ind.chromosomes == []
+
+
+def test_chromosomes_summary(l):
+    ind = l.create_individual('', 'tests/input/chromosomes.csv')
+    assert ind.chromosomes_summary == '1-3, 5, MT'
+
+
+def test_chromosomes_summary_None(l):
+    ind = l.create_individual('')
+    assert ind.chromosomes_summary == ''
+
+
 def test_assembly(l):
     ind = l.create_individual('', 'tests/input/NCBI36.csv')
     assert ind.assembly == 36
+    assert ind.assembly_name == 'NCBI36'
 
 
 def test_load_snps_list(l, snps_GRCh37):
     ind = l.create_individual('')
     ind.load_snps(['tests/input/GRCh37.csv', 'tests/input/GRCh37.csv'])
     pd.testing.assert_frame_equal(ind.snps, snps_GRCh37)
+    assert ind.source == 'generic, generic'
 
 
 def test_load_snps_None(l):
@@ -156,7 +208,7 @@ def test_load_snps_assembly_mismatch_exceed_discrepant_genotypes_threshold(l):
 
 def test_save_snps(l, snps_GRCh37):
     ind = l.create_individual('test save snps', 'tests/input/GRCh37.csv')
-    assert ind.save_snps()
+    assert os.path.relpath(ind.save_snps()) == 'output/test_save_snps.csv'
     ind_saved_snps = l.create_individual('', 'output/test_save_snps.csv')
     pd.testing.assert_frame_equal(ind_saved_snps.snps, snps_GRCh37)
 
@@ -187,6 +239,7 @@ def test_remap_snps_36_to_37(l, snps_GRCh37):
     ind = l.create_individual('', 'tests/input/NCBI36.csv')
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(37)
     assert ind.assembly == 37  # TODO: handle partial remapping; see #24
+    assert ind.assembly_name == 'GRCh37'
     if len(chromosomes_remapped) == 2:
         assert len(chromosomes_not_remapped) == 0
         pd.testing.assert_frame_equal(ind.snps, snps_GRCh37)
@@ -196,6 +249,7 @@ def test_remap_snps_37_to_36(l, snps_NCBI36):
     ind = l.create_individual('', 'tests/input/GRCh37.csv')
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(36)
     assert ind.assembly == 36  # TODO: handle partial remapping; see #24
+    assert ind.assembly_name == 'NCBI36'
     if len(chromosomes_remapped) == 2:
         assert len(chromosomes_not_remapped) == 0
         pd.testing.assert_frame_equal(ind.snps, snps_NCBI36)
@@ -205,6 +259,7 @@ def test_remap_snps_37_to_38(l, snps_GRCh38):
     ind = l.create_individual('', 'tests/input/GRCh37.csv')
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(38)
     assert ind.assembly == 38  # TODO: handle partial remapping; see #24
+    assert ind.assembly_name == 'GRCh38'
     if len(chromosomes_remapped) == 2:
         assert len(chromosomes_not_remapped) == 0
         pd.testing.assert_frame_equal(ind.snps, snps_GRCh38)
@@ -214,6 +269,7 @@ def test_remap_snps_37_to_37(l, snps_GRCh37):
     ind = l.create_individual('', 'tests/input/GRCh37.csv')
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(37)
     assert ind.assembly == 37
+    assert ind.assembly_name == 'GRCh37'
     assert len(chromosomes_remapped) == 0
     assert len(chromosomes_not_remapped) == 2
     pd.testing.assert_frame_equal(ind.snps, snps_GRCh37)
@@ -224,6 +280,7 @@ def test_remap_snps_no_EnsemblRestClient(l):
     ind._ensembl_rest_client = None
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(38)
     assert ind.assembly == 37
+    assert ind.assembly_name == 'GRCh37'
     assert len(chromosomes_remapped) == 0
     assert len(chromosomes_not_remapped) == 2
 
@@ -240,35 +297,11 @@ def test_remap_snps_invalid_assembly(l):
     ind = l.create_individual('', 'tests/input/GRCh37.csv')
     chromosomes_remapped, chromosomes_not_remapped = ind.remap_snps(-1)
     assert ind.assembly == 37
+    assert ind.assembly_name == 'GRCh37'
     assert len(chromosomes_remapped) == 0
     assert len(chromosomes_not_remapped) == 2
 
 
-def test__read_23andme_None(l):
-    ind = l.create_individual('')
-    assert ind._read_23andme(None) is None
-
-
-def test__read_ftdna_None(l):
-    ind = l.create_individual('')
-    assert ind._read_ftdna(None) is None
-
-
-def test__read_ancestry_None(l):
-    ind = l.create_individual('')
-    assert ind._read_ancestry(None) is None
-
-
-def test__read_generic_csv_None(l):
-    ind = l.create_individual('')
-    assert ind._read_generic_csv(None) is None
-
-
-def test__lookup_assembly_with_snp_pos_None(l):
-    ind = l.create_individual('')
-    assert ind._lookup_assembly_with_snp_pos(None, None) is None
-
-
 def test___repr__(l):
     ind = l.create_individual('test')
     assert "Individual('test')" == ind.__repr__()
diff --git a/tests/test_snps.py b/tests/test_snps.py
new file mode 100644
index 0000000..df4e070
--- /dev/null
+++ b/tests/test_snps.py
@@ -0,0 +1,81 @@
+"""
+Copyright (C) 2018 Andrew Riha
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+import pandas as pd
+import pytest
+
+
+def create_snp_df(rsid, chrom, pos, genotype):
+    df = pd.DataFrame({'rsid': rsid, 'chrom': chrom, 'pos': pos, 'genotype': genotype},
+                      columns=['rsid', 'chrom', 'pos', 'genotype'])
+    df = df.set_index('rsid')
+    return df
+
+
+@pytest.fixture(scope='module')
+def snps_discrepant_pos():
+    return create_snp_df(rsid=['rs3094315'], chrom=['1'], pos=[1], genotype=['AA'])
+
+
+@pytest.fixture(scope='module')
+def snps_GRCh38():
+    from lineage.snps import SNPs
+    return SNPs('tests/input/GRCh38.csv')
+
+
+@pytest.fixture(scope='module')
+def snps():
+    from lineage.snps import SNPs
+    return SNPs('tests/input/chromosomes.csv')
+
+
+@pytest.fixture(scope='module')
+def snps_none():
+    from lineage.snps import SNPs
+    return SNPs(None)
+
+
+def test_assembly_name(snps_GRCh38):
+    assert snps_GRCh38.assembly_name == 'GRCh38'
+
+
+def test_snp_count(snps):
+    assert snps.snp_count == 5
+
+
+def test_chromosomes(snps):
+    assert snps.chromosomes == ['1', '2', '3', '5', 'MT']
+
+
+def test_chromosomes_summary(snps):
+    assert snps.chromosomes_summary == '1-3, 5, MT'
+
+
+def test__read_raw_data(snps_none):
+    assert snps_none.snps is None
+    assert snps_none.source == ''
+
+
+def test__lookup_assembly_with_snp_pos_None(snps_discrepant_pos):
+    from lineage.snps import detect_assembly
+    assert detect_assembly(snps_discrepant_pos) is None
+
+
+def test_get_assembly_name_None():
+    from lineage.snps import get_assembly_name
+    assert get_assembly_name(None) is ''