Better cut dataframes (cut at Rfam mapping)

Louis BECQUEY
Commit c21f55fc9d38da93a83baa819981dbb141cc8783 c21f55fc 1 parent 98c2907a
Showing 4 changed files with 128 additions and 50 deletions
.gitignore
README.md
RNAnet.py
statistics.py
--- a/.gitignore
View file @c21f55f
+++ b/.gitignore
View file @c21f55f
@@ -7,6 +7,7 @@ results/
 # temporary results files
 data/
+esl*
 # environment stuff
 .vscode/
--- a/README.md
View file @c21f55f
+++ b/README.md
View file @c21f55f
@@ -49,6 +49,19 @@ Other folders are created and not deleted, which you might want to conserve to a
 * `path-to-3D-folder-you-passed-in-option/annotations/` contains the raw JSON annotation files of the previous mmCIF structures. You may find additional information into them which is not properly supported by RNANet yet.
 # How to run
+
+## Required computational resources
+- CPU: no requirements. The program is optimized for multi-core CPUs, you might want to use Intel Xeons, AMD Ryzens, etc.
+- GPU: not required
+- RAM: 16 GB with a large swap partition is okay. 32 GB is recommended (usage peaks at ~27 GB)
+- Storage: to date, it takes 60 GB for the 3D data (36 GB if you don't use the --extract option), 11 GB for the sequence data, and 7GB for the outputs (5.6 GB database, 1 GB archive of CSV files). You need to add a few more for the dependencies. Go for a 100GB partition and you are good to go. The computation speed is really decreased if you use a fast storage device (e.g. SSD instead of hard drive, or even better, a NVMe SSD) because of permanent I/O with the SQlite database.
+- Network : We query the Rfam public MySQL server on port 4497. Make sure your network enables communication (there should not be any issue on private networks, but maybe you company/university closes ports by default). You will get an error message if the port is not open. Around 30 GB of data is downloaded.
+
+To give you an estimation, our last full run took exactly 12h, excluding the time to download the MMCIF files containing RNA (around 25GB to download) and the time to compute statistics.
+Measured the 23rd of June 2020 on a 16-core AMD Ryzen 7 3700X CPU @3.60GHz, plus 32 Go RAM, and a 7200rpm Hard drive. Total CPU time spent: 135 hours (user+kernel modes), corresponding to 12h (actual time spent with the 16-core CPU). 
+
+Update runs are much quicker, around 3 hours. It depends mostly on what RNA families are concerned by the update.
+
 ## Dependencies
 You need to install:
 - DSSR, you need to register to the X3DNA forum [here](http://forum.x3dna.org/site-announcements/download-instructions/) and then download the DSSR binary [on that page](http://forum.x3dna.org/downloads/3dna-download/). 
@@ -91,6 +104,8 @@ The detailed list of options is below:
 ## Post-computation task: estimate quality
 The file statistics.py is supposed to give a summary on the produced dataset. See the results/ folder. It can be run automatically after RNANet if you pass the `-s` option.
+
+
 # How to further filter the dataset
 You may want to build your own sub-dataset by querying the results/RNANet.db file. Here are quick examples using Python3 and its sqlite3 package.
@@ -108,6 +123,7 @@ Step 1 : We first get a list of chains that are below our favorite resolution th
 with sqlite3.connect("results/RNANet.db) as connection:
     chain_list = pd.read_sql("""SELECT chain_id, structure_id, chain_name
                                 FROM chain JOIN structure 
+                                ON chain.structure_id = structure.pdb_id
                                 WHERE resolution < 4.0 
                                 ORDER BY structure_id ASC;""",
                             con=connection)
@@ -146,6 +162,7 @@ We will simply modify the Step 1 above:
 with sqlite3.connect("results/RNANet.db) as connection:
     chain_list = pd.read_sql("""SELECT chain_id, structure_id, chain_name
                                 FROM chain JOIN structure 
+                                ON chain.structure_id = structure.pdb_id
                                 WHERE date < "2018-06-01" 
                                 ORDER BY structure_id ASC;""",
                             con=connection)
@@ -160,6 +177,7 @@ If you want just one example of each RNA 3D chain, use in Step 1:
 with sqlite3.connect("results/RNANet.db) as connection:
     chain_list = pd.read_sql("""SELECT UNIQUE chain_id, structure_id, chain_name
                                 FROM chain JOIN structure
+                                ON chain.structure_id = structure.pdb_id
                                 ORDER BY structure_id ASC;""",
                             con=connection)
 ```
--- a/RNAnet.py
View file @c21f55f
+++ b/RNAnet.py
View file @c21f55f
@@ -6,14 +6,15 @@ from Bio import AlignIO, SeqIO
 from Bio.PDB import MMCIFParser
 from Bio.PDB.mmcifio import MMCIFIO
 from Bio.PDB.MMCIF2Dict import MMCIF2Dict 
-from Bio.PDB.PDBExceptions import PDBConstructionWarning
+from Bio.PDB.PDBExceptions import PDBConstructionWarning, BiopythonWarning
+from Bio.PDB.Dice import ChainSelector
 from Bio._py3k import urlretrieve as _urlretrieve
 from Bio._py3k import urlcleanup as _urlcleanup
 from Bio.Alphabet import generic_rna
 from Bio.Seq import Seq
 from Bio.SeqRecord import SeqRecord
 from Bio.Align import MultipleSeqAlignment, AlignInfo
-from collections import OrderedDict
+from collections import OrderedDict, defaultdict
 from functools import partial, wraps
 from os import path, makedirs
 from multiprocessing import Pool, Manager, set_start_method
@@ -40,6 +41,7 @@ errsymb = '\U0000274C'
 LSU_set = {"RF00002", "RF02540", "RF02541", "RF02543", "RF02546"}   # From Rfam CLAN 00112
 SSU_set = {"RF00177", "RF02542",  "RF02545", "RF01959", "RF01960"}  # From Rfam CLAN 00111
+no_nts_set = set()
 class NtPortionSelector(object):
     """Class passed to MMCIFIO to select some chain portions in an MMCIF file.
@@ -170,29 +172,40 @@ class Chain:
         with warnings.catch_warnings():
             # Ignore the PDB problems. This mostly warns that some chain is discontinuous.
             warnings.simplefilter('ignore', PDBConstructionWarning)  
+            warnings.simplefilter('ignore', BiopythonWarning)  
             # Load the whole mmCIF into a Biopython structure object:
             mmcif_parser = MMCIFParser()
             s = mmcif_parser.get_structure(self.pdb_id, path_to_3D_data + "RNAcifs/"+self.pdb_id+".cif")
+            
             # Extract the desired chain
             c = s[model_idx][self.pdb_chain_id]
             if (self.pdb_end - self.pdb_start):
-                # Pay attention to residue numbering
+                # # Pay attention to residue numbering 
-                first_number = c.child_list[0].get_id()[1]          # the chain's first residue is numbered 'first_number'
+                # first_number = c.child_list[0].get_id()[1]          # the chain's first residue is numbered 'first_number'
+                # if self.pdb_start < self.pdb_end:                             
+                #     start = self.pdb_start + first_number - 1       # shift our start_position by 'first_number'
+                #     end = self.pdb_end + first_number - 1           # same for the end position
+                # else:
+                #     self.reversed = True                            # the 3D chain is numbered backwards compared to the Rfam family
+                #     end = self.pdb_start + first_number - 1
+                #     start = self.pdb_end + first_number - 1
+
                 if self.pdb_start < self.pdb_end:                             
-                    start = self.pdb_start + first_number - 1       # shift our start_position by 'first_number'
+                    start = self.pdb_start        # shift our start_position by 'first_number'
-                    end = self.pdb_end + first_number - 1           # same for the end position
+                    end = self.pdb_end            # same for the end position
                 else:
-                    self.reversed = True                            # the 3D chain is numbered backwards compared to the Rfam family
+                    self.reversed = True          # the 3D chain is numbered backwards compared to the Rfam family
-                    end = self.pdb_start + first_number - 1
+                    end = self.pdb_start
-                    start = self.pdb_end + first_number - 1
+                    start = self.pdb_end
             else:
                 start = c.child_list[0].get_id()[1]  # the chain's first residue is numbered 'first_number'
                 end = c.child_list[-1].get_id()[1]   # the chain's last residue number
             # Define a selection
+            # sel = ChainSelector(self.pdb_chain_id, start, end, model_id = model_idx)
             sel = NtPortionSelector(model_idx, self.pdb_chain_id, start, end, khetatm)
             # Save that selection on the mmCIF object s to file
@@ -205,7 +218,9 @@ class Chain:
     def extract_3D_data(self):
         """ Maps DSSR annotations to the chain. """
+        ############################################
         # Load the mmCIF annotations from file
+        ############################################
         try:
             with open(path_to_3D_data + "annotations/" + self.pdb_id + ".json", 'r') as json_file:
                 json_object = json.load(json_file)
@@ -219,39 +234,53 @@ class Chain:
         # Print eventual warnings given by DSSR, and abort if there are some
         if "warning" in json_object.keys():
             warn(f"found DSSR warning in annotation {self.pdb_id}.json: {json_object['warning']}. Ignoring {self.chain_label}.")
+            if "no nucleotides" in json_object['warning']:
+                no_nts_set.add(self.pdb_id)
             self.delete_me = True
             self.error_messages = f"DSSR warning {self.pdb_id}.json: {json_object['warning']}. Ignoring {self.chain_label}."
             return 1
+        ############################################
+        # Create the data-frame
+        ############################################
         try:
-            # Prepare a data structure (Pandas DataFrame) for the nucleotides
+            # Create the Pandas DataFrame for the nucleotides of the right chain
             nts = json_object["nts"]                         # sub-json-object
             df = pd.DataFrame(nts)                           # conversion to dataframe
             df = df[ df.chain_name == self.pdb_chain_id ]    # keeping only this chain's nucleotides
+
+            # Assert nucleotides of the chain are found
             if df.empty:
                 warn(f"Could not find nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. Ignoring chain {self.chain_label}.", error=True)
+                no_nts_set.add(self.pdb_id)
                 self.delete_me = True
                 self.error_messages = f"Could not find nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. We expect a problem with {self.pdb_id} mmCIF download. Delete it and retry."
                 return 1
-            # remove low pertinence or undocumented descriptors
+            # Remove low pertinence or undocumented descriptors, convert angles values
             cols_we_keep = ["index_chain", "nt_resnum", "nt_name", "nt_code", "nt_id", "dbn", "alpha", "beta", "gamma", "delta", "epsilon", "zeta",
                 "epsilon_zeta", "bb_type", "chi", "glyco_bond", "form", "ssZp", "Dp", "eta", "theta", "eta_prime", "theta_prime", "eta_base", "theta_base",
                 "v0", "v1", "v2", "v3", "v4", "amplitude", "phase_angle", "puckering" ]
             df = df[cols_we_keep]
-
+            df.loc[:,['alpha', 'beta','gamma','delta','epsilon','zeta','epsilon_zeta','chi','v0', 'v1', 'v2', 'v3', 'v4', # Conversion to radians
-            # Convert angles to radians
-            df.loc[:,['alpha', 'beta','gamma','delta','epsilon','zeta','epsilon_zeta','chi','v0', 'v1', 'v2', 'v3', 'v4',
                         'eta','theta','eta_prime','theta_prime','eta_base','theta_base', 'phase_angle']] *= np.pi/180.0
-            # mapping [-pi, pi] into [0, 2pi]
+            df.loc[:,['alpha', 'beta','gamma','delta','epsilon','zeta','epsilon_zeta','chi','v0', 'v1', 'v2', 'v3', 'v4', # mapping [-pi, pi] into [0, 2pi]
-            df.loc[:,['alpha', 'beta','gamma','delta','epsilon','zeta','epsilon_zeta','chi','v0', 'v1', 'v2', 'v3', 'v4',
                         'eta','theta','eta_prime','theta_prime','eta_base','theta_base', 'phase_angle']] %= (2.0*np.pi)
+
         except KeyError as e:
             warn(f"Error while parsing DSSR {self.pdb_id}.json output:{e}", error=True)
             self.delete_me = True
             self.error_messages = f"Error while parsing DSSR's json output:\n{e}"
             return 1
+        # Remove nucleotides of the chain that are outside the Rfam mapping, if any
+        if self.pdb_start and self.pdb_end:
+            df = df.drop(df[(df.nt_resnum < self.pdb_start) | (df.nt_resnum > self.pdb_end)].index)
+
+        #############################################
+        # Solve some common issues and drop ligands
+        #############################################
+
         # Shift numbering when duplicate residue numbers are found.
         # Example: 4v9q-DV contains 17 and 17A which are both read 17 by DSSR.
         while True in df.duplicated(['nt_resnum']).values:
@@ -266,11 +295,12 @@ class Chain:
             or (len(df.index_chain) >= 2 and df.iloc[[-1]].nt_resnum.iloc[0] > 50 + df.iloc[[-2]].nt_resnum.iloc[0])):
             df = df.head(-1) 
-        # Assert some nucleotides exist
+        # Assert some nucleotides still exist
         try:
             l = df.iloc[-1,1] - df.iloc[0,1] + 1    # length of chain from nt_resnum point of view
         except IndexError:
             warn(f"Could not find real nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. Ignoring chain {self.chain_label}.", error=True)
+            no_nts_set.add(self.pdb_id)
             self.delete_me = True
             self.error_messages = f"Could not find nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. We expect a problem with {self.pdb_id} mmCIF download. Delete it and retry."
             return 1
@@ -281,14 +311,15 @@ class Chain:
             df.iloc[:, 0] -= st
         df = df.drop(df[df.index_chain < 0].index)  # drop eventual ones with index_chain < the first residue (usually, ligands)
-        # Add a sequence column just for the alignments
+        # Re-Assert some nucleotides still exist
-        df['nt_align_code'] = [ str(x).upper()
+        try:
-                                        .replace('NAN', '-')      # Unresolved nucleotides are gaps
+            l = df.iloc[-1,1] - df.iloc[0,1] + 1    # length of chain from nt_resnum point of view
-                                        .replace('?', '-')        # Unidentified residues, let's delete them
+        except IndexError:
-                                        .replace('T', 'U')        # 5MU are modified to t, which gives T
+            warn(f"Could not find real nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. Ignoring chain {self.chain_label}.", error=True)
-                                        .replace('P', 'U')        # Pseudo-uridines, but it is not really right to change them to U, see DSSR paper, Fig 2
+            no_nts_set.add(self.pdb_id)
-                                for x in df['nt_code'] ]
+            self.delete_me = True
-
+            self.error_messages = f"Could not find nucleotides of chain {self.pdb_chain_id} in annotation {self.pdb_id}.json. We expect a problem with {self.pdb_id} mmCIF download. Delete it and retry."
+            return 1
         # Add eventual missing rows because of unsolved residues in the chain.
         # Sometimes, the 3D structure is REALLY shorter than the family it's mapped to,
@@ -304,7 +335,7 @@ class Chain:
         # portion solved in 3D  1 |--------------|79 85|------------| 156
         # Rfam mapping           3 |------------------------------------------ ... -------| 3353 (yes larger, 'cause it could be inferred)
         # nt resnum              3 |--------------------------------|  156
-        # index_chain            1 |-------------|77 83|------------|  149
+        # index_chain            1 |-------------|77 83|------------|  149 (before correction)
         # expected data point    1 |--------------------------------|  154
         #
@@ -314,15 +345,26 @@ class Chain:
             for i in sorted(diff):
                 # Add a row at position i
                 df = pd.concat([    df.iloc[:i], 
-                                    pd.DataFrame({"index_chain": i+1, "nt_resnum": i+resnum_start, 
+                                    pd.DataFrame({"index_chain": i+1, "nt_resnum": i+resnum_start, "nt_code":'-', "nt_name":'-'}, index=[i]), 
-                                                    "nt_code":'-', "nt_name":'-', 'nt_align_code':'-'}, index=[i]), 
+                                    df.iloc[i:]       ])
-                                    df.iloc[i:]
-                                ])
                 # Increase the index_chain of all following lines
                 df.iloc[i+1:, 0] += 1
             df = df.reset_index(drop=True)
         self.full_length = len(df.index_chain)
+        # Add a sequence column just for the alignments
+        df['nt_align_code'] = [ str(x).upper()
+                                    .replace('NAN', '-')      # Unresolved nucleotides are gaps
+                                    .replace('?', '-')        # Unidentified residues, let's delete them
+                                    .replace('T', 'U')        # 5MU are modified to t, which gives T
+                                    .replace('P', 'U')        # Pseudo-uridines, but it is not really right to change them to U, see DSSR paper, Fig 2
+                                for x in df['nt_code'] ]
+
+
+        #######################################
+        # Compute new features
+        #######################################
+
         # One-hot encoding sequence
         df["is_A"] = [ 1 if x=="A" else 0 for x in df["nt_code"] ]
         df["is_C"] = [ 1 if x=="C" else 0 for x in df["nt_code"] ]
@@ -415,7 +457,13 @@ class Chain:
                         newpairs.append(v)
             df['paired'] = newpairs
-        # Saving to database
+        self.seq = "".join(df.nt_code)
+        self.seq_to_align = "".join(df.nt_align_code)
+        self.length = len([ x for x in self.seq_to_align if x != "-" ])
+
+        ####################################
+        # Save everything to database
+        ####################################
         with sqlite3.connect(runDir+"/results/RNANet.db", timeout=10.0) as conn:
             # Register the chain in table chain
             if self.pdb_start is not None:
@@ -444,11 +492,7 @@ class Chain:
                 ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?);""", 
             many=True, data=list(df.to_records(index=False)), warn_every=10)
-        # Now load data from the database
+        
-        self.seq = "".join(df.nt_code)
-        self.seq_to_align = "".join(df.nt_align_code)
-        self.length = len([ x for x in self.seq_to_align if x != "-" ])
-
         # Remove too short chains
         if self.length < 5:
             warn(f"{self.chain_label} sequence is too short, let's ignore it.\t", error=True)
@@ -811,7 +855,7 @@ class Pipeline:
         global path_to_seq_data
         try:
-            opts, args = getopt.getopt( sys.argv[1:], "r:hs", 
+            opts, _ = getopt.getopt( sys.argv[1:], "r:hs", 
                                     [   "help", "resolution=", "keep-hetatm=", "from-scratch",
                                         "fill-gaps=", "3d-folder=", "seq-folder=",
                                         "no-homology", "ignore-issues", "extract",
@@ -822,7 +866,7 @@ class Pipeline:
         for opt, arg in opts:
-            if opt in ["--from-scratch", "--update-mmcifs", "--update-homolgous"] and "tobedefinedbyoptions" in [path_to_3D_data, path_to_seq_data]:
+            if opt in ["--from-scratch", "--update-mmcifs", "--update-homologous"] and "tobedefinedbyoptions" in [path_to_3D_data, path_to_seq_data]:
                 print("Please provide --3d-folder and --seq-folder first, so that we know what to delete and update.")
                 exit()
@@ -890,7 +934,7 @@ class Pipeline:
                 warn("Deleting previous database and recomputing from scratch.")
                 subprocess.run(["rm", "-rf", 
                                 path_to_3D_data + "annotations",
-                                path_to_3D_data + "RNAcifs",
+                                # path_to_3D_data + "RNAcifs",  # DEBUG : keep the cifs !
                                 path_to_3D_data + "rna_mapped_to_Rfam",
                                 path_to_3D_data + "rnaonly",
                                 path_to_seq_data + "realigned",
@@ -899,9 +943,6 @@ class Pipeline:
                                 runDir + "/known_issues_reasons.txt", 
                                 runDir + "/results/RNANet.db"])
             elif opt == "--update-homologous":
-                if "tobedefinedbyoptions" == path_to_seq_data:
-                    warn("Please provide --seq-folder before --update-homologous in the list of options.", error=True)
-                    exit(1)
                 warn("Deleting previous sequence files and recomputing alignments.")
                 subprocess.run(["rm", "-rf", 
                                 path_to_seq_data + "realigned",
@@ -942,8 +983,9 @@ class Pipeline:
             print("> Building list of structures...", flush=True)
             p = Pool(initializer=init_worker, initargs=(tqdm.get_lock(),), processes=ncores)
             try:
+
                 pbar = tqdm(full_structures_list, maxinterval=1.0, miniters=1, bar_format="{percentage:3.0f}%|{bar}|")
-                for i, newchains in enumerate(p.imap_unordered(partial(work_infer_mappings, not self.REUSE_ALL, allmappings), full_structures_list)): 
+                for _, newchains in enumerate(p.imap_unordered(partial(work_infer_mappings, not self.REUSE_ALL, allmappings), full_structures_list)): 
                     self.update += newchains
                     pbar.update(1) # Everytime the iteration finishes, update the global progress bar
@@ -1053,7 +1095,8 @@ class Pipeline:
                     warn(f"Adding {c[1].chain_label} to known issues.")
                     ki.write(c[1].chain_label + '\n')
                     kir.write(c[1].chain_label + '\n' + c[1].error_messages + '\n\n')
-                    sql_execute(conn, f"UPDATE chain SET issue = 1 WHERE chain_id = ?;", data=(c[1].db_chain_id,))
+                    with sqlite3.connect(runDir+"/results/RNANet.db") as conn:
+                        sql_execute(conn, f"UPDATE chain SET issue = 1 WHERE chain_id = ?;", data=(c[1].db_chain_id,))
         ki.close()
         kir.close()
@@ -1194,7 +1237,7 @@ class Pipeline:
         p = Pool(initializer=init_worker, initargs=(tqdm.get_lock(),), processes=3)
         try:
             pbar = tqdm(total=len(self.loaded_chains), desc="Saving chains to CSV", position=0, leave=True) 
-            for i, _ in enumerate(p.imap_unordered(work_save, self.loaded_chains)):
+            for _, _2 in enumerate(p.imap_unordered(work_save, self.loaded_chains)):
                 pbar.update(1) 
             pbar.close()
             p.close()
@@ -2158,6 +2201,8 @@ if __name__ == "__main__":
         pp.dl_and_annotate(retry=True, coeff_ncores=0.3)  #
         pp.build_chains(retry=True, coeff_ncores=1.0)     # Use half the cores to reduce required amount of memory
     print(f"> Loaded {len(pp.loaded_chains)} RNA chains ({len(pp.update) - len(pp.loaded_chains)} errors).")
+    if len(no_nts_set):
+        print(f"Among errors, {len(no_nts_set)} structures seem to contain RNA chains without defined nucleotides:", no_nts_set, flush=True)
     pp.checkpoint_save_chains()
     if not pp.HOMOLOGY:
@@ -2165,7 +2210,7 @@ if __name__ == "__main__":
         for c in pp.loaded_chains:
             work_save(c, homology=False)
         print("Completed.")
-        exit()
+        exit(0)
     # At this point, structure, chain and nucleotide tables of the database are up to date.
     # (Modulo some statistics computed by statistics.py)
--- a/statistics.py
View file @c21f55f
+++ b/statistics.py
View file @c21f55f
@@ -168,6 +168,8 @@ def stats_len():
     lengths = []
     conn = sqlite3.connect("results/RNANet.db")
     for i,f in enumerate(fam_list):
+
+        # Define a color for that family in the plot
         if f in LSU_set:
             cols.append("red") # LSU
         elif f in SSU_set:
@@ -178,11 +180,15 @@ def stats_len():
             cols.append("orange")
         else:
             cols.append("grey")
+
+        # Get the lengths of chains
         l = [ x[0] for x in sql_ask_database(conn, f"SELECT COUNT(index_chain) FROM (SELECT chain_id FROM chain WHERE rfam_acc='{f}') NATURAL JOIN nucleotide GROUP BY chain_id;") ]
         lengths.append(l)
+
         notify(f"[{i+1}/{len(fam_list)}] Computed {f} chains lengths")
     conn.close()
+    # Plot the figure
     fig = plt.figure(figsize=(10,3))
     ax = fig.gca()
     ax.hist(lengths, bins=100, stacked=True, log=True, color=cols, label=fam_list)
@@ -191,6 +197,8 @@ def stats_len():
     ax.set_xlim(left=-150)
     ax.tick_params(axis='both', which='both', labelsize=8)
     fig.tight_layout()
+
+    # Draw the legend
     fig.subplots_adjust(right=0.78)
     filtered_handles = [mpatches.Patch(color='red'), mpatches.Patch(color='white'), mpatches.Patch(color='white'), mpatches.Patch(color='white'),
                         mpatches.Patch(color='blue'), mpatches.Patch(color='white'), mpatches.Patch(color='white'),
@@ -204,6 +212,8 @@ def stats_len():
                        'Other']
     ax.legend(filtered_handles, filtered_labels, loc='right', 
                 ncol=1, fontsize='small', bbox_to_anchor=(1.3, 0.5))
+
+    # Save the figure
     fig.savefig("results/figures/lengths.png")
     notify("Computed sequence length statistics and saved the figure.")
@@ -224,10 +234,12 @@ def stats_freq():
     Outputs results/frequencies.csv
     REQUIRES tables chain, nucleotide up to date."""
+    # Initialize a Counter object for each family
     freqs = {}
     for f in fam_list:
         freqs[f] = Counter()
+    # List all nt_names happening within a RNA family and store the counts in the Counter
     conn = sqlite3.connect("results/RNANet.db")
     for i,f in enumerate(fam_list):
         counts = dict(sql_ask_database(conn, f"SELECT nt_name, COUNT(nt_name) FROM (SELECT chain_id from chain WHERE rfam_acc='{f}') NATURAL JOIN nucleotide GROUP BY nt_name;"))
@@ -235,6 +247,7 @@ def stats_freq():
         notify(f"[{i+1}/{len(fam_list)}] Computed {f} nucleotide frequencies.")
     conn.close()
+    # Create a pandas DataFrame, and save it to CSV.
     df = pd.DataFrame()
     for f in fam_list:
         tot = sum(freqs[f].values())
@@ -347,8 +360,8 @@ def stats_pairs():
             fam_pbar = tqdm(total=len(fam_list), desc="Pair-types in families", position=0, leave=True) 
             results = []
             allpairs = []
-            for i, _ in enumerate(p.imap_unordered(parallel_stats_pairs, fam_list)):
+            for _, newp_famdf in enumerate(p.imap_unordered(parallel_stats_pairs, fam_list)):
-                newpairs, fam_df = _
+                newpairs, fam_df = newp_famdf
                 fam_pbar.update(1)
                 results.append(fam_df)
                 allpairs.append(newpairs)
@@ -432,13 +445,14 @@ def seq_idty():
     Creates temporary results files in data/*.npy
     REQUIRES tables chain, family un to date."""
+    # List the families for which we will compute sequence identity matrices
     conn = sqlite3.connect("results/RNANet.db")
     famlist = [ x[0] for x in sql_ask_database(conn, "SELECT rfam_acc from (SELECT rfam_acc, COUNT(chain_id) as n_chains FROM family NATURAL JOIN chain GROUP BY rfam_acc) WHERE n_chains > 1 ORDER BY rfam_acc ASC;") ]
     ignored = [ x[0] for x in sql_ask_database(conn, "SELECT rfam_acc from (SELECT rfam_acc, COUNT(chain_id) as n_chains FROM family NATURAL JOIN chain GROUP BY rfam_acc) WHERE n_chains < 2 ORDER BY rfam_acc ASC;") ]
     if len(ignored):
         print("Idty matrices: Ignoring families with only one chain:", " ".join(ignored)+'\n')
-    # compute distance matrices
+    # compute distance matrices (or ignore if data/RF0****.npy exists)
     p = Pool(processes=8)
     p.map(to_dist_matrix, famlist)
     p.close()
@@ -502,7 +516,7 @@ def per_chain_stats():
         # Set the values
         sql_execute(conn, "UPDATE chain SET chain_freq_A = ?, chain_freq_C = ?, chain_freq_G = ?, chain_freq_U = ?, chain_freq_other = ? WHERE chain_id= ?;",
-                       many=True, data=list(df.to_records(index=False)), warn_every=10)
+                          many=True, data=list(df.to_records(index=False)), warn_every=10)
     notify("Updated the database with per-chain base frequencies")
 if __name__ == "__main__":