Repository: jjmccollum/teiphy

Issue:
Should teiphy convert all variation units, including trivial ones?
Since we decided in #56 that ascertainment bias correction should be handled by the user's phylogenetic software rather by `teiphy`, it is worth asking if `teiphy` should convert _all_ variation units, and not just those with more than one substantive variant reading (after collapsing readings of ignored types with their "parent" readings), to sites in the output file, in case these are needed for ascertainment bias correction.

If we decide this is necessary, then this change should be tested in a dedicated branch, as whether it is worth implementing will depend on whether it causes any of the software tested in the GitHub workflows to throw any errors. (I know that some programs will at least issue warnings about constant sites.)

Discussion:

Patch:
diff --git a/.github/workflows/iqtree.yml b/.github/workflows/iqtree.yml
index 9016ce9..a4048d2 100644
--- a/.github/workflows/iqtree.yml
+++ b/.github/workflows/iqtree.yml
@@ -26,5 +26,5 @@ jobs:
         sudo apt-get install -y iqtree
     - name: Phylogenetics Run
       run: |
-        iqtree2 -s iqtree_example.nexus -m MK+ASC -bb 1000
+        iqtree2 -s iqtree_example.nexus -m MK -bb 1000
       
diff --git a/docs/advanced.rst b/docs/advanced.rst
index d6482f7..32e87ca 100644
--- a/docs/advanced.rst
+++ b/docs/advanced.rst
@@ -342,6 +342,23 @@ attestations, then you can do so via
 
 assuming that there are at most three multiple attestations in any unit.
 
+Ascertainment Bias
+------------------
+
+To facilitate accurate branch length estimation, ``teiphy`` includes all constant sites (i.e., variation units with only one substantive reading after trivial readings have been merged with their parents) in its output by default.
+(The only exception is output for the STEMMA program, which is parsimony-based and does not estimate branch lengths.)
+If you want only informative sites to be included in your output, then you can ensure that this happens by including the ``--drop-constant`` flag in your command.
+
+Note that if you do exclude constant sites, or if your TEI XML collation only includes locations in the textual tradition where the readings differ, 
+then it will be important to correct for ascertainment bias when performing phylogenetic analysis which utilizes the branch lengths 
+(such as maximum-likelihood or Bayesian techniques).
+This is because when the data is filtered to only include locations where changes have occurred, the analysis can be biased towards longer branch lengths.
+There are methods to correct for this, for example, as discussed in Lewis, 2001 :cite:p:`lewis_likelihood_2001`.
+Because this correction is amply supported by existing phylogenetic software packages, ``teiphy`` will not perform this correction itself.
+You should consult the documentation for your phylogenetic software package to determine how it applies ascertainment bias correction.
+For example, `IQ-TREE <http://www.iqtree.org/doc/Substitution-Models#ascertainment-bias-correction>`_ allows users to correct for ascertainment bias by adding ``+ASC`` to the model name.
+(Just keep in mind that your input must be free of constant sites for this to work!)
+
 Supported Output Formats and Options
 ------------------------------------
 
@@ -413,15 +430,4 @@ To run this script with the example input in verbose mode with the settings desc
 
    teiphy -t reconstructed -t defective -t orthographic -t subreading -m lac -m overlap -s"*" -s T -s /1 -s /2 -s /3 --fill-correctors --verbose example\ubs_ephesians.xml ubs_ephesians.nxs
 
-from the command line.
-
-Ascertainment Bias
-------------------
-
-If the apparatus in TEI only includes locations in the textual tradition where the readings differ, 
-then it will be important to correct for ascertainment bias when performing phylogenetic analysis which utilizes the branch lengths 
-(such as when using maximum likelihood or Bayesian techniques). 
-This is because the data is filtered to only include locations where changes have occurred and this leads the analysis to have a bias towards longer branch lengths.
-There are methods to correct for this, for example, as discussed in Lewis, 2001 :cite:p:`lewis_likelihood_2001`. 
-``teiphy`` will not perform this correction itself and you should consult the documentation for the phylogenetic software package being used to 
-determine how. For example, `IQ-TREE <http://www.iqtree.org/doc/Substitution-Models#ascertainment-bias-correction>`_ allows users to correct for ascertainment bias by adding ``+ASC`` to the model name.
\ No newline at end of file
+from the command line.
\ No newline at end of file
diff --git a/example/ubs_ephesians.xml b/example/ubs_ephesians.xml
index 173bab8..92b95d7 100644
--- a/example/ubs_ephesians.xml
+++ b/example/ubs_ephesians.xml
@@ -606,7 +606,6 @@
                     </listRelation>
                 </note>
             </app>
-            <!-- this unit should be ignored -->
             <app xml:id="B10K4V28U24-26">
                 <lem><w>το</w><w>αγαθον</w></lem>
                 <rdg n="1" wit="UBS P46 01 02 03 06 010 012 18 33 35 88 424 606 915 1175 1505 1611 1739 1881 1910 2495 vg syrp syrh copbo Ambrosiaster Chrysostom Jerome MariusVictorinus Pelagius TheodoreOfMopsuestia"><w>το</w><w>αγαθον</w></rdg>
@@ -921,7 +920,6 @@
                     </listRelation>
                 </note>
             </app>
-            <!-- this unit should always be ignored -->
             <app xml:id="B10K6V20U12">
                 <lem><w>ινα</w></lem>
                 <rdg n="1" wit="UBS P46 01 02 03 06 010 012 18 33 35 88 424 606 915 1175 1505 1611 1739 1881 2495 vg syrp syrh copsa copbo Ambrosiaster Chrysostom Jerome MariusVictorinus Pelagius TheodoreOfMopsuestia"><w>ινα</w></rdg>
diff --git a/teiphy/collation.py b/teiphy/collation.py
index 5e7fef0..47baddd 100644
--- a/teiphy/collation.py
+++ b/teiphy/collation.py
@@ -66,8 +66,9 @@ class Collation:
         self.witness_index_by_id = {}
         self.variation_units = []
         self.readings_by_witness = {}
-        self.substantive_variation_unit_ids = []
+        self.variation_unit_ids = []
         self.substantive_variation_unit_reading_tuples = []
+        self.substantive_readings_by_variation_unit_id = {}
         # Now parse the XML tree to populate these data structures:
         if self.verbose:
             print("Initializing collation...")
@@ -108,7 +109,7 @@ class Collation:
 
     def parse_list_wit(self, xml: et.ElementTree):
         """Given an XML tree for a collation, populates its list of witnesses from its listWit element.
-        If the XML tree does not contain a listWit element, then an Exception is thrown listing all distinct witness sigla encountered in the collation.
+        If the XML tree does not contain a listWit element, then a ParsingException is thrown listing all distinct witness sigla encountered in the collation.
 
         Args:
             xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
@@ -213,38 +214,35 @@ class Collation:
         """
         # In a first pass, populate lists of substantive (variation unit ID, reading ID) tuples and reading labels
         # and a map from reading IDs to the indices of their parent substantive reading in this unit:
-        substantive_variation_unit_reading_tuples = []
         reading_id_to_index = {}
+        self.substantive_readings_by_variation_unit_id[vu.id] = []
         for rdg in vu.readings:
             # If this reading is missing (e.g., lacunose or inapplicable due to an overlapping variant) or targets another reading, then skip it:
             if rdg.type in self.missing_reading_types or len(rdg.certainties) > 0:
                 continue
             # If this reading is trivial, then map it to the last substantive index:
             if rdg.type in self.trivial_reading_types:
-                reading_id_to_index[rdg.id] = len(substantive_variation_unit_reading_tuples) - 1
+                reading_id_to_index[rdg.id] = len(self.substantive_readings_by_variation_unit_id[vu.id]) - 1
                 continue
             # Otherwise, the reading is substantive: add it to the map and update the last substantive index:
-            substantive_variation_unit_reading_tuples.append(tuple([vu.id, rdg.id]))
-            reading_id_to_index[rdg.id] = len(substantive_variation_unit_reading_tuples) - 1
+            self.substantive_readings_by_variation_unit_id[vu.id].append(rdg.id)
+            self.substantive_variation_unit_reading_tuples.append(tuple([vu.id, rdg.id]))
+            reading_id_to_index[rdg.id] = len(self.substantive_readings_by_variation_unit_id[vu.id]) - 1
         # If the list of substantive readings only contains one entry, then this variation unit is not informative;
         # return an empty dictionary and add nothing to the list of substantive reading labels:
         if self.verbose:
             print(
                 "Variation unit %s has %d substantive readings."
-                % (vu.id, len(substantive_variation_unit_reading_tuples))
+                % (vu.id, len(self.substantive_readings_by_variation_unit_id[vu.id]))
             )
         readings_by_witness_for_unit = {}
-        if len(substantive_variation_unit_reading_tuples) <= 1:
-            return readings_by_witness_for_unit
-        # Otherwise, add these substantive reading ID tuples and labels to their corresponding list and initialize the output dictionary with empty sets for all base witnesses:
-        for substantive_variation_unit_reading_tuple in substantive_variation_unit_reading_tuples:
-            self.substantive_variation_unit_reading_tuples.append(substantive_variation_unit_reading_tuple)
+        # Initialize the output dictionary with empty sets for all base witnesses:
         for wit in self.witnesses:
-            readings_by_witness_for_unit[wit.id] = [0] * len(substantive_variation_unit_reading_tuples)
+            readings_by_witness_for_unit[wit.id] = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
         # In a second pass, assign each base witness a set containing the readings it supports in this unit:
         for rdg in vu.readings:
             # Initialize the dictionary indicating support for this reading (or its disambiguations):
-            rdg_support = [0] * len(substantive_variation_unit_reading_tuples)
+            rdg_support = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
             # If this is a missing reading (e.g., a lacuna or an overlap), then we can skip this reading, as its corresponding set will be empty:
             if rdg.type in self.missing_reading_types:
                 continue
@@ -297,14 +295,14 @@ class Collation:
         t0 = time.time()
         # Initialize the data structures to be populated here:
         self.readings_by_witness = {}
-        self.substantive_variation_unit_ids = []
+        self.variation_unit_ids = []
         for wit in self.witnesses:
             self.readings_by_witness[wit.id] = []
         # Populate them for each variation unit:
         for vu in self.variation_units:
             readings_by_witness_for_unit = self.get_readings_by_witness_for_unit(vu)
             if len(readings_by_witness_for_unit) > 0:
-                self.substantive_variation_unit_ids.append(vu.id)
+                self.variation_unit_ids.append(vu.id)
             for wit in readings_by_witness_for_unit:
                 self.readings_by_witness[wit].append(readings_by_witness_for_unit[wit])
         # Optionally, fill the lacunae of the correctors:
@@ -326,7 +324,7 @@ class Collation:
         if self.verbose:
             print(
                 "Populated dictionary for %d witnesses over %d substantive variation units in %0.4fs."
-                % (len(self.witnesses), len(self.substantive_variation_unit_ids), t1 - t0)
+                % (len(self.witnesses), len(self.variation_unit_ids), t1 - t0)
             )
         return
 
@@ -355,6 +353,7 @@ class Collation:
     def to_nexus(
         self,
         file_addr: Union[Path, str],
+        drop_constant: bool = False,
         char_state_labels: bool = True,
         frequency: bool = False,
         ambiguous_as_missing: bool = False,
@@ -365,6 +364,7 @@ class Collation:
 
         Args:
             file_addr: A string representing the path to an output NEXUS file; the file type should be .nex, .nexus, or .nxs.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             char_state_labels: An optional flag indicating whether or not to include the CharStateLabels block.
             frequency: An optional flag indicating whether to use the StatesFormat=Frequency setting
                 instead of the StatesFormat=StatesPresent setting
@@ -378,16 +378,24 @@ class Collation:
             mrbayes: An optional flag indicating whether to add a MrBayes block that specifies model settings and age calibrations for witnesses.
                 This option is intended for inputs to MrBayes.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Start by calculating the values we will be using here:
         ntax = len(self.witnesses)
-        nchar = (
-            len(self.readings_by_witness[self.witnesses[0].id]) if ntax > 0 else 0
-        )  # if the number of taxa is 0, then the number of characters is irrelevant
+        nchar = len(substantive_variation_unit_ids)
         taxlabels = [slugify(wit.id, lowercase=False, separator='_') for wit in self.witnesses]
         max_taxlabel_length = max(
             [len(taxlabel) for taxlabel in taxlabels]
         )  # keep track of the longest taxon label for tabular alignment purposes
-        charlabels = [slugify(vu_id, lowercase=False, separator='_') for vu_id in self.substantive_variation_unit_ids]
+        charlabels = [slugify(vu_id, lowercase=False, separator='_') for vu_id in substantive_variation_unit_ids]
         missing_symbol = '?'
         symbols = self.get_nexus_symbols()
         with open(file_addr, "w", encoding="utf-8") as f:
@@ -407,8 +415,6 @@ class Collation:
             # If the char_state_labels is set, then write the labels for character-state labels, with each on its own line:
             if char_state_labels:
                 f.write("\tCharStateLabels")
-                substantive_variation_unit_ids_set = set(self.substantive_variation_unit_ids)
-                substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
                 vu_ind = 1
                 for vu in self.variation_units:
                     if vu.id not in substantive_variation_unit_ids_set:
@@ -438,7 +444,10 @@ class Collation:
                 taxlabel = taxlabels[i]
                 if frequency:
                     sequence = "\n\t\t" + taxlabel
-                    for rdg_support in self.readings_by_witness[wit.id]:
+                    for j, vu_id in enumerate(self.variation_unit_ids):
+                        if vu_id not in substantive_variation_unit_ids_set:
+                            continue
+                        rdg_support = self.readings_by_witness[wit.id][j]
                         sequence += "\n\t\t\t"
                         # If this reading is lacunose in this witness, then use the missing character:
                         if sum(rdg_support) == 0:
@@ -455,13 +464,16 @@ class Collation:
                     sequence = "\n\t\t" + taxlabel
                     # Add enough space after this label ensure that all sequences are nicely aligned:
                     sequence += " " * (max_taxlabel_length - len(taxlabel) + 1)
-                    for rdg_support in self.readings_by_witness[wit.id]:
+                    for j, vu_id in enumerate(self.variation_unit_ids):
+                        if vu_id not in substantive_variation_unit_ids_set:
+                            continue
+                        rdg_support = self.readings_by_witness[wit.id][j]
                         # If this reading is lacunose in this witness, then use the missing character:
                         if sum(rdg_support) == 0:
                             sequence += missing_symbol
                             continue
                         rdg_inds = [
-                            i for i, w in enumerate(rdg_support) if w > 0
+                            k for k, w in enumerate(rdg_support) if w > 0
                         ]  # the index list consists of the indices of all readings with any degree of certainty assigned to them
                         # For singleton readings, just print the symbol:
                         if len(rdg_inds) == 1:
@@ -588,18 +600,27 @@ class Collation:
         hennig86_symbols = possible_symbols[:nsymbols]
         return hennig86_symbols
 
-    def to_hennig86(self, file_addr: Union[Path, str]):
+    def to_hennig86(self, file_addr: Union[Path, str], drop_constant: bool = False):
         """Writes this Collation to a file in Hennig86 format with the given address.
         Note that because Hennig86 format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
 
         Args:
             file_addr: A string representing the path to an output file.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Start by calculating the values we will be using here:
         ntax = len(self.witnesses)
-        nchar = (
-            len(self.readings_by_witness[self.witnesses[0].id]) if ntax > 0 else 0
-        )  # if the number of taxa is 0, then the number of characters is irrelevant
+        nchar = len(substantive_variation_unit_ids)
         taxlabels = []
         for wit in self.witnesses:
             taxlabel = wit.id
@@ -626,13 +647,16 @@ class Collation:
                 taxlabel = taxlabels[i]
                 # Add enough space after this label ensure that all sequences are nicely aligned:
                 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel) + 1))
-                for rdg_support in self.readings_by_witness[wit.id]:
+                for j, vu_id in enumerate(self.variation_unit_ids):
+                    if vu_id not in substantive_variation_unit_ids_set:
+                        continue
+                    rdg_support = self.readings_by_witness[wit.id][j]
                     # If this reading is lacunose in this witness, then use the missing character:
                     if sum(rdg_support) == 0:
                         sequence += missing_symbol
                         continue
                     rdg_inds = [
-                        i for i, w in enumerate(rdg_support) if w > 0
+                        k for k, w in enumerate(rdg_support) if w > 0
                     ]  # the index list consists of the indices of all readings with any degree of certainty assigned to them
                     # For singleton readings, just print the symbol:
                     if len(rdg_inds) == 1:
@@ -666,18 +690,27 @@ class Collation:
         phylip_symbols = possible_symbols[:nsymbols]
         return phylip_symbols
 
-    def to_phylip(self, file_addr: Union[Path, str]):
+    def to_phylip(self, file_addr: Union[Path, str], drop_constant: bool = False):
         """Writes this Collation to a file in PHYLIP format with the given address.
         Note that because PHYLIP format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
 
         Args:
             file_addr: A string representing the path to an output file.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Start by calculating the values we will be using here:
         ntax = len(self.witnesses)
-        nchar = (
-            len(self.readings_by_witness[self.witnesses[0].id]) if ntax > 0 else 0
-        )  # if the number of taxa is 0, then the number of characters is irrelevant
+        nchar = len(substantive_variation_unit_ids)
         taxlabels = []
         for wit in self.witnesses:
             taxlabel = wit.id
@@ -697,13 +730,16 @@ class Collation:
                 taxlabel = taxlabels[i]
                 # Add enough space after this label ensure that all sequences are nicely aligned:
                 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel))) + "\t"
-                for rdg_support in self.readings_by_witness[wit.id]:
+                for j, vu_id in enumerate(self.variation_unit_ids):
+                    if vu_id not in substantive_variation_unit_ids_set:
+                        continue
+                    rdg_support = self.readings_by_witness[wit.id][j]
                     # If this reading is lacunose in this witness, then use the missing character:
                     if sum(rdg_support) == 0:
                         sequence += missing_symbol
                         continue
                     rdg_inds = [
-                        i for i, w in enumerate(rdg_support) if w > 0
+                        k for k, w in enumerate(rdg_support) if w > 0
                     ]  # the index list consists of the indices of all readings with any degree of certainty assigned to them
                     # For singleton readings, just print the symbol:
                     if len(rdg_inds) == 1:
@@ -736,18 +772,27 @@ class Collation:
         fasta_symbols = possible_symbols[:nsymbols]
         return fasta_symbols
 
-    def to_fasta(self, file_addr: Union[Path, str]):
+    def to_fasta(self, file_addr: Union[Path, str], drop_constant: bool = False):
         """Writes this Collation to a file in FASTA format with the given address.
         Note that because FASTA format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
 
         Args:
             file_addr: A string representing the path to an output file.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Start by calculating the values we will be using here:
         ntax = len(self.witnesses)
-        nchar = (
-            len(self.readings_by_witness[self.witnesses[0].id]) if ntax > 0 else 0
-        )  # if the number of taxa is 0, then the number of characters is irrelevant
+        nchar = len(substantive_variation_unit_ids)
         taxlabels = []
         for wit in self.witnesses:
             taxlabel = wit.id
@@ -765,13 +810,16 @@ class Collation:
                 taxlabel = taxlabels[i]
                 # Add enough space after this label ensure that all sequences are nicely aligned:
                 sequence = ">%s\n" % taxlabel
-                for rdg_support in self.readings_by_witness[wit.id]:
+                for j, vu_id in enumerate(self.variation_unit_ids):
+                    if vu_id not in substantive_variation_unit_ids_set:
+                        continue
+                    rdg_support = self.readings_by_witness[wit.id][j]
                     # If this reading is lacunose in this witness, then use the missing character:
                     if sum(rdg_support) == 0:
                         sequence += missing_symbol
                         continue
                     rdg_inds = [
-                        i for i, w in enumerate(rdg_support) if w > 0
+                        k for k, w in enumerate(rdg_support) if w > 0
                     ]  # the index list consists of the indices of all readings with any degree of certainty assigned to them
                     # For singleton readings, just print the symbol:
                     if len(rdg_inds) == 1:
@@ -782,10 +830,11 @@ class Collation:
                 f.write("%s\n" % (sequence))
         return
 
-    def to_numpy(self, split_missing: bool = True):
+    def to_numpy(self, drop_constant: bool = False, split_missing: bool = True):
         """Returns this Collation in the form of a NumPy array, along with arrays of its row and column labels.
 
         Args:
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             split_missing: An optional flag indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings; if False, then missing data is ignored (i.e., all states are 0). Default value is True.
 
         Returns:
@@ -793,10 +842,18 @@ class Collation:
             A list of substantive reading ID strings.
             A list of witness ID strings.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Initialize the output array with the appropriate dimensions:
         reading_labels = []
-        substantive_variation_unit_ids_set = set(self.substantive_variation_unit_ids)
-        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         for vu in self.variation_units:
             if vu.id not in substantive_variation_unit_ids_set:
                 continue
@@ -807,9 +864,13 @@ class Collation:
         witness_labels = [wit.id for wit in self.witnesses]
         matrix = np.zeros((len(reading_labels), len(witness_labels)), dtype=float)
         # Then populate it with the appropriate values:
-        for col_ind, wit_id in enumerate(witness_labels):
+        col_ind = 0
+        for i, wit in enumerate(self.witnesses):
             row_ind = 0
-            for rdg_support in self.readings_by_witness[wit_id]:
+            for j, vu_id in enumerate(self.variation_unit_ids):
+                if vu_id not in substantive_variation_unit_ids_set:
+                    continue
+                rdg_support = self.readings_by_witness[wit.id][j]
                 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
                 if sum(rdg_support) == 0:
                     if split_missing:
@@ -823,19 +884,31 @@ class Collation:
                     for i in range(len(rdg_support)):
                         matrix[row_ind, col_ind] = rdg_support[i]
                         row_ind += 1
+            col_ind += 1
         return matrix, reading_labels, witness_labels
 
-    def to_distance_matrix(self, proportion=False):
+    def to_distance_matrix(self, drop_constant: bool = False, proportion=False):
         """Transforms this Collation into a NumPy distance matrix between witnesses, along with an array of its labels for the witnesses.
         Distances can be computed either as counts of disagreements (the default setting), or as proportions of disagreements over all variation units where both witnesses have singleton readings.
 
         Args:
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
 
         Returns:
             A NumPy distance matrix with a row and column for each witness.
             A list of witness ID strings.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Initialize the output array with the appropriate dimensions:
         witness_labels = [wit.id for wit in self.witnesses]
         matrix = np.zeros((len(witness_labels), len(witness_labels)), dtype=float)
@@ -851,11 +924,13 @@ class Collation:
                 if i > j:
                     matrix[i, j] = matrix[j, i]
                     continue
-                for vu_ind in range(len(self.substantive_variation_unit_ids)):
-                    wit_1_rdg_support = self.readings_by_witness[wit_1][vu_ind]
-                    wit_2_rdg_support = self.readings_by_witness[wit_2][vu_ind]
-                    wit_1_rdg_inds = [k for k, w in enumerate(wit_1_rdg_support) if w > 0]
-                    wit_2_rdg_inds = [k for k, w in enumerate(wit_2_rdg_support) if w > 0]
+                for k, vu_id in enumerate(self.variation_unit_ids):
+                    if vu_id not in substantive_variation_unit_ids_set:
+                        continue
+                    wit_1_rdg_support = self.readings_by_witness[wit_1][k]
+                    wit_2_rdg_support = self.readings_by_witness[wit_2][k]
+                    wit_1_rdg_inds = [l for l, w in enumerate(wit_1_rdg_support) if w > 0]
+                    wit_2_rdg_inds = [l for l, w in enumerate(wit_2_rdg_support) if w > 0]
                     if len(wit_1_rdg_inds) != 1 or len(wit_2_rdg_inds) != 1:
                         continue
                     extant_units += 1
@@ -869,63 +944,73 @@ class Collation:
                     matrix[i, j] = disagreements
         return matrix, witness_labels
 
-    def to_long_table(self):
+    def to_long_table(self, drop_constant: bool = False):
         """Returns this Collation in the form of a long table with columns for taxa, characters, reading indices, and reading values.
         Note that this method treats ambiguous readings as missing data.
 
+        Args:
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
+
         Returns:
             A NumPy array with columns for taxa, characters, reading indices, and reading values, and rows for each combination of these values in the matrix.
             A list of column label strings.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment:
+        substantive_variation_unit_ids = self.variation_unit_ids
+        if drop_constant:
+            substantive_variation_unit_ids = [
+                vu_id
+                for vu_id in self.variation_unit_ids
+                if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+            ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # Initialize the outputs:
         column_labels = ["taxon", "character", "state", "value"]
         long_table_list = []
-        # First, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
-        # to the readings' texts:
+        # Populate a dictionary mapping (variation unit index, reading index) tuples to reading texts:
         reading_texts_by_indices = {}
-        substantive_variation_unit_ids_set = set(self.substantive_variation_unit_ids)
-        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
-        vu_ind = 0
-        for vu in self.variation_units:
+        for j, vu in enumerate(self.variation_units):
             if vu.id not in substantive_variation_unit_ids_set:
                 continue
-            rdg_ind = 0
+            k = 0
             for rdg in vu.readings:
                 key = tuple([vu.id, rdg.id])
                 if key not in substantive_variation_unit_reading_tuples_set:
                     continue
-                indices = tuple([vu_ind, rdg_ind])
+                indices = tuple([j, k])
                 reading_texts_by_indices[indices] = rdg.text
-                rdg_ind += 1
-            if rdg_ind > 0:
-                vu_ind += 1
+                k += 1
         # Then populate the output list with the appropriate values:
         witness_labels = [wit.id for wit in self.witnesses]
         missing_symbol = '?'
-        for wit_id in witness_labels:
-            for vu_ind, rdg_support in enumerate(self.readings_by_witness[wit_id]):
-                # Get the variation unit label for this unit:
-                vu_label = self.substantive_variation_unit_ids[vu_ind]
+        for i, wit in enumerate(self.witnesses):
+            row_ind = 0
+            for j, vu_id in enumerate(self.variation_unit_ids):
+                if vu_id not in substantive_variation_unit_ids_set:
+                    continue
+                rdg_support = self.readings_by_witness[wit.id][j]
                 # Populate a list of nonzero coefficients for this reading support vector:
                 rdg_inds = [k for k, w in enumerate(rdg_support) if w > 0]
                 # If this list does not consist of exactly one reading, then treat it as missing data:
                 if len(rdg_inds) != 1:
-                    long_table_list.append([wit_id, vu_label, missing_symbol, missing_symbol])
+                    long_table_list.append([wit.id, vu_id, missing_symbol, missing_symbol])
                     continue
-                rdg_ind = rdg_inds[0]
-                rdg_text = reading_texts_by_indices[(vu_ind, rdg_ind)]
+                k = rdg_inds[0]
+                rdg_text = reading_texts_by_indices[(j, k)]
                 # Replace empty reading texts with the omission placeholder:
                 if rdg_text == "":
                     rdg_text = "om."
-                long_table_list.append([wit_id, vu_label, rdg_ind, rdg_text])
+                long_table_list.append([wit.id, vu_id, k, rdg_text])
         # Then convert the long table entries list to a NumPy array:
         long_table = np.array(long_table_list)
         return long_table, column_labels
 
-    def to_dataframe(self, long_table: bool = False, split_missing: bool = True):
+    def to_dataframe(self, drop_constant: bool = False, long_table: bool = False, split_missing: bool = True):
         """Returns this Collation in the form of a Pandas DataFrame array, including the appropriate row and column labels.
 
         Args:
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             long_table: An optional flag indicating whether or not to generate a long table with columns for taxa, characters, reading indices, and reading values.
             Note that if this option is set, ambiguous readings will be treated as missing data, and the split_missing option will be ignored.
             split_missing: An optional flag indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings; if False, then missing data is ignored (i.e., all states are 0). Default value is True.
@@ -937,46 +1022,63 @@ class Collation:
         # Proceed based on whether the long_table option is set:
         if long_table:
             # Convert the collation to a long table and get its column labels first:
-            long_table, column_labels = self.to_long_table()
+            long_table, column_labels = self.to_long_table(drop_constant=drop_constant)
             df = pd.DataFrame(long_table, columns=column_labels)
         else:
             # Convert the collation to a NumPy array and get its row and column labels first:
-            matrix, reading_labels, witness_labels = self.to_numpy(split_missing)
+            matrix, reading_labels, witness_labels = self.to_numpy(
+                drop_constant=drop_constant, split_missing=split_missing
+            )
             df = pd.DataFrame(matrix, index=reading_labels, columns=witness_labels)
         return df
 
-    def to_csv(self, file_addr: Union[Path, str], long_table: bool = False, split_missing: bool = True, **kwargs):
+    def to_csv(
+        self,
+        file_addr: Union[Path, str],
+        drop_constant: bool = False,
+        long_table: bool = False,
+        split_missing: bool = True,
+        **kwargs
+    ):
         """Writes this Collation to a comma-separated value (CSV) file with the given address.
 
         If your witness IDs are numeric (e.g., Gregory-Aland numbers), then they will be written in full to the CSV file, but Excel will likely interpret them as numbers and truncate any leading zeroes!
 
         Args:
             file_addr: A string representing the path to an output CSV file; the file type should be .csv.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             long_table: An optional flag indicating whether or not to generate a long table with columns for taxa, characters, reading indices, and reading values.
             Note that if this option is set, ambiguous readings will be treated as missing data, and the split_missing option will be ignored.
             split_missing: An optional flag indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings; if False, then missing data is ignored (i.e., all states are 0). Default value is True.
             **kwargs: Keyword arguments for pandas.DataFrame.to_csv.
         """
         # Convert the collation to a Pandas DataFrame first:
-        df = self.to_dataframe(long_table, split_missing)
+        df = self.to_dataframe(drop_constant=drop_constant, long_table=long_table, split_missing=split_missing)
         # If this is a long table, then do not include row indices:
         if long_table:
             return df.to_csv(file_addr, encoding="utf-8", index=False, **kwargs)
         return df.to_csv(file_addr, encoding="utf-8", **kwargs)
 
-    def to_excel(self, file_addr: Union[Path, str], long_table: bool = False, split_missing: bool = True):
+    def to_excel(
+        self,
+        file_addr: Union[Path, str],
+        drop_constant: bool = False,
+        long_table: bool = False,
+        split_missing: bool = True,
+    ):
         """Writes this Collation to an Excel (.xlsx) file with the given address.
 
         Since Pandas is deprecating its support for xlwt, specifying an output in old Excel (.xls) output is not recommended.
 
         Args:
             file_addr: A string representing the path to an output Excel file; the file type should be .xlsx.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             long_table: An optional flag indicating whether or not to generate a long table with columns for taxa, characters, reading indices, and reading values.
             Note that if this option is set, ambiguous readings will be treated as missing data, and the split_missing option will be ignored.
             split_missing: An optional flag indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings; if False, then missing data is ignored (i.e., all states are 0). Default value is True.
         """
         # Convert the collation to a Pandas DataFrame first:
-        df = self.to_dataframe(long_table, split_missing)
+        df = self.to_dataframe(drop_constant=drop_constant, long_table=long_table, split_missing=split_missing)
         # If this is a long table, then do not include row indices:
         if long_table:
             return df.to_excel(file_addr, index=False)
@@ -990,39 +1092,45 @@ class Collation:
         Args:
             file_addr: A string representing the path to an output STEMMA prep file; the file should have no extension.
             The accompanying chron file will match this file name, except that it will have "_chron" appended to the end.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
         """
+        # Populate a list of sites that will correspond to columns of the sequence alignment
+        # (by default, constant sites are dropped):
+        substantive_variation_unit_ids = [
+            vu_id for vu_id in self.variation_unit_ids if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
+        ]
+        substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
+        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
         # In a first pass, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
         # to the readings' texts:
         reading_texts_by_indices = {}
-        substantive_variation_unit_ids_set = set(self.substantive_variation_unit_ids)
-        substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
-        vu_ind = 0
-        for vu in self.variation_units:
+        for j, vu in enumerate(self.variation_units):
             if vu.id not in substantive_variation_unit_ids_set:
                 continue
-            rdg_ind = 0
+            k = 0
             for rdg in vu.readings:
                 key = tuple([vu.id, rdg.id])
                 if key not in substantive_variation_unit_reading_tuples_set:
                     continue
-                indices = tuple([vu_ind, rdg_ind])
+                indices = tuple([j, k])
                 reading_texts_by_indices[indices] = rdg.text
-                rdg_ind += 1
-            if rdg_ind > 0:
-                vu_ind += 1
+                k += 1
         # In a second pass, populate another dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
         # to the witnesses exclusively supporting those readings:
         reading_wits_by_indices = {}
         for indices in reading_texts_by_indices:
             reading_wits_by_indices[indices] = []
-        for wit in self.readings_by_witness:
-            for vu_ind, rdg_support in enumerate(self.readings_by_witness[wit]):
+        for i, wit in enumerate(self.witnesses):
+            for j, vu_id in enumerate(self.variation_unit_ids):
+                if vu_id not in substantive_variation_unit_ids_set:
+                    continue
+                rdg_support = self.readings_by_witness[wit.id][j]
                 # If this witness does not exclusively support exactly one reading at this unit, then treat it as lacunose:
-                if len([i for i, w in enumerate(rdg_support) if w > 0]) != 1:
+                if len([k for k, w in enumerate(rdg_support) if w > 0]) != 1:
                     continue
-                rdg_ind = rdg_support.index(1)
-                indices = tuple([vu_ind, rdg_ind])
-                reading_wits_by_indices[indices].append(wit)
+                k = rdg_support.index(1)
+                indices = tuple([j, k])
+                reading_wits_by_indices[indices].append(wit.id)
         # In a third pass, write to the STEMMA file:
         chron_file_addr = str(file_addr) + "_chron"
         with open(file_addr, "w", encoding="utf-8") as f:
@@ -1032,14 +1140,16 @@ class Collation:
             # Then add a line indicating that all witnesses are lacunose unless they are specified explicitly:
             f.write("= $? $* ;\n\n")
             # Then proceed for each variation unit:
-            for vu_ind, vu_id in enumerate(self.substantive_variation_unit_ids):
+            for j, vu_id in enumerate(self.variation_unit_ids):
+                if vu_id not in substantive_variation_unit_ids_set:
+                    continue
                 # Print the variation unit ID first:
                 f.write("@ %s\n" % vu_id)
                 # In a first pass, print the texts of all readings enclosed in brackets:
                 f.write("[ ")
-                rdg_ind = 0
+                k = 0
                 while True:
-                    indices = tuple([vu_ind, rdg_ind])
+                    indices = tuple([j, k])
                     if indices not in reading_texts_by_indices:
                         break
                     text = slugify(
@@ -1049,29 +1159,29 @@ class Collation:
                     if text == "":
                         text = "\u2013"
                     # The first reading should not be preceded by anything:
-                    if rdg_ind == 0:
+                    if k == 0:
                         f.write(text)
                         f.write(" |")
                     # Every subsequent reading should be preceded by a space:
-                    elif rdg_ind > 0:
+                    elif k > 0:
                         f.write(" %s" % text)
-                    rdg_ind += 1
+                    k += 1
                 f.write(" ]\n")
                 # In a second pass, print the indices and witnesses for all readings enclosed in angle brackets:
-                rdg_ind = 0
+                k = 0
                 f.write("\t< ")
                 while True:
-                    indices = tuple([vu_ind, rdg_ind])
+                    indices = tuple([j, k])
                     if indices not in reading_wits_by_indices:
                         break
                     wits = " ".join(reading_wits_by_indices[indices])
                     # Open the variant reading support block with an angle bracket:
-                    if rdg_ind == 0:
-                        f.write("%d %s" % (rdg_ind, wits))
+                    if k == 0:
+                        f.write("%d %s" % (k, wits))
                     # Open all subsequent variant reading support blocks with pipes on the next line:
                     else:
-                        f.write("\n\t| %d %s" % (rdg_ind, wits))
-                    rdg_ind += 1
+                        f.write("\n\t| %d %s" % (k, wits))
+                    k += 1
                 f.write(" >\n")
         # In a fourth pass, write to the chron file:
         max_id_length = max(
@@ -1122,6 +1232,7 @@ class Collation:
         self,
         file_addr: Union[Path, str],
         format: Format = None,
+        drop_constant: bool = False,
         long_table: bool = False,
         split_missing: bool = True,
         char_state_labels: bool = True,
@@ -1137,6 +1248,7 @@ class Collation:
             format (Format, optional): The desired output format.
                 If None then it is infered from the file suffix.
                 Defaults to None.
+            drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             long_table (bool, optional): An optional flag indicating whether or not to generate a long table
                 with columns for taxa, characters, reading indices, and reading values.
                 Not applicable for NEXUS, HENNIG86, PHYLIP, FASTA, or STEMMA format.
@@ -1174,6 +1286,7 @@ class Collation:
         if format == Format.NEXUS:
             return self.to_nexus(
                 file_addr,
+                drop_constant=drop_constant,
                 char_state_labels=char_state_labels,
                 frequency=frequency,
                 ambiguous_as_missing=ambiguous_as_missing,
@@ -1182,22 +1295,28 @@ class Collation:
             )
 
         if format == format.HENNIG86:
-            return self.to_hennig86(file_addr)
+            return self.to_hennig86(file_addr, drop_constant=drop_constant)
 
         if format == format.PHYLIP:
-            return self.to_phylip(file_addr)
+            return self.to_phylip(file_addr, drop_constant=drop_constant)
 
         if format == format.FASTA:
-            return self.to_fasta(file_addr)
+            return self.to_fasta(file_addr, drop_constant=drop_constant)
 
         if format == Format.CSV:
-            return self.to_csv(file_addr, long_table=long_table, split_missing=split_missing)
+            return self.to_csv(
+                file_addr, drop_constant=drop_constant, long_table=long_table, split_missing=split_missing
+            )
 
         if format == Format.TSV:
-            return self.to_csv(file_addr, long_table=long_table, split_missing=split_missing, sep="\t")
+            return self.to_csv(
+                file_addr, drop_constant=drop_constant, long_table=long_table, split_missing=split_missing, sep="\t"
+            )
 
         if format == Format.EXCEL:
-            return self.to_excel(file_addr, long_table=long_table, split_missing=split_missing)
+            return self.to_excel(
+                file_addr, drop_constant=drop_constant, long_table=long_table, split_missing=split_missing
+            )
 
         if format == Format.STEMMA:
             return self.to_stemma(file_addr)
diff --git a/teiphy/main.py b/teiphy/main.py
index 91c494c..bbecfd3 100644
--- a/teiphy/main.py
+++ b/teiphy/main.py
@@ -47,6 +47,10 @@ def to_file(
         False,
         help="Use the StatesFormat=Frequency setting instead of the StatesFormat=StatesPresent setting (and thus represent all states with frequency vectors rather than symbols) in NEXUS output.",
     ),
+    drop_constant: bool = typer.Option(
+        False,
+        help="If set, do not write constant sites (i.e., variation units with one substantive reading) to output.",
+    ),
     ambiguous_as_missing: bool = typer.Option(
         False,
         help="Use the missing symbol instead of multistate symbols (and thus treat all ambiguities as missing data) in NEXUS output; this option is only applied if the --frequency option is not set.",
@@ -111,6 +115,7 @@ def to_file(
     coll.to_file(
         output,
         format=format,
+        drop_constant=drop_constant,
         char_state_labels=labels,
         frequency=frequency,
         ambiguous_as_missing=ambiguous_as_missing,


Test patch:
diff --git a/tests/test_collation.py b/tests/test_collation.py
index 24d63a5..7a7630f 100644
--- a/tests/test_collation.py
+++ b/tests/test_collation.py
@@ -5,7 +5,7 @@ from pathlib import Path
 import numpy as np
 from lxml import etree as et
 
-from teiphy import Collation
+from teiphy import tei_ns, Collation
 
 root_dir = Path("__file__").parent.parent
 input_example = root_dir / "example/ubs_ephesians.xml"
@@ -15,99 +15,85 @@ class CollationDefaultTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_witnesses = xml.xpath("//tei:listWit/tei:witness", namespaces={"tei": tei_ns})
+        self.xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+        self.xml_readings = xml.xpath("//tei:rdg", namespaces={"tei": tei_ns})
         self.collation = Collation(xml)
 
     def test_witnesses(self):
-        self.assertEqual(len(self.collation.witnesses), 38)
+        self.assertEqual(len(self.collation.witnesses), len(self.xml_witnesses))
 
     def test_witness_index_by_id(self):
-        self.assertEqual(len(self.collation.witness_index_by_id), 38)
+        self.assertEqual(len(self.collation.witness_index_by_id), len(self.xml_witnesses))
 
     def test_variation_units(self):
-        self.assertEqual(len(self.collation.variation_units), 42)
+        self.assertEqual(len(self.collation.variation_units), len(self.xml_variation_units))
 
     def test_readings_by_witness(self):
-        self.assertEqual(len(self.collation.readings_by_witness), 38)
+        self.assertEqual(len(self.collation.readings_by_witness), len(self.xml_witnesses))
 
     def test_substantive_variation_unit_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_ids), 40
-        )  # apps B10K4V28U24-26 and B10K6V20U12 should always be ignored
+        self.assertEqual(len(self.collation.variation_unit_ids), len(self.xml_variation_units))
 
-    def test_substantive_reading_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_reading_tuples), 160
-        )  # all readings except for the ones in apps B10K4V28U24-26 and B10K6V20U12
+    def test_substantive_readings_by_variation_unit_id(self):
+        self.assertEqual(len(self.collation.substantive_readings_by_variation_unit_id), len(self.xml_variation_units))
+
+    def test_substantive_variation_unit_reading_tuples(self):
+        self.assertEqual(len(self.collation.substantive_variation_unit_reading_tuples), len(self.xml_readings))
 
 
 class CollationTrivialReconstructedTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_readings = xml.xpath("//tei:rdg[not(@type) or @type!=\"reconstructed\"]", namespaces={"tei": tei_ns})
         self.collation = Collation(xml, trivial_reading_types=["reconstructed"])
 
-    def test_substantive_variation_unit_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_ids), 40
-        )  # apps B10K4V28U24-26 and B10K6V20U12 should always be ignored
-
-    def test_substantive_reading_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_reading_tuples), 147
-        )  # all readings except for the ones in apps B10K4V28U24-26 and B10K6V20U12 and all reconstructed readings
+    def test_substantive_variation_unit_reading_tuples(self):
+        self.assertEqual(len(self.collation.substantive_variation_unit_reading_tuples), len(self.xml_readings))
 
 
 class CollationTrivialDefectiveTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_readings = xml.xpath(
+            "//tei:rdg[not(@type) or (@type!=\"reconstructed\" and @type!=\"defective\")]", namespaces={"tei": tei_ns}
+        )
         self.collation = Collation(xml, trivial_reading_types=["reconstructed", "defective"])
 
-    def test_substantive_variation_unit_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_ids), 40
-        )  # apps B10K4V28U24-26 and B10K6V20U12 should always be ignored
-
-    def test_substantive_reading_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_reading_tuples), 108
-        )  # all readings except for the ones in apps B10K4V28U24-26 and B10K6V20U12 and all reconstructed and defective readings
+    def test_substantive_variation_unit_reading_tuples(self):
+        self.assertEqual(len(self.collation.substantive_variation_unit_reading_tuples), len(self.xml_readings))
 
 
 class CollationTrivialOrthographicTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_readings = xml.xpath(
+            "//tei:rdg[not(@type) or (@type!=\"reconstructed\" and @type!=\"defective\" and @type !=\"orthographic\")]",
+            namespaces={"tei": tei_ns},
+        )
         self.collation = Collation(xml, trivial_reading_types=["reconstructed", "defective", "orthographic"])
 
-    def test_substantive_variation_unit_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_ids), 40
-        )  # all units except B10K6V20U12 and B10K4V28U24-26
-
-    def test_substantive_reading_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_reading_tuples), 105
-        )  # all readings except for the ones in app B10K6V20U12 and B10K4V28U24-26 and all reconstructed, defective, and orthographic readings
+    def test_substantive_variation_unit_reading_tuples(self):
+        self.assertEqual(len(self.collation.substantive_variation_unit_reading_tuples), len(self.xml_readings))
 
 
 class CollationTrivialSubreadingTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_readings = xml.xpath(
+            "//tei:rdg[not(@type) or (@type!=\"reconstructed\" and @type!=\"defective\" and @type !=\"orthographic\" and @type !=\"subreading\")]",
+            namespaces={"tei": tei_ns},
+        )
         self.collation = Collation(
             xml, trivial_reading_types=["reconstructed", "defective", "orthographic", "subreading"]
         )
 
-    def test_substantive_variation_unit_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_ids), 40
-        )  # all units except B10K6V20U12 and B10K4V28U24-26
-
-    def test_substantive_reading_ids(self):
-        self.assertEqual(
-            len(self.collation.substantive_variation_unit_reading_tuples), 98
-        )  # all readings except for the ones in app B10K6V20U12 and B10K4V28U24-26 and all readings that have non-substantive types
+    def test_substantive_variation_unit_reading_tuples(self):
+        self.assertEqual(len(self.collation.substantive_variation_unit_reading_tuples), len(self.xml_readings))
 
 
 class CollationMissingTestCase(unittest.TestCase):
@@ -117,14 +103,14 @@ class CollationMissingTestCase(unittest.TestCase):
         self.collation = Collation(xml, missing_reading_types=["lac", "overlap"])
 
     def test_missing_lac(self):
-        vu_ind = self.collation.substantive_variation_unit_ids.index("B10K1V1U24-26")
+        vu_ind = self.collation.variation_unit_ids.index("B10K1V1U24-26")
         rdg_support = self.collation.readings_by_witness["04"][vu_ind]
         self.assertEqual(
             sum(rdg_support), 0
         )  # all entries in the reading support vector for this lacunose witness should be 0
 
     def test_missing_overlap(self):
-        vu_ind = self.collation.substantive_variation_unit_ids.index("B10K3V20U8-10")
+        vu_ind = self.collation.variation_unit_ids.index("B10K3V20U8-10")
         rdg_support = self.collation.readings_by_witness["606"][vu_ind]
         self.assertEqual(
             sum(rdg_support), 0
@@ -160,7 +146,7 @@ class CollationManuscriptSuffixesTestCase(unittest.TestCase):
         self.assertEqual(self.collation.get_base_wit("424T*"), "424")
 
     def test_merged_attestations(self):
-        vu_ind = self.collation.substantive_variation_unit_ids.index("B10K3V14U14-18")
+        vu_ind = self.collation.variation_unit_ids.index("B10K3V14U14-18")
         rdg_support = self.collation.readings_by_witness["1910"][vu_ind]
         self.assertEqual(
             rdg_support, [0.5, 0.5, 0]
@@ -176,14 +162,14 @@ class CollationFillCorrectorLacunaeTestCase(unittest.TestCase):
         )
 
     def test_inactive_corrector(self):
-        vu_ind = self.collation.substantive_variation_unit_ids.index("B10K4V8U12-18")
+        vu_ind = self.collation.variation_unit_ids.index("B10K4V8U12-18")
         rdg_support = self.collation.readings_by_witness["06C1"][vu_ind]
         self.assertEqual(
             rdg_support, [0, 0, 1, 0, 0, 0, 0, 0, 0]
         )  # this corrector is inactive in this unit and should default to the first-hand reading
 
     def test_active_corrector(self):
-        vu_ind = self.collation.substantive_variation_unit_ids.index("B10K4V8U12-18")
+        vu_ind = self.collation.variation_unit_ids.index("B10K4V8U12-18")
         rdg_support = self.collation.readings_by_witness["06C2"][vu_ind]
         self.assertEqual(
             rdg_support, [0, 0, 0, 0, 0, 1, 0, 0, 0]
@@ -203,6 +189,7 @@ class CollationOutputTestCase(unittest.TestCase):
     def setUp(self):
         parser = et.XMLParser(remove_comments=True)
         xml = et.parse(input_example, parser=parser)
+        self.xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
         self.collation = Collation(
             xml,
             trivial_reading_types=["reconstructed", "defective", "orthographic", "subreading"],
@@ -271,13 +258,13 @@ class CollationOutputTestCase(unittest.TestCase):
     def test_to_numpy_ignore_missing(self):
         matrix, reading_labels, witness_labels = self.collation.to_numpy(split_missing=False)
         self.assertTrue(
-            matrix.sum(axis=0)[5] < len(self.collation.substantive_variation_unit_ids)
+            matrix.sum(axis=0)[5] < len(self.collation.variation_unit_ids)
         )  # lacuna in the first witness should result in its column summing to less than the total number of substantive variation units
 
     def test_to_numpy_split_missing(self):
         matrix, reading_labels, witness_labels = self.collation.to_numpy(split_missing=True)
         self.assertTrue(
-            abs(matrix.sum(axis=0)[5] - len(self.collation.substantive_variation_unit_ids) < 1e-4)
+            abs(matrix.sum(axis=0)[5] - len(self.collation.variation_unit_ids)) < 1e-4
         )  # the column for the first witness should sum to the total number of substantive variation units (give or take some rounding error)
 
     def test_to_distance_matrix(self):
@@ -286,7 +273,15 @@ class CollationOutputTestCase(unittest.TestCase):
         self.assertTrue(np.all(matrix == matrix.T))  # matrix should be symmetrical
         self.assertEqual(
             matrix[0, 1], 13
-        )  # entry for UBS and P46 should be 13 (remember not to count P46 lacunae and ambiguities, P46 defective readings under the UBS reading, and non-substantive variation units)
+        )  # entry for UBS and P46 should be 13 (remember not to count P46 lacunae and ambiguities and P46 defective readings under the UBS reading)
+
+    def test_to_distance_matrix_drop_constant(self):
+        matrix, witness_labels = self.collation.to_distance_matrix(drop_constant=True)
+        self.assertEqual(np.trace(matrix), 0)  # diagonal entries should be 0
+        self.assertTrue(np.all(matrix == matrix.T))  # matrix should be symmetrical
+        self.assertEqual(
+            matrix[0, 1], 13
+        )  # entry for UBS and P46 should be 13 (remember not to count P46 lacunae and ambiguities and P46 defective readings under the UBS reading)
 
     def test_to_distance_matrix_proportion(self):
         matrix, witness_labels = self.collation.to_distance_matrix(proportion=True)
@@ -294,8 +289,17 @@ class CollationOutputTestCase(unittest.TestCase):
         self.assertTrue(np.all(matrix == matrix.T))  # matrix should be symmetrical
         self.assertTrue(np.all(matrix >= 0.0) and np.all(matrix <= 1.0))  # all elements should be between 0 and 1
         self.assertTrue(
-            abs(matrix[0, 1] - 13 / 38) < 1e-4
-        )  # entry for UBS and P46 should be close to 13/38 (of 40 substantive variation units, P46 is lacunose at one and ambiguous at another)
+            abs(matrix[0, 1] - 13 / (len(self.xml_variation_units) - 2)) < 1e-4
+        )  # entry for UBS and P46 should be 13 divided by the number of variation units where neither witness is lacunose or ambiguous
+
+    def test_to_distance_matrix_drop_constant_proportion(self):
+        matrix, witness_labels = self.collation.to_distance_matrix(drop_constant=True, proportion=True)
+        self.assertEqual(np.trace(matrix), 0)  # diagonal entries should be 0
+        self.assertTrue(np.all(matrix == matrix.T))  # matrix should be symmetrical
+        self.assertTrue(np.all(matrix >= 0.0) and np.all(matrix <= 1.0))  # all elements should be between 0 and 1
+        self.assertTrue(
+            abs(matrix[0, 1] - 13 / (len(self.xml_variation_units) - 2 - 2)) < 1e-4
+        )  # entry for UBS and P46 should be 13 divided by the number of non-constant variation units where neither witness is lacunose or ambiguous
 
     def test_to_long_table(self):
         long_table, column_labels = self.collation.to_long_table()
diff --git a/tests/test_main.py b/tests/test_main.py
index 454cd2f..88f8726 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -1,9 +1,11 @@
 from pathlib import Path
 import tempfile
 from typer.testing import CliRunner
+from lxml import etree as et
 
 from teiphy.main import app
 from teiphy.collation import ParsingException
+from teiphy.common import tei_ns
 
 runner = CliRunner()
 
@@ -66,6 +68,20 @@ def test_to_nexus():
         text = output.read_text(encoding="utf-8")
         assert text.startswith("#NEXUS")
         assert "CharStateLabels" in text
+        assert "22 B10K4V28U24_26" in text
+        assert "StatesFormat=Frequency" not in text
+
+
+def test_to_nexus_drop_constant():
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.nexus"
+        result = runner.invoke(app, ["--drop-constant", str(input_example), str(output)])
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="utf-8")
+        assert text.startswith("#NEXUS")
+        assert "CharStateLabels" in text
+        assert "22 B10K4V28U24_26" not in text
         assert "StatesFormat=Frequency" not in text
 
 
@@ -89,6 +105,19 @@ def test_to_nexus_frequency():
         assert output.exists()
         text = output.read_text(encoding="utf-8")
         assert text.startswith("#NEXUS")
+        assert "22 B10K4V28U24_26" in text
+        assert "StatesFormat=Frequency" in text
+
+
+def test_to_nexus_drop_constant_frequency():
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.nexus"
+        result = runner.invoke(app, ["--drop-constant", "--frequency", str(input_example), str(output)])
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="utf-8")
+        assert text.startswith("#NEXUS")
+        assert "22 B10K4V28U24_26" not in text
         assert "StatesFormat=Frequency" in text
 
 
@@ -186,6 +215,11 @@ def test_to_nexus_mrbayes_some_dates():
 
 
 def test_to_hennig86():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_witnesses = xml.xpath("//tei:listWit/tei:witness", namespaces={"tei": tei_ns})
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
     with tempfile.TemporaryDirectory() as tmp_dir:
         output = Path(tmp_dir) / "test.tnt"
         result = runner.invoke(
@@ -208,9 +242,76 @@ def test_to_hennig86():
         text = output.read_text(encoding="utf-8")
         assert text.startswith("nstates")
         assert "xread" in text
+        assert "%d %d" % (len(xml_variation_units), len(xml_witnesses)) in text
+
+
+def test_to_hennig86_drop_constant():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_witnesses = xml.xpath("//tei:listWit/tei:witness", namespaces={"tei": tei_ns})
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.tnt"
+        result = runner.invoke(
+            app,
+            [
+                "-t reconstructed",
+                "-t defective",
+                "-t orthographic",
+                "-m lac",
+                "-m overlap",
+                "-s *",
+                "-s T",
+                "--fill-correctors",
+                "--drop-constant",
+                str(input_example),
+                str(output),
+            ],
+        )
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="utf-8")
+        assert text.startswith("nstates")
+        assert "xread" in text
+        assert "%d %d" % (len(xml_variation_units) - 2, len(xml_witnesses)) in text
 
 
 def test_to_phylip():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_witnesses = xml.xpath("//tei:listWit/tei:witness", namespaces={"tei": tei_ns})
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.phy"
+        result = runner.invoke(
+            app,
+            [
+                "-t reconstructed",
+                "-t defective",
+                "-t orthographic",
+                "-m lac",
+                "-m overlap",
+                "-s *",
+                "-s T",
+                "--fill-correctors",
+                str(input_example),
+                str(output),
+            ],
+        )
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="ascii")
+        assert text.startswith("%d %d" % (len(xml_witnesses), len(xml_variation_units)))
+
+
+def test_to_phylip_drop_constant():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_witnesses = xml.xpath("//tei:listWit/tei:witness", namespaces={"tei": tei_ns})
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
     with tempfile.TemporaryDirectory() as tmp_dir:
         output = Path(tmp_dir) / "test.phy"
         result = runner.invoke(
@@ -224,6 +325,7 @@ def test_to_phylip():
                 "-s *",
                 "-s T",
                 "--fill-correctors",
+                "--drop-constant",
                 str(input_example),
                 str(output),
             ],
@@ -231,10 +333,14 @@ def test_to_phylip():
         assert result.exit_code == 0
         assert output.exists()
         text = output.read_text(encoding="ascii")
-        assert text.startswith("38 40")
+        assert text.startswith("%d %d" % (len(xml_witnesses), len(xml_variation_units) - 2))
 
 
 def test_to_fasta():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
     with tempfile.TemporaryDirectory() as tmp_dir:
         output = Path(tmp_dir) / "test.fa"
         result = runner.invoke(
@@ -256,6 +362,37 @@ def test_to_fasta():
         assert output.exists()
         text = output.read_text(encoding="ascii")
         assert text.startswith(">UBS")
+        assert "0" * len(xml_variation_units) in text
+
+
+def test_to_fasta_drop_constant():
+    parser = et.XMLParser(remove_comments=True)
+    xml = et.parse(input_example, parser=parser)
+    xml_variation_units = xml.xpath("//tei:app", namespaces={"tei": tei_ns})
+
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.fa"
+        result = runner.invoke(
+            app,
+            [
+                "-t reconstructed",
+                "-t defective",
+                "-t orthographic",
+                "-m lac",
+                "-m overlap",
+                "-s *",
+                "-s T",
+                "--fill-correctors",
+                "--drop-constant",
+                str(input_example),
+                str(output),
+            ],
+        )
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="ascii")
+        assert text.startswith(">UBS")
+        assert "0" * (len(xml_variation_units) - 2) in text
 
 
 def test_to_csv():
@@ -266,6 +403,18 @@ def test_to_csv():
         assert output.exists()
         text = output.read_text(encoding="utf-8")
         assert text.startswith(",UBS,P46,01,02,03,04,06")
+        assert "B10K4V28U24-26," in text
+
+
+def test_to_csv_drop_constant():
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.csv"
+        result = runner.invoke(app, ["--drop-constant", str(input_example), str(output)])
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="utf-8")
+        assert text.startswith(",UBS,P46,01,02,03,04,06")
+        assert "B10K4V28U24-26," not in text
 
 
 def test_to_csv_long_table():
@@ -276,6 +425,18 @@ def test_to_csv_long_table():
         assert output.exists()
         text = output.read_text(encoding="utf-8")
         assert text.startswith("taxon,character,state,value")
+        assert "B10K4V28U24-26," in text
+
+
+def test_to_csv_drop_constant_long_table():
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        output = Path(tmp_dir) / "test.csv"
+        result = runner.invoke(app, ["--drop-constant", "--long-table", str(input_example), str(output)])
+        assert result.exit_code == 0
+        assert output.exists()
+        text = output.read_text(encoding="utf-8")
+        assert text.startswith("taxon,character,state,value")
+        assert "B10K4V28U24-26," not in text
 
 
 def test_to_tsv():