Repository: vrmarcelino/CCMetagen (language: python)

Issue:
Tidy up repository
This looks like a big PR, but it's mostly housekeeping to improve repository readability.

**Summary of changes:**
- Remove unnecessary `.Rhistory` files
- Move `benchmarking` and `tutorial` directories to `docs` directory
- Update README.md
  - Use conda installation instructions as recommended
  - Back up old README file in `docs` directory
  - Update links to reflect new repository structure
- Rename `contributions.txt` file
- Apply Black formatting


Discussion / hints:

Fix patch:
diff --git a/.Rhistory b/.Rhistory
deleted file mode 100644
index e69de29..0000000
diff --git a/Contributions.txt b/AUTHORS.md
similarity index 73%
rename from Contributions.txt
rename to AUTHORS.md
index d201bd6..1bf580c 100644
--- a/Contributions.txt
+++ b/AUTHORS.md
@@ -4,6 +4,7 @@ CCMetagen code:
 Vanessa R. Marcelino
 Jan P. Buchmann
 Andreas Sjödin
+Vini Salazar (maintenance, database development)
 
 CCMetagen web server and KMA:
 Philip T.L.C. Clausen
diff --git a/CCMetagen.py b/CCMetagen.py
index 7fa956f..4ed79f3 100755
--- a/CCMetagen.py
+++ b/CCMetagen.py
@@ -1,27 +1,24 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+
 """
 CCMetagen main script
 
 @ V.R.Marcelino
-Created on Wed Jul 25 17:13:10 2018
-
+Created on Wed Jul 25 17:13:10 2018.
 """
 
 # imports
 import sys
 import pandas as pd
-from argparse import ArgumentParser
+import argparse
 import subprocess
 import os
 from ete3 import NCBITaxa
 import re
 
 # local imports
-from ccmetagen import __version__
-from ccmetagen import fParseKMA
-from ccmetagen import cTaxInfo # needed in fParseKMA
-from ccmetagen import fNCBItax # needed in fParseKMA
+from ccmetagen import __version__, _ARGPARSE_DEFAULTS, fParseKMA
 
 # version formatting
 version_numb = "v" + __version__
@@ -29,97 +26,209 @@
 # help
 if __name__ == "__main__":
     if len(sys.argv) == 1:
-        print ("")
-        print ("CCMetagen - Identify species in metagenome datasets")
-        print (version_numb)
-        print ("To be used with KMA")
-        print ("")
-        print ("Usage: CCMetagen.py <options> ")
-        print ("Ex: CCMetagen.py -i KMA_out/2_mtg.res -o 2_mtg_result")
-        print ("")
-        print ("")
-        print ("""When running CCMetagen on multiple files in a folder:
+        print("")
+        print("CCMetagen - Identify species in metagenome datasets")
+        print(version_numb)
+        print("To be used with KMA")
+        print("")
+        print("Usage: CCMetagen.py <options> ")
+        print("Ex: CCMetagen.py -i KMA_out/2_mtg.res -o 2_mtg_result")
+        print("")
+        print("")
+        print(
+            """When running CCMetagen on multiple files in a folder:
     input_dir=KMA_out
     output_dir=CCMetagen_results
     mkdir $output_dir
     for f in $input_dir/*.res; do
         out=$output_dir/${f/$input_dir\\/}
         CCMetagen.py -i $f -o $out
-    done""")
-        print ("")
-        print ("For help and options, type: CCMetagen.py -h")
-        print ("")
-        print ("")
+    done"""
+        )
+        print("")
+        print("For help and options, type: CCMetagen.py -h")
+        print("")
+        print("")
         sys.exit()
 
-
-    parser = ArgumentParser()
-
-    parser.add_argument('-m', '--mode', default = 'both',
-                        help="""what do you want CCMetagen to do?
-                        Valid options are 'visual', 'text' or 'both':
-                            text: parses kma, filters based on quality and output a text file with taxonomic information and detailed mapping information
-                            visual: parses kma, filters based on quality and output a simplified text file and a krona html file for visualization
-                            both: outputs both text and visual file formats. Default = both""", required=False)
-
-    parser.add_argument('-i', '--res_fp', help='Path to the KMA result (.res file)', required=True)
-    parser.add_argument('-o', '--output_fp', default = 'CCMetagen_out',
-                        help='Path to the output file. Default = CCMetagen_out', required=False)
-    parser.add_argument('-r', '--reference_database', default = 'nt',
-                        help='Which reference database was used. Options: UNITE, RefSeq or nt. Default = nt', required=False)
-    parser.add_argument('-ef', '--extended_output_file', default = 'n',
-                        help="""Produce an extended output file that includes the percentage of classified reads.
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+
+    parser.add_argument(
+        "-m",
+        "--mode",
+        default=_ARGPARSE_DEFAULTS.get("mode"),
+        help="""
+        What would you like CCMetagen to do?\n
+        Valid options are 'visual', 'text' or 'both':\n
+            \ttext: parses kma, filters based on quality and output a text file with taxonomic information and detailed mapping information.\n
+            \tvisual: parses kma, filters based on quality and output a simplified text file and a krona html file for visualization.\n
+            \tboth: outputs both text and visual file formats.
+            """,
+        required=False,
+    )
+
+    parser.add_argument(
+        "-i", "--res_fp", help="Path to the KMA result (.res file)", required=True
+    )
+    parser.add_argument(
+        "-o",
+        "--output_fp",
+        default=_ARGPARSE_DEFAULTS.get("output_fp"),
+        help="Path to the output file.",
+        required=False,
+    )
+    parser.add_argument(
+        "-r",
+        "--reference_database",
+        default=_ARGPARSE_DEFAULTS.get("reference_database"),
+        help="Which reference database was used. Options: UNITE, RefSeq or nt.",
+        required=False,
+    )
+    parser.add_argument(
+        "-ef",
+        "--extended_output_file",
+        default=_ARGPARSE_DEFAULTS.get("extended_output_file"),
+        help="""Produce an extended output file that includes the percentage of classified reads.
                         Options: y or n. To use this featire, you need to generate the mapstat file when
                         required unning KMA (use flag -ef), and use it as input in CCMetagen (flag --mapstat). 
-                        Default = n""", required=False)
-
-    parser.add_argument('-du', '--depth_unit', default = 'kma',
-                        help="""Desired unit for Depth(abundance) measurements.
+                        """,
+        required=False,
+    )
+
+    parser.add_argument(
+        "-du",
+        "--depth_unit",
+        default=_ARGPARSE_DEFAULTS.get("depth_unit"),
+        help="""Desired unit for Depth(abundance) measurements.
                         Default = kma (KMA default depth, which is the number of nucleotides overlapping each template,
                         divided by the lengh of the template).
                         Alternatively, you can have abundance calculated in Reads Per Million (RPM, option 'rpm'), 
                         in number of nucleotides overlaping the template (option 'nc') or in number of fragments (i.e. PE reads, option 'fr').
                         If you use the 'nc', 'rpm' or 'fr' options, remember to change the default --depth parameter accordingly.
-                        Valid options are nc, rpm, fr and kma""", required=False)
-    parser.add_argument('-map', '--mapstat', help="""Path to the mapstat file produced with KMA when using the -ef flag (.mapstat).
+                        Valid options are 'nc', 'rpm', 'fr' and 'kma'.""",
+        required=False,
+    )
+    parser.add_argument(
+        "-map",
+        "--mapstat",
+        default=_ARGPARSE_DEFAULTS.get("mapstat"),
+        help="""Path to the mapstat file produced with KMA when using the -ef flag (.mapstat).
                         Required when calculating abundances in RPM or in number of fragments, visualing the abundances in readCounts or readContAln
-                        for the krona graph or when producing the extended_output_file""", required = False)
-    parser.add_argument('-d', '--depth', default = 0.2,
-                        help="""Minimum sequencing depth. Default = 0.2. The unit corresponds to the one used with --depth_unit
+                        for the krona graph or when producing the extended_output_file.""",
+        required=False,
+    )
+    parser.add_argument(
+        "-d",
+        "--depth",
+        default=_ARGPARSE_DEFAULTS.get("depth"),
+        help="""Minimum sequencing depth. The unit corresponds to the one used with --depth_unit
                         If you use --depth_unit different from the default, change this accordingly.
-                        """,type=float, required=False)
-    parser.add_argument('-c', '--coverage', default = 20,
-                        help="""Minimum coverage. Default = 20 (i.e. 20%% of the reference sequence)""",type=float, required=False)
-    parser.add_argument('-q', '--query_identity', default = 50,
-                        help='Minimum query identity (Phylum level). Default = 50', type=float, required=False)
-    parser.add_argument('-p', '--pvalue', default = 0.05,
-                        help='Minimum p-value. Default = 0.05.',type=float, required=False)
-    parser.add_argument('-k', '--krona_mode', default = 'Depth',
-                        help="""Abundance measure for the Krona graph. Default = Depth. 
-                        Alternatively, you can choose readCount ('rc') or readCountAln ('rca').
+                        """,
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-c",
+        "--coverage",
+        default=_ARGPARSE_DEFAULTS.get("coverage"),
+        help="""Percetange of minimum coverage, e.g. 20 means 20%% of the reference sequence.""",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-q",
+        "--query_identity",
+        default=_ARGPARSE_DEFAULTS.get("query_identity"),
+        help="Minimum query identity (Phylum level).",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-p",
+        "--pvalue",
+        default=_ARGPARSE_DEFAULTS.get("pvalue"),
+        help="Minimum p-value.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-k",
+        "--krona_mode",
+        default=_ARGPARSE_DEFAULTS.get("krona_mode"),
+        help="""Abundance measure for the Krona graph. 
+                        You can choose depth ('Depth'), readCount ('rc') or readCountAln ('rca').
                         'rc' and 'rca' require a specified --mapstat file.
-                        """, required=False)                 
-    parser.add_argument('-tax', '--local_taxfile', default = None,
-                        help='Use if a local taxdump file wants to be used. Default = None', required=False)
-
+                        """,
+        required=False,
+    )
+    parser.add_argument(
+        "-tax",
+        "--local_taxfile",
+        default=_ARGPARSE_DEFAULTS.get("local_taxfile"),
+        help="Use if a local taxdump file wants to be used.",
+        required=False,
+    )
 
     # similarity thresholds:
-    parser.add_argument('-st', '--species_threshold', default = 98.41,
-                        help='Species-level similarity threshold. Default = 98.41',type=float, required=False)
-    parser.add_argument('-gt', '--genus_threshold', default = 96.31,
-                        help='Genus-level similarity threshold. Default = 96.31',type=float, required=False)
-    parser.add_argument('-ft', '--family_threshold', default = 88.51,
-                        help='Family-level similarity threshold. Default = 88.51',type=float, required=False)
-    parser.add_argument('-ot', '--order_threshold', default = 81.21,
-                        help='Order-level similarity threshold. Default = 81.21',type=float, required=False)
-    parser.add_argument('-ct', '--class_threshold', default = 80.91,
-                        help='Class-level similarity threshold. Default = 80.91',type=float, required=False)
-    parser.add_argument('-pt', '--phylum_threshold', default = 0,
-                        help='Phylum-level similarity threshold. Default = 0 - not applied',type=float, required=False)
-    parser.add_argument('-off', '--turn_off_sim_thresholds', default = 'n',
-                        help='Turns simularity-based filtering off. Options = y or n. Default = n', required=False)
-
-    parser.add_argument('-v', '--version', action='version', version=version_numb)
+    parser.add_argument(
+        "-st",
+        "--species_threshold",
+        default=_ARGPARSE_DEFAULTS.get("species_threshold"),
+        help="Species-level similarity threshold.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-gt",
+        "--genus_threshold",
+        default=_ARGPARSE_DEFAULTS.get("genus_threshold"),
+        help="Genus-level similarity threshold.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-ft",
+        "--family_threshold",
+        default=_ARGPARSE_DEFAULTS.get("family_threshold"),
+        help="Family-level similarity threshold.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-ot",
+        "--order_threshold",
+        default=_ARGPARSE_DEFAULTS.get("order_threshold"),
+        help="Order-level similarity threshold.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-ct",
+        "--class_threshold",
+        default=_ARGPARSE_DEFAULTS.get("class_threshold"),
+        help="Class-level similarity threshold.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-pt",
+        "--phylum_threshold",
+        default=_ARGPARSE_DEFAULTS.get("phylum_threshold"),
+        help="Phylum-level similarity threshold. 0 means not applied.",
+        type=float,
+        required=False,
+    )
+    parser.add_argument(
+        "-off",
+        "--turn_off_sim_thresholds",
+        default=_ARGPARSE_DEFAULTS.get("turn_off_sim_thresholds"),
+        help="Turns simularity-based filtering off. Options = 'y' or 'n'.",
+        required=False,
+    )
+
+    parser.add_argument("-v", "--version", action="version", version=version_numb)
 
     args = parser.parse_args()
 
@@ -138,11 +247,10 @@
     if taxfile is not None and not os.path.isfile(taxfile):
         raise ValueError("If used, the argument taxfile has to point to a taxdump file")
 
-
     # taxonomic thresholds:
     off = args.turn_off_sim_thresholds
 
-    if off == 'y':
+    if off == "y":
         st = 0
         gt = 0
         ft = 0
@@ -150,7 +258,7 @@
         ct = 0
         pt = 0
 
-    elif off == 'n':
+    elif off == "n":
         st = args.species_threshold
         gt = args.genus_threshold
         ft = args.family_threshold
@@ -159,188 +267,240 @@
         pt = args.phylum_threshold
 
     else:
-        print ("-off argument should be either y or n. ")
+        print("-off argument should be either y or n. ")
         sys.exit("Try again.")
 
     # developing and debugging:
-    #args.output_fp = "CCMetagen_nt_results"
-    #f = "sim_metagen.res"
-    #mapstat="sim_metagen.mapstat"
-    #ref_database = "nt"
-    #mode = 'text'
-    #c = 20
-    #q = 50
-    #d = 0.2
-    #p = 0.05
-    #st = 99
-    #gt = 98
-    #ft = 95
-    #ot = 80
-    #ct = 0
-    #pt = 0
-    #du = 'kma'
-    #ef = 'y' # extended format - add a flag, default = 'n'
-    #k = 'Depth'
-
+    # args.output_fp = "CCMetagen_nt_results"
+    # f = "sim_metagen.res"
+    # mapstat="sim_metagen.mapstat"
+    # ref_database = "nt"
+    # mode = 'text'
+    # c = 20
+    # q = 50
+    # d = 0.2
+    # p = 0.05
+    # st = 99
+    # gt = 98
+    # ft = 95
+    # ot = 80
+    # ct = 0
+    # pt = 0
+    # du = 'kma'
+    # ef = 'y' # extended format - add a flag, default = 'n'
+    # k = 'Depth'
 
     ##### Checks:
 
     # Run implicitly ete3.NCBITaxa.__init__() to check for valid taxonomy database
-    #NCBITaxa()
+    # NCBITaxa()
     if taxfile is not None:
         NCBITaxa(taxfile)
     else:
         NCBITaxa()
-        
 
     # Warning if RefDatabase is unknown
-    if ref_database not in ("UNITE", "RefSeq","nt"):
-        print (""" Reference database (-r) must be either UNITE, RefSeq or nt.
-            the input is case sensitive and the default is nt.""")
+    if ref_database not in ("UNITE", "RefSeq", "nt"):
+        print(
+            """ Reference database (-r) must be either UNITE, RefSeq or nt.
+            the input is case sensitive and the default is nt."""
+        )
         sys.exit("Try again.")
 
-
     # check if ef flag is correct
     if ef not in ("y", "n"):
-        print ("Unrecognized argument %s. Use either '-ef n' (default) or '-ef y'" %(ef))
+        print(
+            "Unrecognized argument %s. Use either '-ef n' (default) or '-ef y'" % (ef)
+        )
         sys.exit("Try again.")
 
-
     ##### Read input files and output a pandas dataframe
-    print ("")
-    print ("Reading file %s" %(f))
-    print ("")
+    print("")
+    print("Reading file %s" % (f))
+    print("")
 
-    df = pd.read_csv(f, sep='\t', index_col=0, encoding='latin1')
+    df = pd.read_csv(f, sep="\t", index_col=0, encoding="latin1")
 
     # Rename headers:
     df.index.name = "Closest_match"
 
-
     ##### Adjust depth to reflect number of bases or RPM if needed:
 
     # number of nucleotides:
-    if du == 'nc':
-        df['Depth'] = df.Depth * df.Template_length
-        print ("Calculating depth as number of nucleotides, ignoring template length.")
-        print ("""Remember to adjust minimum depth value (ex: -d 200) to filter low abundance hits.""")
-
-    # RPM:   
-    elif du == 'rpm':
-        print ("Calculating RPM...")
-        print ("""
+    if du == "nc":
+        df["Depth"] = df.Depth * df.Template_length
+        print("Calculating depth as number of nucleotides, ignoring template length.")
+        print(
+            """Remember to adjust minimum depth value (ex: -d 200) to filter low abundance hits."""
+        )
+
+    # RPM:
+    elif du == "rpm":
+        print("Calculating RPM...")
+        print(
+            """
             Note 1: to calculate RPM, you need to generate the mapstat file when
             running KMA (flag -ef), and use it as input in CCMetagen (flag --mapstat).
             
             Note 2: you might want to adjust the minimum depth (-d) value accordingly.
             The default minimum depth is 0.2.
-            """)
+            """
+        )
 
         with open(mapstat) as mapfile:
-            fragments_line=mapfile.readlines()[3]
-        total_frags = re.split(r'(\t|\n)',fragments_line)[2]
-        df_stats = pd.read_csv(mapstat, sep='\t', index_col=0, header = 6, encoding='latin1')
-        df['Depth'] = 1000000 * df_stats['fragmentCount'] / int(total_frags)
-
-
-    # number of PE reads (frags):   
-    elif du == 'fr':
-        print ("Calculating number of PE reads (fragments)")
-        print ("""
+            fragments_line = mapfile.readlines()[3]
+        total_frags = re.split(r"(\t|\n)", fragments_line)[2]
+        df_stats = pd.read_csv(
+            mapstat, sep="\t", index_col=0, header=6, encoding="latin1"
+        )
+        df["Depth"] = 1000000 * df_stats["fragmentCount"] / int(total_frags)
+
+    # number of PE reads (frags):
+    elif du == "fr":
+        print("Calculating number of PE reads (fragments)")
+        print(
+            """
             Note 1: to calculate number of fragments, you need to generate the mapstat file when
             running KMA (flag -ef), and use it as input in CCMetagen (flag --mapstat).
             
             Note 2: you might want to adjust the minimum depth (-d) value accordingly.
             Ex: Use -d 2 to only consider matches with 2 fragments or more.
-            """)
+            """
+        )
 
         with open(mapstat) as mapfile:
-            fragments_line=mapfile.readlines()[3]
-        total_frags = re.split(r'(\t|\n)',fragments_line)[2]
-        df_stats = pd.read_csv(mapstat, sep='\t', index_col=0, header = 6, encoding='latin1')
-        df['Depth'] = df_stats['fragmentCount']
+            fragments_line = mapfile.readlines()[3]
+        total_frags = re.split(r"(\t|\n)", fragments_line)[2]
+        df_stats = pd.read_csv(
+            mapstat, sep="\t", index_col=0, header=6, encoding="latin1"
+        )
+        df["Depth"] = df_stats["fragmentCount"]
 
-    elif du == 'kma':
-        print ("")
+    elif du == "kma":
+        print("")
 
     else:
-        print ("""Warning: the depth unit you specified makes no sense.
-            --depth_unit option must be nc, rpm, fr or kma. Using 'kma'.""")
-        print ("")
+        print(
+            """Warning: the depth unit you specified makes no sense.
+            --depth_unit option must be nc, rpm, fr or kma. Using 'kma'."""
+        )
+        print("")
 
     ##### Quality control + taxonomic assignments
-        
+
     # quality filter (coverage, query identity, Depth and p-value)
     df = fParseKMA.res_filter(df, ref_database, c, q, d, p)
 
     # add tax info
     df = fParseKMA.populate_w_tax(df, ref_database, st, gt, ft, ot, ct, pt)
 
-
     ##### Output a file with tax info
-    if (mode == 'text') or (mode == 'both'):
-
+    if (mode == "text") or (mode == "both"):
         # save to file
         out = args.output_fp + ".ccm.csv"
         pd.DataFrame.to_csv(df, out)
 
-        print ("csv file saved as %s" %(out))
-        print ("")
+        print("csv file saved as %s" % (out))
+        print("")
 
     ##### Output a Krona file based on specified mode
-    if (mode == 'visual') or (mode == 'both'):
-        
+    if (mode == "visual") or (mode == "both"):
         # ReadCount (number of reads mapped to the template)
-        if (k == 'rc'):
+        if k == "rc":
             if mapstat is not None:
-                print ("Extracting read count from mapstat file")
+                print("Extracting read count from mapstat file")
                 # load mapstat file if needed
-                if du == 'kma' or du == 'nc':
+                if du == "kma" or du == "nc":
                     print(mapstat)
-                    df_stats = pd.read_csv(mapstat, sep='\t', index_col=0, header = 6, encoding='latin1')
-                    
-                df['readCount'] = df_stats['readCount']
-                krona_info = df[['readCount','Superkingdom','Kingdom','Phylum','Class','Order','Family','Genus','Species']]
+                    df_stats = pd.read_csv(
+                        mapstat, sep="\t", index_col=0, header=6, encoding="latin1"
+                    )
+
+                df["readCount"] = df_stats["readCount"]
+                krona_info = df[
+                    [
+                        "readCount",
+                        "Superkingdom",
+                        "Kingdom",
+                        "Phylum",
+                        "Class",
+                        "Order",
+                        "Family",
+                        "Genus",
+                        "Species",
+                    ]
+                ]
             else:
                 print("""Error, no mapstat file specified""")
-                print ("""To calculate the read count for this krona mode, you need to generate the mapstat file when
+                print(
+                    """To calculate the read count for this krona mode, you need to generate the mapstat file when
                 running KMA (flag -ef), and use it as input in CCMetagen (flag --mapstat).
-                """)
-                print ("")
-                
+                """
+                )
+                print("")
+
         # ReadCountAln (number of reads aligned)
-        elif (k == 'rca'):
+        elif k == "rca":
             if mapstat is not None:
-                print ("Extracting read count alignment from mapstat file")
+                print("Extracting read count alignment from mapstat file")
                 # load mapstat file if needed
-                if du == 'kma' or du == 'nc':
-                    df_stats = pd.read_csv(mapstat, sep='\t', index_col=0, header = 6, encoding='latin1')
-                
-                df['readCountAln'] = df_stats['readCountAln']
-                krona_info = df[['readCountAln','Superkingdom','Kingdom','Phylum','Class','Order','Family','Genus','Species']]
+                if du == "kma" or du == "nc":
+                    df_stats = pd.read_csv(
+                        mapstat, sep="\t", index_col=0, header=6, encoding="latin1"
+                    )
+
+                df["readCountAln"] = df_stats["readCountAln"]
+                krona_info = df[
+                    [
+                        "readCountAln",
+                        "Superkingdom",
+                        "Kingdom",
+                        "Phylum",
+                        "Class",
+                        "Order",
+                        "Family",
+                        "Genus",
+                        "Species",
+                    ]
+                ]
             else:
                 print("""Error, no mapstat file specified""")
-                print ("""To calculate the read count alignment for this krona mode, you need to generate the mapstat file when
+                print(
+                    """To calculate the read count alignment for this krona mode, you need to generate the mapstat file when
                     running KMA (flag -ef), and use it as input in CCMetagen (flag --mapstat).
-                    """)
-                print ("")
+                    """
+                )
+                print("")
 
         # Depth (default)
-        elif k == 'Depth':
-            krona_info = df[['Depth','Superkingdom','Kingdom','Phylum','Class','Order','Family','Genus','Species']]
-        
+        elif k == "Depth":
+            krona_info = df[
+                [
+                    "Depth",
+                    "Superkingdom",
+                    "Kingdom",
+                    "Phylum",
+                    "Class",
+                    "Order",
+                    "Family",
+                    "Genus",
+                    "Species",
+                ]
+            ]
+
         else:
-            print ("""Error: Invalid --krona_mode argument. 
-            Possible arguments: rc, rca, or Depth.""")
-            print ("")
-        
+            print(
+                """Error: Invalid --krona_mode argument. 
+            Possible arguments: rc, rca, or Depth."""
+            )
+            print("")
 
         # remove the unk_xx for better krona representation
-        krona_info = krona_info.replace('unk_.*$', value = '',regex=True)
+        krona_info = krona_info.replace("unk_.*$", value="", regex=True)
 
         # save dataframe to file
         out1 = args.output_fp + ".tsv"
-        pd.DataFrame.to_csv(krona_info, out1, sep='\t', index=False, header=False)
+        pd.DataFrame.to_csv(krona_info, out1, sep="\t", index=False, header=False)
 
         # save krona file
         out2 = args.output_fp + ".html"
@@ -348,48 +508,55 @@
         shell_command = "ktImportText " + out1 + " -o " + out2
         subprocess.run(shell_command, shell=True)
 
-        print ("krona file saved as %s" %(out2))
-        print ("")
-        
+        print("krona file saved as %s" % (out2))
+        print("")
 
     ##### Extended format - calculate read mapping stats
-    if ef == 'y':
-        print ("Calculating read mapping stats...")
-            
+    if ef == "y":
+        print("Calculating read mapping stats...")
+
         with open(mapstat) as mapfile:
-            fragments_line=mapfile.readlines()[3]
-        total_frags = re.split(r'(\t|\n)',fragments_line)[2]
+            fragments_line = mapfile.readlines()[3]
+        total_frags = re.split(r"(\t|\n)", fragments_line)[2]
 
         ### check if the input was single-end or paired end ###
         with open(mapstat) as mapfile:
             kma_command_line = mapfile.readlines()[5]
-        input_file_command = re.split(r'(\t|\n| )',kma_command_line)[6]
+        input_file_command = re.split(r"(\t|\n| )", kma_command_line)[6]
 
         if input_file_command == "-i":
-            print ("""Parsing results based on single-end sequences.""")
+            print("""Parsing results based on single-end sequences.""")
             total_reads = int(total_frags)
         elif input_file_command == "-ipe":
-            print ("""Parsing results based on paired-end sequences.
-            Your data should have the same number of reads in the fwd and rev files.""")
+            print(
+                """Parsing results based on paired-end sequences.
+            Your data should have the same number of reads in the fwd and rev files."""
+            )
 
             total_reads = 2 * int(total_frags)
         else:
-            print ("""Could not identify whether your input files were single- or paired-end. 
-                    Treating as paired-end.""")
+            print(
+                """Could not identify whether your input files were single- or paired-end. 
+                    Treating as paired-end."""
+            )
             total_reads = 2 * int(total_frags)
         ###
 
-        df_stats = pd.read_csv(mapstat, sep='\t', index_col=0, header = 6, encoding='latin1')
-        
+        df_stats = pd.read_csv(
+            mapstat, sep="\t", index_col=0, header=6, encoding="latin1"
+        )
+
         # delete species in df_stats that are not in the CCM result dataframe:
-        df_stats_filt = df_stats[df_stats.index.isin(df.index)].copy() # the copy handles the pandas warning
-        df_stats_filt['perc_map'] = df_stats_filt['readCount'] / total_reads * 100
-        total_mapped = sum(df_stats_filt['perc_map'])
-        
+        df_stats_filt = df_stats[
+            df_stats.index.isin(df.index)
+        ].copy()  # the copy handles the pandas warning
+        df_stats_filt["perc_map"] = df_stats_filt["readCount"] / total_reads * 100
+        total_mapped = sum(df_stats_filt["perc_map"])
+
         stats_out = args.output_fp + "_stats.csv"
-        pd.DataFrame.to_csv(df_stats_filt, stats_out )
-        
-        print ("\nStats file saved as %s" %(stats_out))
-        print ("""\nProportion of reads mapped to the database: %f%%\n""" %(total_mapped))
+        pd.DataFrame.to_csv(df_stats_filt, stats_out)
 
-        
+        print("\nStats file saved as %s" % (stats_out))
+        print(
+            """\nProportion of reads mapped to the database: %f%%\n""" % (total_mapped)
+        )
diff --git a/CCMetagen_extract_seqs.py b/CCMetagen_extract_seqs.py
index 2907564..e801f83 100755
--- a/CCMetagen_extract_seqs.py
+++ b/CCMetagen_extract_seqs.py
@@ -1,15 +1,16 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+
 """
-Script to extract sequence reads (or contigs) assigned to a given taxon
+Script to extract sequence reads (or contigs) assigned to a given taxon.
 
-Created on Wed Jun 10 19:37:36 2020
+Created on Wed Jun 10 19:37:36 2020.
 
 @author: V.R.Marcelino
 """
 
 
-#imports
+# imports
 
 import sys
 import pandas as pd
@@ -19,25 +20,52 @@
 
 # help
 if len(sys.argv) == 1:
-    print ("")
-    print ("""Produces a fasta file containing the sequences assigned to a given taxon""")   
-    print ("For help and options, type: CCMetagen_extract_seqs.py -h")
-    print ("")
+    print("")
+    print(
+        """Produces a fasta file containing the sequences assigned to a given taxon"""
+    )
+    print("For help and options, type: CCMetagen_extract_seqs.py -h")
+    print("")
 
 
 parser = ArgumentParser()
-parser.add_argument('-ifrag', '--input_frags', help='The path to the file containing frags (decompressed)', required=True)
-parser.add_argument('-iccm', '--input_ccmetagen', help='The path to the ccmetagen result csv file', required=True)
-parser.add_argument('-l', '--taxonomic_level', help="""Taxonomic level to merge the results. Options:
-                    Species, Genus, Family, Order, Class, Phylum, Kingdom and Superkingdom""", required=True)
-parser.add_argument('-t', '--taxon', help="""Taxon for which you want to extract sequences. 
-                    Use quotation marks to specify species (e.g. -t 'Escherichia coli')""", required=True)
-parser.add_argument('-o', '--output_fp', default = 'wanted_taxon_seqs', 
-                    help='Path to the output file. Default = wanted_taxon_seqs', required=False)
+parser.add_argument(
+    "-ifrag",
+    "--input_frags",
+    help="The path to the file containing frags (decompressed)",
+    required=True,
+)
+parser.add_argument(
+    "-iccm",
+    "--input_ccmetagen",
+    help="The path to the ccmetagen result csv file",
+    required=True,
+)
+parser.add_argument(
+    "-l",
+    "--taxonomic_level",
+    help="""Taxonomic level to merge the results. Options:
+                    Species, Genus, Family, Order, Class, Phylum, Kingdom and Superkingdom""",
+    required=True,
+)
+parser.add_argument(
+    "-t",
+    "--taxon",
+    help="""Taxon for which you want to extract sequences. 
+                    Use quotation marks to specify species (e.g. -t 'Escherichia coli')""",
+    required=True,
+)
+parser.add_argument(
+    "-o",
+    "--output_fp",
+    default="wanted_taxon_seqs",
+    help="Path to the output file. Default = wanted_taxon_seqs",
+    required=False,
+)
 
 
 args = parser.parse_args()
-frags_fp  = args.input_frags
+frags_fp = args.input_frags
 iccm = args.input_ccmetagen
 tax_level = args.taxonomic_level
 tax_name = args.taxon
@@ -45,14 +73,14 @@
 
 
 # developing and debugging:
-#frags_fp = "test2.frag"
-#iccm = "/test2_CCMetagen_nt_results.csv"
-#tax_name = "Aspergillaceae"
-#tax_level = "Family"
-#o_file = "wanted_seq.fas"
+# frags_fp = "test2.frag"
+# iccm = "/test2_CCMetagen_nt_results.csv"
+# tax_name = "Aspergillaceae"
+# tax_level = "Family"
+# o_file = "wanted_seq.fas"
 
 ## Get a list for all Closest_match where <tax_rank> == <tax_name>:
-df = pd.read_csv(iccm, sep=',', index_col=0)
+df = pd.read_csv(iccm, sep=",", index_col=0)
 
 wanted_df = df[df[tax_level].str.match(tax_name).fillna(False)]
 
@@ -60,24 +88,23 @@
 # Then extract sequences for that taxon.
 ofile_good = open(o_file, "w")
 
-cg = 0 #count good reads
+cg = 0  # count good reads
 with open(frags_fp) as frags:
-    for line in csv.reader(frags, delimiter = '\t'):
+    for line in csv.reader(frags, delimiter="\t"):
         match = line[5]
         seq_id = line[6]
         seq_read = line[0]
 
-        
         if match in wanted_df.index.values:
-            ofile_good.write(">"+ seq_id + "\n" + seq_read + "\n")
-            cg += 1  
+            ofile_good.write(">" + seq_id + "\n" + seq_read + "\n")
+            cg += 1
 
 ofile_good.close()
 
-print ("")
-print ("Done.")
-print ("%i sequences classified as %s were saved in the %s file." %(cg, tax_name, o_file))
-print ("Forward and reverse reads were saved to the same file")
-print ("")
-
-
+print("")
+print("Done.")
+print(
+    "%i sequences classified as %s were saved in the %s file." % (cg, tax_name, o_file)
+)
+print("Forward and reverse reads were saved to the same file")
+print("")
diff --git a/CCMetagen_merge.py b/CCMetagen_merge.py
index cecea71..2998651 100755
--- a/CCMetagen_merge.py
+++ b/CCMetagen_merge.py
@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+
 """
 CCMetagen_merge.py
 
@@ -15,36 +16,67 @@
 CCMetagen_merge.py -i CCMetagen_folder -kr k -l Genus -tlist Candida,Cryptococcus
 
 @ V.R.Marcelino
-Created on Tue Dec 18 10:00:48 2018
+Created on Tue Dec 18 10:00:48 2018.
 """
 
 import sys
-import pandas as pd 
+import pandas as pd
 from argparse import ArgumentParser
 import os
 
 
 parser = ArgumentParser()
-parser.add_argument('-i', '--input_fp', help="""Path to the folder containing CCMetagen text results.
+parser.add_argument(
+    "-i",
+    "--input_fp",
+    help="""Path to the folder containing CCMetagen text results.
                     Note that files must end with ".csv"
-                    and the folder should not contain other .csv files""", required=True)
-parser.add_argument('-t', '--tax_level', default = 'Species',
-                    help="""Taxonomic level to merge the results. Options:
+                    and the folder should not contain other .csv files""",
+    required=True,
+)
+parser.add_argument(
+    "-t",
+    "--tax_level",
+    default="Species",
+    help="""Taxonomic level to merge the results. Options:
                     Closest_match (includes different genes for the same species),
                     Species (Default), Genus, Family, Order, Class, Phylum, Kingdom and Superkingdom
-                    """, required=False)
-parser.add_argument('-o', '--output_fp', default = 'merged_samples', 
-                    help='Path to the output file. Default = merged_samples', required=False)
+                    """,
+    required=False,
+)
+parser.add_argument(
+    "-o",
+    "--output_fp",
+    default="merged_samples",
+    help="Path to the output file. Default = merged_samples",
+    required=False,
+)
 
 # arguments to filter by taxonomy:
-parser.add_argument('-kr', '--keep_or_remove', default = 'n', 
-                    help='keep or remove taxa. Options = k (keep), r (remove) and n (none, default)', required=False)
-
-parser.add_argument('-l', '--filtering_tax_level', default = 'none', 
-                    help='level to perform taxonomic filtering, default = none', required=False)
-
-parser.add_argument('-tlist', '--taxa_list', default = [], type=str, 
-                    help='list taxon names (comma-separated) that you want to keep or exclude', required=False)
+parser.add_argument(
+    "-kr",
+    "--keep_or_remove",
+    default="n",
+    help="keep or remove taxa. Options = k (keep), r (remove) and n (none, default)",
+    required=False,
+)
+
+parser.add_argument(
+    "-l",
+    "--filtering_tax_level",
+    default="none",
+    help="level to perform taxonomic filtering, default = none",
+    required=False,
+)
+
+parser.add_argument(
+    "-tlist",
+    "--taxa_list",
+    default=[],
+    type=str,
+    help="list taxon names (comma-separated) that you want to keep or exclude",
+    required=False,
+)
 
 args = parser.parse_args()
 in_folder = args.input_fp
@@ -54,30 +86,40 @@
 
 
 # debugging:
-#in_folder = "/Users/vmar0011/Documents/Programs_dev/03_CCM_nt"
-#tax_level = "Species"
-#output = "merged_samples_depth_species_no_ecoli"
-#kr = "r" # k for keep, r for remove, and n for none
-#level = "Species"
-#taxa = ["Escherichia coli"]
+# in_folder = "/Users/vmar0011/Documents/Programs_dev/03_CCM_nt"
+# tax_level = "Species"
+# output = "merged_samples_depth_species_no_ecoli"
+# kr = "r" # k for keep, r for remove, and n for none
+# level = "Species"
+# taxa = ["Escherichia coli"]
 
 
 # Set the taxonomic levels to retain info:
-l = ['Superkingdom','Kingdom','Phylum','Class','Order','Family','Genus','Species','Closest_match']
+l = [
+    "Superkingdom",
+    "Kingdom",
+    "Phylum",
+    "Class",
+    "Order",
+    "Family",
+    "Genus",
+    "Species",
+    "Closest_match",
+]
 
 f4agg = {}
 
-f4agg['Depth']= 'sum'
+f4agg["Depth"] = "sum"
 for i in l:
     if i == tax_level:
-        f4agg[i] = 'first'
+        f4agg[i] = "first"
         break
     else:
-        f4agg[i] = 'first'
-    
+        f4agg[i] = "first"
+
 # list of wanted taxa (no depth):
 wanted_taxa = list(f4agg.keys())
-wanted_taxa.remove('Depth')
+wanted_taxa.remove("Depth")
 
 # create new dataframe:
 all_samples = pd.DataFrame()
@@ -86,68 +128,65 @@
 # read input files and merge results
 for file in os.listdir(in_folder):
     if file.endswith(".ccm.csv"):
-        
         sample_name = file.split(".res.ccm.csv")[0]
         result_fp = os.path.join(in_folder, file)
 
-        df = pd.read_csv(result_fp, sep=',', index_col=0)
+        df = pd.read_csv(result_fp, sep=",", index_col=0)
 
-        # this is needed because groupby excludes rows with NAs 
+        # this is needed because groupby excludes rows with NAs
         df = df.fillna("NA")
-        
+
         # keep or remove taxa:
         if kr == "k":
             level = args.filtering_tax_level
             taxa = args.taxa_list.split(",")
-            filt_df_str = "df[df['{}'].isin({})]".format(level,taxa)
+            filt_df_str = "df[df['{}'].isin({})]".format(level, taxa)
             df = eval(filt_df_str)
-            
+
         elif kr == "r":
             level = args.filtering_tax_level
             taxa = args.taxa_list.split(",")
-            filt_df_str = "df[~df['{}'].isin({})]".format(level,taxa)
+            filt_df_str = "df[~df['{}'].isin({})]".format(level, taxa)
             df = eval(filt_df_str)
         elif kr == "n":
-            pass        
+            pass
         else:
-            print ("kr must be k (keep), r (remove) or n (none).")
+            print("kr must be k (keep), r (remove) or n (none).")
             sys.exit("Try again.")
- 
+
         # if tax_level = closest_match, we need an extra column:
-        
-        if tax_level == 'Closest_match':
-            df['Closest_match'] = df.index
-            df.index.name = None # This is needed in modern versions of Pandas (1.0.1)
+
+        if tax_level == "Closest_match":
+            df["Closest_match"] = df.index
+            df.index.name = None  # This is needed in modern versions of Pandas (1.0.1)
 
         depth_by_tax = df.groupby(by=wanted_taxa).agg(f4agg)
-        depth_by_tax.rename(columns={'Depth':sample_name}, inplace=True)
+        depth_by_tax.rename(columns={"Depth": sample_name}, inplace=True)
 
         all_samples = pd.concat([all_samples, depth_by_tax], sort=True, axis=1)
-    
 
-# Group taxon ranks 
-all_samples = all_samples.groupby(by=all_samples.columns, axis = 1).first()
+
+# Group taxon ranks
+all_samples = all_samples.groupby(by=all_samples.columns, axis=1).first()
 
 # add taxon info at the end of the table:
 tax_cols_l = list(f4agg.keys())
-tax_cols_l.remove('Depth')
+tax_cols_l.remove("Depth")
 
 tax_cols = all_samples[tax_cols_l]
 sample_cols = all_samples.drop(columns=tax_cols_l)
 
-all_samples = pd.merge(sample_cols,tax_cols,left_index=True,right_index=True)
+all_samples = pd.merge(sample_cols, tax_cols, left_index=True, right_index=True)
 
 # Fill NaN with zeros:
 all_samples = all_samples.fillna(0)
 
 # Remove the artificially added NAs:
-all_samples = all_samples.replace(["NA"], value = "")
+all_samples = all_samples.replace(["NA"], value="")
 
 
 ### Save
 out = output + ".csv"
 pd.DataFrame.to_csv(all_samples, out, index=False)
 
-print ("Done. Abundance estimates for all samples saved as %s" %(out))
-
-
+print("Done. Abundance estimates for all samples saved as %s" % (out))
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f24a55a..1e2fc38 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,7 +3,12 @@
 This project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html)
 (or tries as best as it can)
 
-## [Unreleased](https://github.com/vrmarcelino/CCMetagen/compare/v1.2.2...master)
+## [v1.4.2](https://github.com/vrmarcelino/CCMetagen/releases/tag/v1.4.2)
+
+- We've recently resumed CCMetagen development. `v1.4.2` is mostly a documental update,
+preparing the software for updates to follow.
+
+# Archive
 
 ## [v1.2.2](https://github.com/vrmarcelino/CCMetagen/compare/v1.1.2...v1.2.2) - 16.06.2020
 
diff --git a/README.md b/README.md
index f0829c9..2bdb3ef 100644
--- a/README.md
+++ b/README.md
@@ -8,81 +8,54 @@ If you this tool, please cite CCMetagen and KMA:
 
   * [Clausen PT, Aarestrup FM, Lund O. 2018. Rapid and precise alignment of raw reads against redundant databases with KMA. BMC bioinformatics. 2018 Dec;19(1):307.](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2336-6)
 
-Besides the guidelines below, we also provide a tutorial to reproduce our metagenome clasisfication analyses of the microbiome of wild birds [here](https://github.com/vrmarcelino/CCMetagen/tree/master/tutorial).
+Besides the guidelines below, we also provide a tutorial to reproduce our metagenome clasisfication analyses of the microbiome of wild birds [here](https://github.com/vrmarcelino/CCMetagen/tree/master/docs/tutorial).
 
 The guidelines below will guide you in using the command-line version of the CCMetagen pipeline.
 
 CCMetagen is also available as a web service at https://cge.food.dtu.dk/services/CCMetagen/.
 Note that we recommend using this command-line version to analyze data exceeding 1.5Gb.
 
-## Requirements and Installation
 
-Make sure you have the dependencies below installed and accessible in your $PATH.
-The guidelines below are for Unix systems.
+## Installation
 
-  * If you do not have it already, download and install [Python 3.6](https://www.python.org/downloads/)
-CCMetagen requires the Python modules [pandas (>0.23)](https://pandas.pydata.org/) and [ETE3](http://etetoolkit.org/).
-The easiest way to install these modules is via conda or pip:
+We recommend installing CCMetagen through [conda](https://docs.conda.io/en/latest/). This will install CCMetagen along with all of its required dependencies.
 
-`conda install pandas`
+After installing conda, you can create a new environment with CCMetagen through the following command:
 
-  * You need a C-compiler and zlib development files to install KMA:
-
-`sudo apt-get install libz-dev`
-
-  * Download and install [KMA](https://bitbucket.org/genomicepidemiology/kma):
 ```
-git clone https://bitbucket.org/genomicepidemiology/kma.git
-cd kma && make
+conda create -n ccmetagen ccmetagen -c bioconda -c conda-forge
 ```
-**Note - a new version of KMA - v1.3.0 – has been released, featuring higher speed and precision. We recommend that you update KMA to v.1.3.0**
 
+You can then activate your environment with:
 
-  * [Krona](https://github.com/marbl/Krona) is required for graphs. To install Krona it in the local folder:
 ```
-wget https://github.com/marbl/Krona/releases/download/v2.7/KronaTools-2.7.tar
-tar xvf KronaTools-2.7.tar 
-cd  KronaTools-2.7
-./install.pl --prefix . 
+conda activate ccmetagen
 ```
 
-  * Then download CCMetagen and add it to your path. You have two options:
+The `-n ccmetagen` flag will name the environment as `ccmetagen`, but you can choose any different name that you'd like. The `ccmetagen` after that specifies that you'd like the CCMetagen package installed into that environment.
+Finally, the `-c bioconda -c conda-forge` specifies that you'd like to use the Bioconda and Conda-Forge channels, which host CCMetagen and its dependencies.
 
-**Install CCMetagen via git:**
-```
-git clone https://github.com/vrmarcelino/CCMetagen
-```
-This will download CCMetagen and the tutorial files.
-You can also just download the python files from this github directory (CCMetagen.py, CCMetagen_merge.py) and the ones in the ccmetagen folder if you rather avoid downloading all other files.
-
-Then add the CCMetagen python scripts to the path, temporarily or permanently. For example:
-`PATH=$PATH<your_folder>/CCMetagen`
-
-To update CCMetagen, go to the CCMetagen folder and type: `git pull`
-
-**Or install CCMetagen via pip:**
-```
-pip install CCMetagen
-```
-This will automatically install the necessary python packages (pandas and ete3), so you can skip that step if you use pip.
+You can also install CCMetagen from source, or using pip, the Python package manager. For that, follow the installation instructions in the deprecated README file in the [docs](https://github.com/vrmarcelino/CCMetagen/tree/master/docs/old_repository_README.md).
 
+Check your CCMetagen installation by running `CCMetagen.py --version` on your command-line.
 
 
 ## Databases
 
-**Option 1** Download the indexed (ready-to-go) nt or RefSeq database [here](http://dx.doi.org/10.25910/5cc7cd40fca8e).
-Download the ncbi_nt_kma.zip file (96GB zipped file, 165GB uncompressed) or the RefSeq_bf.zip (90GB zipped file)
+After installing CCMetagen, you will need a reference database to perform taxonomic classification. There are two ways to obtain this:
+
+**Option 1** Download the indexed (ready-to-go) nt from [here](FIXME).
+Download the ncbi_nt_Dec_2023.zip file (XXGB zipped file, XXXGB uncompressed) or the RefSeq_bf_2024.zip (90GB zipped file)
 Unzip the database, e.g.: `unzip ncbi_nt_kma`.
-The nt database contains the whole in NCBI nucleotide collection (of of Jan 2018), and therefore is suitable to identify a range of microorganisms, including prokaryotes and eukaryotes.
-The RefSeq_bf database contains complete reference bacterial and fungal genomes, suitable for better known habitats such as the human gut or when trying to detect well known species.
+The nt database contains the whole in NCBI nucleotide collection (as of Aug 2023), and therefore is suitable to identify a range of microorganisms, including prokaryotes and eukaryotes.
 
-**Option 2** We have indexed a more recent version of the ncbi nucleotide collection (June 2019) that does not contain environemntal or artificial sequences. The file ncbi_nt_no_env_11jun2019.zip can be found [here](http://dx.doi.org/10.25910/5cc7cd40fca8e) and contains all ncbi nt entries excluding the descendants of environmental eukaryotes (taxid 61964), environmental prokaryotes (48479), unclassified sequences (12908) and artificial sequences (28384).
+There are two versions of the nt database, the one previously mentioned, and another one that does not contain environemntal or artificial sequences. The file ncbi_nt_no_env_Dec2023.zip contains all ncbi nt entries excluding the descendants of environmental eukaryotes (taxid 61964), environmental prokaryotes (48479), unclassified sequences (12908) and artificial sequences (28384).
 
-**Option 3:** Build your own reference database (recommended!)
+**Option 2:** Build your own reference database (recommended!)
 Follow the instructions in the [KMA website](https://bitbucket.org/genomicepidemiology/kma) to index the database.
 It is important that taxids are incorporated in sequence headers for processing with CCMetagen. Sequence headers should look like 
 `>1234|sequence_description`, where 1234 is the taxid. 
-We provide scripts to rename sequences in the nt database [here](https://github.com/vrmarcelino/CCMetagen/tree/master/benchmarking/rename_nt).
+We provide scripts to rename sequences in the nt database [here](https://github.com/vrmarcelino/CCMetagen/tree/master/docs/benchmarking/rename_nt).
 
 If you want to use the RefSeq database, the format is similar to the one required for Kraken. The [Opiniomics blog](http://www.opiniomics.org/building-a-kraken-database-with-new-ftp-structure-and-no-gi-numbers/) describes how to download sequences in an adequate format. Note that you still need to build the index with KMA: `kma_index -i refseq.fna -o refseq_indexed -NI -Sparse -` or `kma_index -i refseq.fna -o refseq_indexed -NI -Sparse TG` for faster analysis.
 
@@ -92,6 +65,7 @@ If you want to use the RefSeq database, the format is similar to the one require
   * First map sequence reads (or contigs) to the database with **KMA**.
 
 For paired-end files:
+
 ```
 kma -ipe $SAMPLE_R1 $SAMPLE_R2 -o sample_out_kma -t_db $db -t $th -1t1 -mem_mode -and -apm f
 ```
@@ -108,32 +82,38 @@ kma -ipe $SAMPLE_R1 $SAMPLE_R2 -o sample_out_kma -t_db $db -t $th -1t1 -mem_mode
 
 Where:
 
-$db is the path to the reference database
-$th is the number of threads
-$SAMPLE_R1 is the path to the mate1 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
-$SAMPLE_R2 is the path to the mate2 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
-$SAMPLE is the path to a single-end metagenome/metatranscriptome file (reads or contigs)
+- `$db` is the path to the reference database
+- `$th` is the number of threads
+- `$SAMPLE_R1` is the path to the mate1 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
+- `$SAMPLE_R2` is the path to the mate2 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
+- `$SAMPLE` is the path to a single-end metagenome/metatranscriptome file (reads or contigs)
 
+Then run `CCMetagen.py`:
 
-  * Then run **CCMetagen**:
 ```
 CCMetagen.py -i $sample_out_kma.res -o results
 ```
+
 Where $sample_out_kma.res is alignment results produced by KMA.
 
 Note that if you are running CCMetagen from the local folder (instead of adding it to your path), you may need to add 'python' before CCMetagen: `python CCMetagen.py -i $sample_out_kma.res -o results`
 
 Done! This will make an additional quality filter and output a text file with ranked taxonomic classifications and a krona graph file for interactive visualization.
 
-An example of the CCMetagen output can be found [here (.csv file)](https://github.com/vrmarcelino/CCMetagen/blob/master/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv) and [here (.html file)](https://htmlpreview.github.io/?https://github.com/vrmarcelino/CCMetagen/blob/master/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html).
+An example of the CCMetagen output can be found [here (.csv file)](https://github.com/vrmarcelino/CCMetagen/blob/master/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv) and [here (.html file)](https://htmlpreview.github.io/?https://github.com/vrmarcelino/CCMetagen/blob/master/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html).
 
-<img src=tutorial/figs_tutorial/krona_photo.png width="500" height="419.64">
+<img src=docs/tutorial/figs_tutorial/krona_photo.png width="500" height="419.64">
 
 In the .csv file, you will find the depth (abundance) of each match.
 
 ## Abundance units
 
-**Depth can be estimated in four ways:** by counting the number of nucleotides matching the reference sequence (use flag --depth_unit nc), by applying an additional correction for template length (default in KMA and CCMetagen), by calculating depth in Reads Per Million (RPM, use flag --depth_unit rpm), or by counting the number of fragments (i.e. number of PE reads matching to teh reference sequence, use flag --depth_unit fr). If you want RPM or fragment units, you will need to suply the .mapstats file generated with KMA (which you get when running kma with the flag '-ef').
+**Depth can be estimated in four ways:** 
+
+1. By applying an additional correction for template length (default in KMA and CCMetagen);
+2. By counting the number of nucleotides matching the reference sequence (use flag --depth_unit nc);
+3. By calculating depth in Reads Per Million (RPM, use flag --depth_unit rpm); or
+4. By counting the number of fragments (i.e. number of PE reads matching to teh reference sequence, use flag --depth_unit fr). If you want RPM or fragment units, you will need to suply the .mapstats file generated with KMA (which you get when running kma with the flag '-ef').
 
 
 ## Balancing sensitivity and specificity
@@ -171,7 +151,7 @@ CCMetagen.py -h
 CCMetagen_merge.py -i $CCMetagen_out
 ```
 
-Where $CCMetagen_out is the folder containing the CCMetagen taxonomic classifications.
+Where `$CCMetagen_out` is the folder containing the CCMetagen taxonomic classifications.
 The results must be in .csv format (default or '--mode text' output of CCMetagen), and these files **must end in ".ccm.csv"**.
 
 The flag '-t' define the taxonomic level to merge the results. The default is species-level.
@@ -214,7 +194,7 @@ Ex: Generate a fasta file containing all sequences that mapped to the genus Esch
 CCMetagen_extract_seqs.py -iccm $CCMetagen_out -ifrag $sample_out_kma.frag -l Genus -t Eschericha
 ```
 
-Where $CCMetagen_out is the .csv file generated with CCMetagen and $sample_out_kma.frag is the .frag file generated with KMA. The frag file needs to be decompressed: `gunzip *.frag.gz`
+Where `$CCMetagen_out` is the .csv file generated with CCMetagen and `$sample_out_kma.frag` is the .frag file generated with KMA. The frag file needs to be decompressed: `gunzip *.frag.gz`
 
 For species-level filtering (where there is a space in taxon names), use quotation marks.
 Ex: Generate a fasta file containing all sequences that mapped to _E. coli_:
@@ -223,7 +203,7 @@ CCMetagen_extract_seqs.py -iccm $CCMetagen_out -ifrag $sample_out_kma.frag -l Sp
 ```
 
 
-**Check out our [tutorial](https://github.com/vrmarcelino/CCMetagen/tree/master/tutorial) for an applied example of the CCMetagen pipeline.**
+**Check out our [tutorial](https://github.com/vrmarcelino/CCMetagen/tree/master/docs/tutorial) for an applied example of the CCMetagen pipeline.**
 
 
 
diff --git a/benchmarking/rename_nt/rename_nt.py b/benchmarking/rename_nt/rename_nt.py
deleted file mode 100644
index 46a7880..0000000
--- a/benchmarking/rename_nt/rename_nt.py
+++ /dev/null
@@ -1,50 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Script to add taxids to nt collection
-
-@ V.R.Marcelino
-Created on Fri Dec 28 10:37:56 2018
-"""
-
-import re
-
-
-# input files:
-acc2taxid_map="accession_taxid_nucl.map"
-nt_in="nt.fa"
-
-
-# function to  get taxids from accession numbers
-def get_tax_id_dic (accession, accession_dic):
-    taxid = accession_dic.get(accession)
-    if taxid == None:
-        taxid = "unk_taxid"
-    return taxid
-
-
-# read acc2taxid.map as dictionary
-with open(acc2taxid_map) as a:
-    acc2tax_dic = dict(x.rstrip().split(None, 1) for x in a)
-    
-
-# read fasta file and output new file on the fly
-# although it says 'open', it reads line by line and won't store the whole file in memory.
-    
-new_nt = open('nt_w_taxid.fas', 'w')
-
-with open(nt_in) as nt:
-    for line in nt:
-        if line.startswith(">"):
-            splited_rec = re.split (r'(>| )', line)
-            accession = splited_rec[2]
-            taxid = get_tax_id_dic(accession,acc2tax_dic)
-            line = ">" + taxid + "|" + "".join(splited_rec[2:])
-        new_nt.write(line)
-            
-new_nt.close()
-
-
-print ("Done!")
-
-
diff --git a/ccmetagen/__init__.py b/ccmetagen/__init__.py
index 7f8c61b..ff8cc66 100755
--- a/ccmetagen/__init__.py
+++ b/ccmetagen/__init__.py
@@ -1,1 +1,29 @@
-__version__ = "1.4.3"
\ No newline at end of file
+__version__ = "1.5.0"
+
+# This can be configured
+_RANK_THRESHOLDS = {
+    "species_threshold": 98.41, # Yeast - Vu et al 2016
+    "genus_threshold": 96.31,   # Yeast - Vu et al 2016
+    "family_threshold": 88.51,  # Filamentous fungi - Vu et al 2019
+    "order_threshold": 81.21,   # Filamentous fungi - Vu et al 2019
+    "class_threshold": 80.91,   # Filamentous fungi - Vu et al 2019
+    "phylum_threshold": 0,      # no data, no filtering    
+}
+_ARGPARSE_DEFAULTS = {
+    "mode": "both",
+    "output_fp": "CCMetagen_out",
+    "reference_database": "nt",
+    "extended_output_file": "n",
+    "depth_unit": "kma",
+    "mapstat": None,
+    "depth": 0.2,
+    "coverage": 20,
+    "query_identity": 50,
+    "pvalue": 0.05,
+    "krona_mode": "Depth",
+    "local_taxfile": None,
+    "turn_off_sim_thresholds": "n",
+    **_RANK_THRESHOLDS,
+}
+
+__all__ = ["__version__", "_ARGPARSE_DEFAULTS", "_RANK_THRESHOLDS"]
diff --git a/ccmetagen/cTaxInfo.py b/ccmetagen/cTaxInfo.py
index c2603d6..c923f0f 100755
--- a/ccmetagen/cTaxInfo.py
+++ b/ccmetagen/cTaxInfo.py
@@ -1,28 +1,45 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+
 """
 Classes that store taxonomic info from matches and NCBI lineage
 
 @ V.R.Marcelino
-Created on 12 Jul 2018
-Last update on 27 March 2019
-
+Created on 12 Jul 2018.
 """
 
+
 ### Stores taxid and other info
-class TaxInfo(): 
-    
-    def __init__(self, TaxId=None, Lineage=None, Sample=None, RefDatabase=None, 
-                 Abundance=None, Superkingdom = None, Superkingdom_TaxId = None, 
-                 Kingdom=None, Kingdom_TaxId=None, Phylum=None, 
-                 Phylum_TaxId=None, Class=None, Class_TaxId=None,
-                 Order=None, Order_TaxId=None, Family=None, Family_TaxId=None,
-                 Genus=None, Genus_TaxId=None, Species=None, Species_TaxId=None,
-                 LCA_TaxId=None, Coverage=None):
-        
+class TaxInfo:
+    def __init__(
+        self,
+        TaxId=None,
+        Lineage=None,
+        Sample=None,
+        RefDatabase=None,
+        Abundance=None,
+        Superkingdom=None,
+        Superkingdom_TaxId=None,
+        Kingdom=None,
+        Kingdom_TaxId=None,
+        Phylum=None,
+        Phylum_TaxId=None,
+        Class=None,
+        Class_TaxId=None,
+        Order=None,
+        Order_TaxId=None,
+        Family=None,
+        Family_TaxId=None,
+        Genus=None,
+        Genus_TaxId=None,
+        Species=None,
+        Species_TaxId=None,
+        LCA_TaxId=None,
+        Coverage=None,
+    ):
         # info from matches
         self.TaxId = TaxId
-        self.Lineage = Lineage # the results of the match, not considering taxIDs
+        self.Lineage = Lineage  # the results of the match, not considering taxIDs
         self.Sample = Sample
         self.RefDatabase = RefDatabase
         self.Abundance = Abundance
@@ -38,13 +55,13 @@ def __init__(self, TaxId=None, Lineage=None, Sample=None, RefDatabase=None,
         self.Class_TaxId = Class_TaxId
         self.Order = Order
         self.Order_TaxId = Order_TaxId
-        self.Family = Order
-        self.Family_TaxId = Order_TaxId
+        self.Family = Family
+        self.Family_TaxId = Family_TaxId
         self.Genus = Genus
         self.Genus_TaxId = Genus_TaxId
         self.Species = Species
         self.Species_TaxId = Species_TaxId
-        
+
         # Taxid of the Lowest Common Ancestor
         self.LCA_TaxId = LCA_TaxId
 
diff --git a/ccmetagen/fNCBItax.py b/ccmetagen/fNCBItax.py
index 7d90a54..3fed8ed 100755
--- a/ccmetagen/fNCBItax.py
+++ b/ccmetagen/fNCBItax.py
@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+
 """
 Function that takes a taxid as input and returns a class containing the taxonomic info and taxids
 for multiple taxonomic ranks.
@@ -7,16 +8,23 @@
 Required to parse the results of KMA and store them in the SQLite3 'bench.db'
 
 @ V.R.Marcelino
-Created on 12 Jul 2018
-
+Created on 12 Jul 2018.
 """
 
 from ete3 import NCBITaxa
-from ccmetagen import cTaxInfo  # where we define classes used here
 
 
 def lineage_extractor(query_taxid, TaxInfo_object, taxfile=None):
-    list_of_taxa_ranks = ['superkingdom', 'kingdom', 'phylum', 'class', 'order', 'family','genus', 'species']
+    list_of_taxa_ranks = [
+        "Superkingdom",
+        "Kingdom",
+        "Phylum",
+        "Class",
+        "Order",
+        "Family",
+        "Genus",
+        "Species",
+    ]
     if taxfile is not None:
         ncbi = NCBITaxa(taxfile)
     else:
@@ -25,66 +33,17 @@ def lineage_extractor(query_taxid, TaxInfo_object, taxfile=None):
     ranks = ncbi.get_rank(lineage)
     names = ncbi.get_taxid_translator(lineage)
 
-# get known data
-
+    # get known data
     for key, val in ranks.items():
-
-        if val == list_of_taxa_ranks[0]:
-            TaxInfo_object.Superkingdom = names[key]
-            TaxInfo_object.Superkingdom_TaxId = key
-
-        elif val == list_of_taxa_ranks[1]:
-            TaxInfo_object.Kingdom = names[key]
-            TaxInfo_object.Kingdom_TaxId = key
-
-        elif val == list_of_taxa_ranks[2]:
-            TaxInfo_object.Phylum = names[key]
-            TaxInfo_object.Phylum_TaxId = key
-
-        elif val == list_of_taxa_ranks[3]:
-            TaxInfo_object.Class = names[key]
-            TaxInfo_object.Class_TaxId = key
-
-        elif val == list_of_taxa_ranks[4]:
-            TaxInfo_object.Order = names[key]
-            TaxInfo_object.Order_TaxId = key
-
-        elif val == list_of_taxa_ranks[5]:
-            TaxInfo_object.Family = names[key]
-            TaxInfo_object.Family_TaxId = key
-
-        elif val == list_of_taxa_ranks[6]:
-            TaxInfo_object.Genus = names[key]
-            TaxInfo_object.Genus_TaxId = key
-
-        elif val == list_of_taxa_ranks[7]:
-            TaxInfo_object.Species = names[key]
-            TaxInfo_object.Species_TaxId = key
-
-# fill in the blanks
-    if TaxInfo_object.Superkingdom is None:
-        TaxInfo_object.Superkingdom = "unk_sk"
-
-    if TaxInfo_object.Kingdom is None:
-        TaxInfo_object.Kingdom = "unk_k"
-
-    if TaxInfo_object.Phylum is None:
-        TaxInfo_object.Phylum = "unk_p"
-
-    if TaxInfo_object.Class is None:
-        TaxInfo_object.Class = "unk_c"
-
-    if TaxInfo_object.Order is None:
-        TaxInfo_object.Order = "unk_o"
-
-    if TaxInfo_object.Family is None:
-        TaxInfo_object.Family = "unk_f"
-
-    if TaxInfo_object.Genus is None:
-        TaxInfo_object.Genus = "unk_g"
-
-    if TaxInfo_object.Species is None:
-        TaxInfo_object.Species = "unk_s"
-
+        for rank in list_of_taxa_ranks:
+            if val == rank:
+                setattr(TaxInfo_object, rank, names[key])
+                setattr(TaxInfo_object, rank + "_TaxId", key)
+
+    # fill in the blanks
+    for attr in list_of_taxa_ranks:
+        if getattr(TaxInfo_object, attr) is None:
+            initial = attr.lower()[0] if attr != "Superkingdom" else "sk"
+            setattr(TaxInfo_object, attr, "unk_" + initial)
 
     return TaxInfo_object
diff --git a/ccmetagen/fParseKMA.py b/ccmetagen/fParseKMA.py
index 0add3df..58a8a43 100755
--- a/ccmetagen/fParseKMA.py
+++ b/ccmetagen/fParseKMA.py
@@ -1,53 +1,61 @@
 #!/usr/bin/env ipython
 # -*- coding: utf-8 -*-
+
 """
 Functions to parse KMA results
 
 @ V.R.Marcelino
-Created on 1 Aug 2018
-
-
+Created on 1 Aug 2018.
 """
+
 import re
 
 # local imports
-from ccmetagen import cTaxInfo
-from ccmetagen import fNCBItax
-
-import subprocess
+from ccmetagen import cTaxInfo, fNCBItax, _RANK_THRESHOLDS
 
 # function to filter a res file in pandas df format:
-def res_filter(df,ref_database, cov,Iden,Depth,p):
+def res_filter(df, ref_database, cov, Iden, Depth, p):
     df = df.drop(df[df.Template_Coverage < cov].index)
 
-    # filter based on identity      
+    # filter based on identity
     df = df.drop(df[df.Query_Identity < Iden].index)
-    
+
     # filter based on depth
     df = df.drop(df[df.Depth < Depth].index)
 
     # filter based on p-values
     df = df.drop(df[df.p_value > p].index)
-     
-    return df
 
+    return df
 
-# function that takes as input a pandas dataframe with KMA results 
-# and add tax information to results 
-def populate_w_tax(in_df, ref_database,species_threshold,genus_threshold,
-                   family_threshold,order_threshold,class_threshold,phylum_threshold, taxfile=None):
-    #defaults:
-    #species_threshold = 98.41 # Yeast - Vu et al 2016
-    #genus_threshold = 96.31 # Yeast - Vu et al 2016
-    #family_threshold = 88.51 # Filamentous fungi - Vu et al 2019
-    #order_threshold = 81.21 # Filamentous fungi - Vu et al 2019
-    #class_threshold = 80.91 # Filamentous fungi - Vu et al 2019
-    #phylum, kingdom and superkingdom = 0  # no data, no filtering
 
+# function that takes as input a pandas dataframe with KMA results
+# and add tax information to results
+def populate_w_tax(
+    in_df,
+    ref_database,
+    species_threshold,
+    genus_threshold,
+    family_threshold,
+    order_threshold,
+    class_threshold,
+    phylum_threshold,
+    taxfile=None,
+):
+    # For default thresholds, see ccmetagen/__init__.py
 
     # Make sure all taxa columns are strings (doesn't automatically happen if the first one is None)
-    in_df = in_df.assign(LCA_TaxId="",Superkingdom="",Kingdom="",Phylum="",Class="",Order="",Family="",Genus="",Species="")
-
+    in_df = in_df.assign(
+        LCA_TaxId="",
+        Superkingdom="",
+        Kingdom="",
+        Phylum="",
+        Class="",
+        Order="",
+        Family="",
+        Genus="",
+        Species="",
+    )
 
     # index == the #template (fungal match)
     for index, row in in_df.iterrows():
@@ -55,102 +63,94 @@ def populate_w_tax(in_df, ref_database,species_threshold,genus_threshold,
 
         # define the tax. rank based on similarity:
         if ref_database == "UNITE":
-            split_match = re.split (r'(\|| )', index)
-            qiden = row['Query_Identity']
+            split_match = re.split(r"(\|| )", index)
+            qiden = row["Query_Identity"]
             match_info.Lineage = split_match[12]
 
             # if taxid is knwon:
-            if split_match[4] != 'unk_taxid':
-                
+            if split_match[4] != "unk_taxid":
                 match_info.TaxId = int(split_match[4])
-                match_info = fNCBItax.lineage_extractor(match_info.TaxId , match_info, taxfile)
-                
-                # Warning about unknown taxids: 
-            else:
-                print ("")
-                print ("WARNING: based on accession number, no taxonomic information was found in NCBI for %s" %(match_info.Lineage))
-                print ("This match will not get NCBItax taxonomic ranks")
-                print ("")
-                match_info.TaxId = split_match[4] # 'unk_taxid'
+                match_info = fNCBItax.lineage_extractor(
+                    match_info.TaxId, match_info, taxfile
+                )
 
+                # Warning about unknown taxids:
+            else:
+                print("")
+                print(
+                    "WARNING: based on accession number, no taxonomic information was found in NCBI for %s"
+                    % (match_info.Lineage)
+                )
+                print("This match will not get NCBItax taxonomic ranks")
+                print("")
+                match_info.TaxId = split_match[4]  # 'unk_taxid'
 
         elif ref_database == "RefSeq":
-            split_match = re.split (r'(\|| )', index)
-            qiden = row['Query_Identity']
+            split_match = re.split(r"(\|| )", index)
+            qiden = row["Query_Identity"]
             match_info.TaxId = int(split_match[4])
             species = split_match[6] + " " + split_match[8]
             match_info.Lineage = species
             # include info from NCBI:
-            match_info = fNCBItax.lineage_extractor(match_info.TaxId, match_info, taxfile)
-
+            match_info = fNCBItax.lineage_extractor(
+                match_info.TaxId, match_info, taxfile
+            )
 
         elif ref_database == "nt":
-            split_match = re.split (r'(\|| )', index)
-            qiden = row['Query_Identity']
+            split_match = re.split(r"(\|| )", index)
+            qiden = row["Query_Identity"]
             match_info.Lineage = split_match[2]
-            
-            #get taxid from accession number
+
+            # get taxid from accession number
             taxid = split_match[0]
 
-            if taxid == 'unk_taxid':
-                # Warning about unknown taxids: 
-                print ("")
-                print ("WARNING: no NCBI's taxid found for accession %s" %(match_info.Lineage))
-                print ("This match will not get taxonomic ranks")
-                print ("")
-                
+            if taxid == "unk_taxid":
+                # Warning about unknown taxids:
+                print("")
+                print(
+                    "WARNING: no NCBI's taxid found for accession %s"
+                    % (match_info.Lineage)
+                )
+                print("This match will not get taxonomic ranks")
+                print("")
+
             else:
-                match_info.TaxId = int(taxid)            
-                match_info = fNCBItax.lineage_extractor(match_info.TaxId, match_info, taxfile)
-   
-            
-        # Populate the df with lineage info and the LCA taxid: 
-        
-        in_df.at[index, 'Superkingdom'] = match_info.Superkingdom
-        in_df.at[index, 'Kingdom'] = match_info.Kingdom
-        
+                match_info.TaxId = int(taxid)
+                match_info = fNCBItax.lineage_extractor(
+                    match_info.TaxId, match_info, taxfile
+                )
+
+        # Populate the df with lineage info and the LCA taxid:
+
+        in_df.at[index, "Superkingdom"] = match_info.Superkingdom
+        in_df.at[index, "Kingdom"] = match_info.Kingdom
+
         # Assign LCA_taxid. Go to Kingdom if possible:
-        in_df.at[index, 'LCA_TaxId'] = match_info.Superkingdom_TaxId
-        
+        in_df.at[index, "LCA_TaxId"] = match_info.Superkingdom_TaxId
+
         if match_info.Kingdom_TaxId is not None:
-            in_df.at[index, 'LCA_TaxId'] = match_info.Kingdom_TaxId
-        
+            in_df.at[index, "LCA_TaxId"] = match_info.Kingdom_TaxId
+
         # if it matches to uncultured or unclassified fungus, use the Fungi LCA itaxid:
-        if match_info.Kingdom == 'Fungi':
-            in_df.at[index, 'LCA_TaxId'] = 4751
-
-        # fill in the rest of the table according to similarity threshold:
-        if qiden >= phylum_threshold:
-            in_df.at[index, 'Phylum'] = match_info.Phylum
-            if match_info.Phylum_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Phylum_TaxId
-
-        if qiden >= class_threshold:
-            in_df.at[index, 'Class'] = match_info.Class
-            if match_info.Class_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Class_TaxId
-    
-        if qiden >= order_threshold:
-            in_df.at[index, 'Order'] = match_info.Order
-            if match_info.Order_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Order_TaxId
-
-        if qiden >= family_threshold:
-            in_df.at[index, 'Family'] = match_info.Family
-            if match_info.Family_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Family_TaxId
-    
-        if qiden >= genus_threshold:
-            in_df.at[index, 'Genus'] = match_info.Genus
-            if match_info.Genus_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Genus_TaxId
-    
-        if qiden >= species_threshold:
-            in_df.at[index, 'Species'] = match_info.Species
-            if match_info.Species_TaxId != None:
-                in_df.at[index, 'LCA_TaxId'] = match_info.Species_TaxId
-    
-    return in_df
+        if match_info.Kingdom == "Fungi":
+            in_df.at[index, "LCA_TaxId"] = 4751
+
+        thresholds = {
+            "Phylum": phylum_threshold,
+            "Class": class_threshold,
+            "Order": order_threshold,
+            "Family": family_threshold,
+            "Genus": genus_threshold,
+            "Species": species_threshold
+        }
+
+        for rank, threshold in thresholds.items():
+            if qiden >= threshold:
+                tax_info_attr = f"{rank}_TaxId"
 
+                in_df.at[index, rank] = getattr(match_info, rank)
 
+                if getattr(match_info, tax_info_attr) is not None:
+                    in_df.at[index, "LCA_TaxId"] = getattr(match_info, tax_info_attr)
 
+        return in_df
diff --git a/benchmarking/README.md b/docs/benchmarking/README.md
similarity index 90%
rename from benchmarking/README.md
rename to docs/benchmarking/README.md
index be353a2..3f7e596 100644
--- a/benchmarking/README.md
+++ b/docs/benchmarking/README.md
@@ -1,5 +1,8 @@
 # Benchmarking
 
+**Notice:** Code in this package was used for the benchmarking of CCMetagen prior to publication. It has now been
+deprecated.
+
 ## Containers
 
 ### Repository layout
diff --git a/benchmarking/convert2sql/CCMetagen2SQL.py b/docs/benchmarking/convert2sql/CCMetagen2SQL.py
similarity index 100%
rename from benchmarking/convert2sql/CCMetagen2SQL.py
rename to docs/benchmarking/convert2sql/CCMetagen2SQL.py
diff --git a/benchmarking/convert2sql/Centrifuge2SQL.py b/docs/benchmarking/convert2sql/Centrifuge2SQL.py
similarity index 100%
rename from benchmarking/convert2sql/Centrifuge2SQL.py
rename to docs/benchmarking/convert2sql/Centrifuge2SQL.py
diff --git a/benchmarking/convert2sql/KMA2SQL.py b/docs/benchmarking/convert2sql/KMA2SQL.py
similarity index 100%
rename from benchmarking/convert2sql/KMA2SQL.py
rename to docs/benchmarking/convert2sql/KMA2SQL.py
diff --git a/benchmarking/convert2sql/Kraken22SQL.py b/docs/benchmarking/convert2sql/Kraken22SQL.py
similarity index 100%
rename from benchmarking/convert2sql/Kraken22SQL.py
rename to docs/benchmarking/convert2sql/Kraken22SQL.py
diff --git a/benchmarking/convert2sql/KrakenUniq2SQL.py b/docs/benchmarking/convert2sql/KrakenUniq2SQL.py
similarity index 100%
rename from benchmarking/convert2sql/KrakenUniq2SQL.py
rename to docs/benchmarking/convert2sql/KrakenUniq2SQL.py
diff --git a/benchmarking/convert2sql/PrecisionRecall.py b/docs/benchmarking/convert2sql/PrecisionRecall.py
similarity index 100%
rename from benchmarking/convert2sql/PrecisionRecall.py
rename to docs/benchmarking/convert2sql/PrecisionRecall.py
diff --git a/benchmarking/convert2sql/PrecisionRecall_bac.py b/docs/benchmarking/convert2sql/PrecisionRecall_bac.py
similarity index 100%
rename from benchmarking/convert2sql/PrecisionRecall_bac.py
rename to docs/benchmarking/convert2sql/PrecisionRecall_bac.py
diff --git a/benchmarking/convert2sql/PrecisionRecall_per_RefDB.py b/docs/benchmarking/convert2sql/PrecisionRecall_per_RefDB.py
similarity index 100%
rename from benchmarking/convert2sql/PrecisionRecall_per_RefDB.py
rename to docs/benchmarking/convert2sql/PrecisionRecall_per_RefDB.py
diff --git a/benchmarking/convert2sql/cResInfo.py b/docs/benchmarking/convert2sql/cResInfo.py
similarity index 100%
rename from benchmarking/convert2sql/cResInfo.py
rename to docs/benchmarking/convert2sql/cResInfo.py
diff --git a/benchmarking/convert2sql/cTaxInfo.py b/docs/benchmarking/convert2sql/cTaxInfo.py
similarity index 100%
rename from benchmarking/convert2sql/cTaxInfo.py
rename to docs/benchmarking/convert2sql/cTaxInfo.py
diff --git a/benchmarking/convert2sql/fNCBItax.py b/docs/benchmarking/convert2sql/fNCBItax.py
similarity index 95%
rename from benchmarking/convert2sql/fNCBItax.py
rename to docs/benchmarking/convert2sql/fNCBItax.py
index 21a6b55..78b1bc1 100755
--- a/benchmarking/convert2sql/fNCBItax.py
+++ b/docs/benchmarking/convert2sql/fNCBItax.py
@@ -1,6 +1,8 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
+UPDATE v1.4.3 - This script has been deprecated. See ccmetagen/fNCBItax.py for the latest version.
+
 Function that takes a taxid as input and returns a class containing the taxonomic info and taxids 
 for multiple taxonomic ranks. 
 
diff --git a/benchmarking/convert2sql/fParseKMA.py b/docs/benchmarking/convert2sql/fParseKMA.py
similarity index 100%
rename from benchmarking/convert2sql/fParseKMA.py
rename to docs/benchmarking/convert2sql/fParseKMA.py
diff --git a/benchmarking/convert2sql/fQueries.py b/docs/benchmarking/convert2sql/fQueries.py
similarity index 100%
rename from benchmarking/convert2sql/fQueries.py
rename to docs/benchmarking/convert2sql/fQueries.py
diff --git a/benchmarking/convert2sql/produce_fungi_table.py b/docs/benchmarking/convert2sql/produce_fungi_table.py
similarity index 100%
rename from benchmarking/convert2sql/produce_fungi_table.py
rename to docs/benchmarking/convert2sql/produce_fungi_table.py
diff --git a/benchmarking/rename_nt/README.md b/docs/benchmarking/rename_nt/README.md
similarity index 100%
rename from benchmarking/rename_nt/README.md
rename to docs/benchmarking/rename_nt/README.md
diff --git a/docs/benchmarking/rename_nt/rename_nt.py b/docs/benchmarking/rename_nt/rename_nt.py
new file mode 100644
index 0000000..262853b
--- /dev/null
+++ b/docs/benchmarking/rename_nt/rename_nt.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script to add taxids to nt collection
+
+@ V.R.Marcelino
+Created on Fri Dec 28 10:37:56 2018
+"""
+import re
+import logging
+
+
+logging.info("Starting script.")
+
+# Input files:
+acc2taxid_map = "accession_taxid_nucl.map"
+nt_in = "nt_sequential.fa"
+
+# Output file:
+nt_out = "nt_w_taxid.fas"
+
+# Function to get taxids from accession numbers
+def get_tax_id_dic(accession, accession_dic):
+    taxid = accession_dic.get(accession, "unk_taxid")
+    return taxid
+
+# Read acc2taxid.map as a dictionary
+logging.info("Reading acc2taxid.map as a dictionary.")
+with open(acc2taxid_map) as a:
+    acc2tax_dic = dict(x.rstrip().split(None, 1) for x in a)
+
+logging.info(f"Dictionary has {len(acc2tax_dic)} entries.")
+
+# Read fasta file and write the new file on the fly
+logging.info("Reading nt file and writing new file.")
+with open(nt_in) as nt, open(nt_out, 'w') as new_nt:
+    ix = 0
+    for line in nt:
+        ix += 1
+        if not ix % 1000000:
+            logging.info("Processed {} lines.".format(ix))
+
+        if line.startswith(">"):
+            splited_rec = re.split(r'(>| )', line)
+            accession = splited_rec[2]
+            taxid = get_tax_id_dic(accession, acc2tax_dic)
+            line = ">" + taxid + "|" + "".join(splited_rec[2:])
+        new_nt.write(line)
+
+logging.info("All done.")
diff --git a/benchmarking/time_GNU/calc_time.py b/docs/benchmarking/time_GNU/calc_time.py
similarity index 100%
rename from benchmarking/time_GNU/calc_time.py
rename to docs/benchmarking/time_GNU/calc_time.py
diff --git a/docs/old_repository_README.md b/docs/old_repository_README.md
new file mode 100644
index 0000000..c60fc5a
--- /dev/null
+++ b/docs/old_repository_README.md
@@ -0,0 +1,345 @@
+# CCMetagen
+
+**Disclaimer:** this is CCMetagen's old repository README file. It has been deprecated as of `v1.4.3`. Although it can be accessed through the Git history, it has been preserved here to facilitate look-up for users.
+
+CCMetagen processes sequence alignments produced with [KMA](https://bitbucket.org/genomicepidemiology/kma), which implements the ConClave sorting scheme to achieve highly accurate read mappings. The pipeline is fast enough to use the whole NCBI nt collection as reference, facilitating the inclusion of understudied organisms, such as microbial eukaryotes, in metagenome surveys. CCMetagen produces ranked taxonomic results in user-friendly formats that are ready for publication or downstream statistical analyses.
+
+If you this tool, please cite CCMetagen and KMA:
+
+  * [Marcelino VR, Clausen PT, Buchman J, Wille M, Iredell JR, Meyer W, Lund O, Sorrell T, Holmes EC. 2019. CCMetagen: comprehensive and accurate identification of eukaryotes and prokaryotes in metagenomic data. Genome Biology. 2020 Dec;21(1):1-5.](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02014-2)
+
+  * [Clausen PT, Aarestrup FM, Lund O. 2018. Rapid and precise alignment of raw reads against redundant databases with KMA. BMC bioinformatics. 2018 Dec;19(1):307.](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2336-6)
+
+Besides the guidelines below, we also provide a tutorial to reproduce our metagenome clasisfication analyses of the microbiome of wild birds [here](https://github.com/vrmarcelino/CCMetagen/tree/master/tutorial).
+
+The guidelines below will guide you in using the command-line version of the CCMetagen pipeline.
+
+CCMetagen is also available as a web service at https://cge.food.dtu.dk/services/CCMetagen/.
+Note that we recommend using this command-line version to analyze data exceeding 1.5Gb.
+
+## Requirements and Installation
+
+Make sure you have the dependencies below installed and accessible in your $PATH.
+The guidelines below are for Unix systems.
+
+  * If you do not have it already, download and install [Python 3.6](https://www.python.org/downloads/)
+CCMetagen requires the Python modules [pandas (>0.23)](https://pandas.pydata.org/) and [ETE3](http://etetoolkit.org/).
+The easiest way to install these modules is via conda or pip:
+
+`conda install pandas`
+
+  * You need a C-compiler and zlib development files to install KMA:
+
+`sudo apt-get install libz-dev`
+
+  * Download and install [KMA](https://bitbucket.org/genomicepidemiology/kma):
+```
+git clone https://bitbucket.org/genomicepidemiology/kma.git
+cd kma && make
+```
+**Note - a new version of KMA - v1.3.0 – has been released, featuring higher speed and precision. We recommend that you update KMA to v.1.3.0**
+
+
+  * [Krona](https://github.com/marbl/Krona) is required for graphs. To install Krona it in the local folder:
+```
+wget https://github.com/marbl/Krona/releases/download/v2.7/KronaTools-2.7.tar
+tar xvf KronaTools-2.7.tar 
+cd  KronaTools-2.7
+./install.pl --prefix . 
+```
+
+  * Then download CCMetagen and add it to your path. You have two options:
+
+**Install CCMetagen via git:**
+```
+git clone https://github.com/vrmarcelino/CCMetagen
+```
+This will download CCMetagen and the tutorial files.
+You can also just download the python files from this github directory (CCMetagen.py, CCMetagen_merge.py) and the ones in the ccmetagen folder if you rather avoid downloading all other files.
+
+Then add the CCMetagen python scripts to the path, temporarily or permanently. For example:
+`PATH=$PATH<your_folder>/CCMetagen`
+
+To update CCMetagen, go to the CCMetagen folder and type: `git pull`
+
+**Or install CCMetagen via pip:**
+```
+pip install CCMetagen
+```
+This will automatically install the necessary python packages (pandas and ete3), so you can skip that step if you use pip.
+
+
+
+## Databases
+
+**Option 1** Download the indexed (ready-to-go) nt or RefSeq database [here](http://dx.doi.org/10.25910/5cc7cd40fca8e).
+Download the ncbi_nt_kma.zip file (96GB zipped file, 165GB uncompressed) or the RefSeq_bf.zip (90GB zipped file)
+Unzip the database, e.g.: `unzip ncbi_nt_kma`.
+The nt database contains the whole in NCBI nucleotide collection (of of Jan 2018), and therefore is suitable to identify a range of microorganisms, including prokaryotes and eukaryotes.
+The RefSeq_bf database contains complete reference bacterial and fungal genomes, suitable for better known habitats such as the human gut or when trying to detect well known species.
+
+**Option 2** We have indexed a more recent version of the ncbi nucleotide collection (June 2019) that does not contain environemntal or artificial sequences. The file ncbi_nt_no_env_11jun2019.zip can be found [here](http://dx.doi.org/10.25910/5cc7cd40fca8e) and contains all ncbi nt entries excluding the descendants of environmental eukaryotes (taxid 61964), environmental prokaryotes (48479), unclassified sequences (12908) and artificial sequences (28384).
+
+**Option 3:** Build your own reference database (recommended!)
+Follow the instructions in the [KMA website](https://bitbucket.org/genomicepidemiology/kma) to index the database.
+It is important that taxids are incorporated in sequence headers for processing with CCMetagen. Sequence headers should look like 
+`>1234|sequence_description`, where 1234 is the taxid. 
+We provide scripts to rename sequences in the nt database [here](https://github.com/vrmarcelino/CCMetagen/tree/master/benchmarking/rename_nt).
+
+If you want to use the RefSeq database, the format is similar to the one required for Kraken. The [Opiniomics blog](http://www.opiniomics.org/building-a-kraken-database-with-new-ftp-structure-and-no-gi-numbers/) describes how to download sequences in an adequate format. Note that you still need to build the index with KMA: `kma_index -i refseq.fna -o refseq_indexed -NI -Sparse -` or `kma_index -i refseq.fna -o refseq_indexed -NI -Sparse TG` for faster analysis.
+
+
+## Quick Start
+
+  * First map sequence reads (or contigs) to the database with **KMA**.
+
+For paired-end files:
+```
+kma -ipe $SAMPLE_R1 $SAMPLE_R2 -o sample_out_kma -t_db $db -t $th -1t1 -mem_mode -and -apm f
+```
+
+For single-end files:
+```
+kma -i $SAMPLE -o sample_out_kma -t_db $db -t $th -1t1 -mem_mode -and
+```
+
+If you want to calculate abundance in reads per million (RPM) or in number of reads (fragments), or if you want to calculate the proportion of mapped reads, add the flag -ef (extended features):
+```
+kma -ipe $SAMPLE_R1 $SAMPLE_R2 -o sample_out_kma -t_db $db -t $th -1t1 -mem_mode -and -apm f -ef
+```
+
+Where:
+
+$db is the path to the reference database
+$th is the number of threads
+$SAMPLE_R1 is the path to the mate1 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
+$SAMPLE_R2 is the path to the mate2 of a paired-end metagenome/metatranscriptome sample (fastq or fasta)
+$SAMPLE is the path to a single-end metagenome/metatranscriptome file (reads or contigs)
+
+
+  * Then run **CCMetagen**:
+```
+CCMetagen.py -i $sample_out_kma.res -o results
+```
+Where $sample_out_kma.res is alignment results produced by KMA.
+
+Note that if you are running CCMetagen from the local folder (instead of adding it to your path), you may need to add 'python' before CCMetagen: `python CCMetagen.py -i $sample_out_kma.res -o results`
+
+Done! This will make an additional quality filter and output a text file with ranked taxonomic classifications and a krona graph file for interactive visualization.
+
+An example of the CCMetagen output can be found [here (.csv file)](https://github.com/vrmarcelino/CCMetagen/blob/master/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv) and [here (.html file)](https://htmlpreview.github.io/?https://github.com/vrmarcelino/CCMetagen/blob/master/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html).
+
+<img src=docs/tutorial/figs_tutorial/krona_photo.png width="500" height="419.64">
+
+In the .csv file, you will find the depth (abundance) of each match.
+
+## Abundance units
+
+**Depth can be estimated in four ways:** by counting the number of nucleotides matching the reference sequence (use flag --depth_unit nc), by applying an additional correction for template length (default in KMA and CCMetagen), by calculating depth in Reads Per Million (RPM, use flag --depth_unit rpm), or by counting the number of fragments (i.e. number of PE reads matching to teh reference sequence, use flag --depth_unit fr). If you want RPM or fragment units, you will need to suply the .mapstats file generated with KMA (which you get when running kma with the flag '-ef').
+
+
+## Balancing sensitivity and specificity
+
+You can **adjust the stringency of the taxonomic assignments** by adjusting the minimum coverage (--coverage), the minimum abundance (--depth), and the minimum level of sequence similarity (--query_identity). Coverage is the percentage of bases in the reference sequence that is covered by the consensus sequence (your query), it can be over 100% when the consensus sequence is larger than the reference (due to insertions for example). You can also adjust the KMA settings to facilitate the identification of more distant-related taxa (see below)
+
+If you change the default depth unit, we recommend adjusting the minimum abundance (--depth) to remove taxa found in low abundance accordingly. For example, you can use -d 200 (200 nucleotides) when using --depth_unit nc, which is similar to -d 0.2 when using the default '--depth_unit kma' option. If you choose to calculate abundances in RPM, you may want to adjust the minimum abundance according to your sequence depth.
+For example, to calculate abundances in RPM, and filter out all matches with less than one read per million:
+
+```
+CCMetagen.py -i $sample_out_kma.res -o results -map $sample_out_kma.mapstat --depth_unit rpm --depth 1
+```
+
+If you would like to know the **proportion of reads mapped** to each template, run kma with the '-ef' flag. This will generate a file with the '.mapstat' extension. Then provide this file to CCMetagen (-map $sample_out_kma.mapstat) and add the flag '-ef y':
+
+```
+CCMetagen.py -i $sample_out_kma.res -o results -map $sample_out_kma.mapstat -ef y
+```
+This will filter the .mapstat file, removing the templates that did not pass CCMetagen's quality control, will add the percentage of mapped reads for each template and will output a file with extension 'stats_csv'. It will also output the overall proportion of reads mapped to these templates in the terminal. For more details about the additional columns of this file, please check [KMA's manual](https://bitbucket.org/genomicepidemiology/kma/src/master/KMAspecification.pdf).
+
+When working with highly complex environemnts for which reference databases are scarce (e.g. many soil and marine metagenomes), it is common to obtain a low proportion of classified reads, especially if the sequencing depth is low. For a more sensitive analysis, besides relaxing the CCMetatgen settings, you can adjust the KMA aligner settings, by for example: removing the `-and` and the `-apm f` flags, so that you can get a match even when the reference sequences are not significantly overrepresented or when only one of the PE reads maps to the template. Check the [KMA manual](https://bitbucket.org/genomicepidemiology/kma/src/master/KMAspecification.pdf) for more details. It can also be useful to build a customized reference database with additional genomes of organisms that are closely related to what you expect to find in your samples.
+
+## Understanding the ranked taxonomic output of CCMetagen:
+The taxonomic classifications reflect the sequence similarity between query and reference sequences, according to default or user-defined similarity thresholds. For example, if a match is 97% similar to the reference sequence, the match will not get a species-level classification. If the match is 85% similar to the reference sequence, then the species, genus and family-level classifications will be 'none'.
+Note that this is different from identifications tagged as unk_x (unknown taxa). These unknowns indicate taxa where higher-rank classifications have not been defined (according to the NCBI taxonomy database), and it is unrelated to sequence similarity.
+
+
+For a list of options to customize your analyze, type:
+```
+CCMetagen.py -h
+```
+
+  * **To get the abundance of each taxon, and/or summarize results for multiple samples, use CCMetagen_merge**:
+```
+CCMetagen_merge.py -i $CCMetagen_out
+```
+
+Where $CCMetagen_out is the folder containing the CCMetagen taxonomic classifications.
+The results must be in .csv format (default or '--mode text' output of CCMetagen), and these files **must end in ".ccm.csv"**.
+
+The flag '-t' define the taxonomic level to merge the results. The default is species-level.
+
+You can also filter out specific taxa, at any taxonomic level:
+
+Use flag -kr to keep (k) or remove (r) taxa.
+Use flag -l to set the taxonomic level for the filtering.
+Use flag -tlist to list the taxa to keep or remove (separated by comma).
+
+EX1: Filter out bacteria: `CCMetagen_merge.py -i $CCMetagen_out -kr r -l Kingdom -tlist Bacteria`
+
+EX2: Filter out bacteria and Metazoa: `CCMetagen_merge.py -i $CCMetagen_out -kr r -l Kingdom -tlist Bacteria, Metazoa`
+
+EX3: Merge results at family-level, and remove Metazoa and Viridiplantae taxa at Kingdom level:
+```
+CCMetagen_merge.py -i $CCMetagen_out -t Family -kr r -l Kingdom -tlist Metazoa,Viridiplantae -o family_table
+```
+
+For species-level filtering (where there is a space in taxa names), use quotation marks.
+Ex 4: Keep only _Escherichia coli_ and _Candida albicans_:
+```
+CCMetagen_merge.py -i 05_KMetagen/ -kr k -l Species -tlist "Escherichia coli,Candida albicans"
+```
+
+If you only have one sample, you can also use CMetagen_merge to get one line per taxa.
+
+To see all options, type:
+```
+CCMetagen_merge.py -h
+```
+This file should look like [this](https://github.com/vrmarcelino/CCMetagen/blob/master/tutorial/figs_tutorial/Bird_family_table_filtered.csv).
+
+
+* **To extract sequences of a given taxon, use CCMetagen_extract_seqs**:
+
+This script will produce a fasta file containing all reads assigned to a taxon of interest. 
+Ex: Generate a fasta file containing all sequences that mapped to the genus Eschericha:
+```
+CCMetagen_extract_seqs.py -iccm $CCMetagen_out -ifrag $sample_out_kma.frag -l Genus -t Eschericha
+```
+
+Where $CCMetagen_out is the .csv file generated with CCMetagen and $sample_out_kma.frag is the .frag file generated with KMA. The frag file needs to be decompressed: `gunzip *.frag.gz`
+
+For species-level filtering (where there is a space in taxon names), use quotation marks.
+Ex: Generate a fasta file containing all sequences that mapped to _E. coli_:
+```
+CCMetagen_extract_seqs.py -iccm $CCMetagen_out -ifrag $sample_out_kma.frag -l Species -t "Escherichia coli"
+```
+
+
+**Check out our [tutorial](https://github.com/vrmarcelino/CCMetagen/tree/master/tutorial) for an applied example of the CCMetagen pipeline.**
+
+
+
+## FAQs
+
+* Error taxid not found.
+  You probably need to update your local ETE3 database, which contains the taxids and lineage information:
+```
+python
+from ete3 import NCBITaxa
+ncbi = NCBITaxa()
+ncbi.update_taxonomy_database()
+quit()
+```
+
+* TypeError: concat() got an unexpected keyword argument 'sort'.
+  If you get this error, please update the python module pandas:
+```
+pip install pandas --upgrade --user
+```
+
+* WARNING: no NCBI's taxid found for accession [something], this match will not get taxonomic ranks
+
+  This is not an error, this is just a warning indicating that one of your query sequences matchs to a genbank record for which the NCBI taxonomic identifier (taxid) is not known. CCMetagen therefore will not be able to assign taxonomic ranks to this match, but you will still be able to see it in the output file.
+
+* KeyError: "['Superkingdom' 'Kingdom' 'Phylum' 'Class' 'Order' 'Family' .... ] not in index"
+  Make sure that the output of CCMetagen ends in '.csv'.
+
+* The results of the CCMetagen_merge.py at different taxonomic levels do not sum up.
+  As explained above, this script merges all unclassified taxa at a given taxonomic level. For example, if you have 20 matches to the genus _Candida_, but only 2 matches were classified at the species level, the output of CCMetagen_merge.py -t Species (default) will only have the abundances of two classified _Candida_ species, while the others will be merged with the "Unclassified" taxa. The output of CCMetagen_merge.py -t Genus however will contain all 20 matches. 
+  If this behaviour is undesirable, one option is to disable the similarity thresholds (use flag -off) - so that all taxonomic levels are reported regardless of their similarity to the reference sequence. Alternatively, you can cluster species at the 'Closest_match' (using the flag --tax_level Closest_match).
+
+
+## Complete option list
+
+CCMetagen:
+```
+usage: CCMetagen.py [-h] [-m MODE] -i RES_FP [-o OUTPUT_FP]
+                    [-r REFERENCE_DATABASE] [-ef EXTENDED_OUTPUT_FILE]
+                    [-du DEPTH_UNIT] [-map MAPSTAT] [-d DEPTH] [-c COVERAGE]
+                    [-q QUERY_IDENTITY] [-p PVALUE] [-st SPECIES_THRESHOLD]
+                    [-gt GENUS_THRESHOLD] [-ft FAMILY_THRESHOLD]
+                    [-ot ORDER_THRESHOLD] [-ct CLASS_THRESHOLD]
+                    [-pt PHYLUM_THRESHOLD] [-off TURN_OFF_SIM_THRESHOLDS]
+                    [--version]
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -m MODE, --mode MODE  what do you want CCMetagen to do? Valid options are
+                        'visual', 'text' or 'both': text: parses kma, filters
+                        based on quality and output a text file with taxonomic
+                        information and detailed mapping information visual:
+                        parses kma, filters based on quality and output a
+                        simplified text file and a krona html file for
+                        visualization both: outputs both text and visual file
+                        formats. Default = both
+  -i RES_FP, --res_fp RES_FP
+                        Path to the KMA result (.res file)
+  -o OUTPUT_FP, --output_fp OUTPUT_FP
+                        Path to the output file. Default = CCMetagen_out
+  -r REFERENCE_DATABASE, --reference_database REFERENCE_DATABASE
+                        Which reference database was used. Options: UNITE,
+                        RefSeq or nt. Default = nt
+  -ef EXTENDED_OUTPUT_FILE, --extended_output_file EXTENDED_OUTPUT_FILE
+                        Produce an extended output file that includes the
+                        percentage of classified reads. Options: y or n. To
+                        use this featire, you need to generate the mapstat
+                        file when required unning KMA (use flag -ef), and use
+                        it as input in CCMetagen (flag --mapstat). Default = n
+  -du DEPTH_UNIT, --depth_unit DEPTH_UNIT
+                        Desired unit for Depth(abundance) measurements.
+                        Default = kma (KMA default depth, which is the number
+                        of nucleotides overlapping each template, divided by
+                        the lengh of the template). Alternatively, you can
+                        have abundance calculated in Reads Per Million (RPM,
+                        option 'rpm'), in number of nucleotides overlaping the
+                        template (option 'nc') or in number of fragments (i.e.
+                        PE reads, option 'fr'). If you use the 'nc', 'rpm' or
+                        'fr' options, remember to change the default --depth
+                        parameter accordingly. Valid options are nc, rpm, fr
+                        and kma
+  -map MAPSTAT, --mapstat MAPSTAT
+                        Path to the mapstat file produced with KMA when using
+                        the -ef flag (.mapstat). Required when calculating
+                        abundances in RPM or in number of fragments, or when
+                        producing the extended_output_file
+  -d DEPTH, --depth DEPTH
+                        minimum sequencing depth. Default = 0.2. The unit
+                        corresponds to the one used with --depth_unit If you
+                        use --depth_unit different from the default, change
+                        this accordingly.
+  -c COVERAGE, --coverage COVERAGE
+                        Minimum coverage. Default = 20 (i.e. 20% of the
+                        reference sequence)
+  -q QUERY_IDENTITY, --query_identity QUERY_IDENTITY
+                        Minimum query identity (Phylum level). Default = 50
+  -p PVALUE, --pvalue PVALUE
+                        Minimum p-value. Default = 0.05.
+  -st SPECIES_THRESHOLD, --species_threshold SPECIES_THRESHOLD
+                        Species-level similarity threshold. Default = 98.41
+  -gt GENUS_THRESHOLD, --genus_threshold GENUS_THRESHOLD
+                        Genus-level similarity threshold. Default = 96.31
+  -ft FAMILY_THRESHOLD, --family_threshold FAMILY_THRESHOLD
+                        Family-level similarity threshold. Default = 88.51
+  -ot ORDER_THRESHOLD, --order_threshold ORDER_THRESHOLD
+                        Order-level similarity threshold. Default = 81.21
+  -ct CLASS_THRESHOLD, --class_threshold CLASS_THRESHOLD
+                        Class-level similarity threshold. Default = 80.91
+  -pt PHYLUM_THRESHOLD, --phylum_threshold PHYLUM_THRESHOLD
+                        Phylum-level similarity threshold. Default = 0 - not
+                        applied
+  -off TURN_OFF_SIM_THRESHOLDS, --turn_off_sim_thresholds TURN_OFF_SIM_THRESHOLDS
+                        Turns simularity-based filtering off. Options = y or
+                        n. Default = n
+  --version             show program's version number and exit
+ ```
+
diff --git a/tutorial/PhyloSeq_graphs_publication.R b/docs/tutorial/PhyloSeq_graphs_publication.R
similarity index 100%
rename from tutorial/PhyloSeq_graphs_publication.R
rename to docs/tutorial/PhyloSeq_graphs_publication.R
diff --git a/tutorial/README.md b/docs/tutorial/README.md
similarity index 100%
rename from tutorial/README.md
rename to docs/tutorial/README.md
diff --git a/tutorial/figs_tutorial/Bird_family_table_filtered.csv b/docs/tutorial/figs_tutorial/Bird_family_table_filtered.csv
similarity index 100%
rename from tutorial/figs_tutorial/Bird_family_table_filtered.csv
rename to docs/tutorial/figs_tutorial/Bird_family_table_filtered.csv
diff --git a/tutorial/figs_tutorial/Most_abudnant_families.png b/docs/tutorial/figs_tutorial/Most_abudnant_families.png
similarity index 100%
rename from tutorial/figs_tutorial/Most_abudnant_families.png
rename to docs/tutorial/figs_tutorial/Most_abudnant_families.png
diff --git a/tutorial/figs_tutorial/Nicer_fam_figure.png b/docs/tutorial/figs_tutorial/Nicer_fam_figure.png
similarity index 100%
rename from tutorial/figs_tutorial/Nicer_fam_figure.png
rename to docs/tutorial/figs_tutorial/Nicer_fam_figure.png
diff --git a/tutorial/figs_tutorial/Subset_bird_data.zip b/docs/tutorial/figs_tutorial/Subset_bird_data.zip
similarity index 100%
rename from tutorial/figs_tutorial/Subset_bird_data.zip
rename to docs/tutorial/figs_tutorial/Subset_bird_data.zip
diff --git a/tutorial/figs_tutorial/Superkingdom.png b/docs/tutorial/figs_tutorial/Superkingdom.png
similarity index 100%
rename from tutorial/figs_tutorial/Superkingdom.png
rename to docs/tutorial/figs_tutorial/Superkingdom.png
diff --git a/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res b/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res
similarity index 100%
rename from tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res
rename to docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res
diff --git a/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv b/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv
similarity index 100%
rename from tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv
rename to docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.csv
diff --git a/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html b/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html
similarity index 100%
rename from tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html
rename to docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.html
diff --git a/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.tsv b/docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.tsv
similarity index 100%
rename from tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.tsv
rename to docs/tutorial/figs_tutorial/Turnstone_Temperate_Flu_Ng.res.tsv
diff --git a/tutorial/figs_tutorial/krona_photo.png b/docs/tutorial/figs_tutorial/krona_photo.png
similarity index 100%
rename from tutorial/figs_tutorial/krona_photo.png
rename to docs/tutorial/figs_tutorial/krona_photo.png
diff --git a/setup.py b/setup.py
index b908cbe..924856b 100644
--- a/setup.py
+++ b/setup.py
@@ -1,3 +1,4 @@
 import setuptools
 
+# See setup.cfg for setup configuration.
 setuptools.setup()
diff --git a/tutorial/.Rhistory b/tutorial/.Rhistory
deleted file mode 100755
index 93b3103..0000000
--- a/tutorial/.Rhistory
+++ /dev/null
@@ -1,512 +0,0 @@
-plot_bar(Merged_taxa, fill = "Family")
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-genus.csv')
-taxa_names(Merged_taxa)
-list_names <- taxa_names(Merged_taxa)
-list_names
-list_names[6]
-list_names[6] <- 'Others'
-list_names
-taxa_names(Merged_taxa) <- list_names
-plot_bar(Merged_taxa, fill = "Family")
-Merged_taxa
-taxa_names(Merged_taxa)
-plot_bar(Merged_taxa)
-plot_bar(Merged_taxa, fill = "Family")
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-Taxa2Merge <- c('Unclassified',Taxa2Merge) # include string 'Unclassified'
-Taxa2Merge
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-ist_names <- taxa_names(Merged_taxa)
-list_names[6] <- 'Others'
-taxa_names(Merged_taxa) <- list_names
-plot_bar(Merged_taxa, fill = "Family")
-list_names
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families = otu_table(families_raw, taxa_are_rows = TRUE)
-families
-CCMeta_physeq = phyloseq(families, tax)
-CCMeta_physeq
-#### Analyse and Plot ####
-# plot bar graph
-# takes a while if the file is large
-#plot_bar(CCMeta_physeq, fill = "Kingdom")
-#plot_bar(CCMeta_physeq, fill = "Family")
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-Taxa2Merge <- c('Unclassified',Taxa2Merge) # include string 'Unclassified'
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-list_names <- taxa_names(Merged_taxa)
-list_names[6] <- 'Others'
-taxa_names(Merged_taxa) <- list_names
-list_names
-list_names <- taxa_names(Merged_taxa)
-list_names
-list_names[13] <- 'Others' # check depending on the number of taxa to merge
-taxa_names(Merged_taxa) <- list_names
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-list_names <- taxa_names(Merged_taxa)
-list_names[13] <- 'Others' # check depending on the number of taxa to merge
-taxa_names(Merged_taxa) <- list_names
-plot_bar(Merged_taxa, fill = "Family")
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-taxa_raw
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-tax = tax_table(taxa_raw)
-families = otu_table(families_raw, taxa_are_rows = TRUE)
-families
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-SortedOTUs
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-Taxa2Merge
-Taxa2Merge <- c('Unclassified',Taxa2Merge) # include string 'Unclassified'
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-Merged_taxa
-plot_bar(Merged_taxa, fill = "Family")
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-SortedOTUs
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-Taxa2Merge
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge) # by default, this will merge all to the taxa with highest abundance (i.e. Unclassified)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-fam_check.csv')
-taxa_names(Merged_taxa)
-taxa_names(Merged_taxa)["Trichomonadidae"]
-taxa_names(Merged_taxa)["Trichomonadidae"] <- "Others"
-list_names <- taxa_names(Merged_taxa)
-list_names
-list_names[-1]
-list_names[:-1]
-list_names[,-1]
-list_names[len(list_names)]
-list_names[length(list_names)]
-list_names[length(list_names)] <- "Others"
-list_names
-taxa_names(Merged_taxa) <- list_names
-taxa_names(Merged_taxa)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-#### Analyse and Plot ####
-# plot bar graph
-# takes a while if the file is large
-#plot_bar(CCMeta_physeq, fill = "Kingdom")
-#plot_bar(CCMeta_physeq, fill = "Family")
-# There are too many taxa to represent, so let's keep the names of the 20 most abundant ones,
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-#Taxa2Merge <- c('Unclassified',Taxa2Merge) # include string 'Unclassified'
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-# by default, this will merge all to the taxa with highest abundance, but in the plots, this is called 'NA'
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-# double-check:
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-fam_check.csv')
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[20:length(SortedOTUs)] # 20th and beyond will be merged
-#Taxa2Merge <- c('Unclassified',Taxa2Merge) # include string 'Unclassified'
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-fam_check.csv')
-SortedOTUs
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[22:length(SortedOTUs)] # 22th and beyond will be merged (so we keep 20 classified taxa, plus Unclassified and 'Others')
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-#### Import data and convert to a phyloseq object ####
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[22:length(SortedOTUs)] # 22th and beyond will be merged (so we keep 20 classified taxa, plus Unclassified and 'Others')
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[22:length(SortedOTUs)] # 22th and beyond will be merged
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-fam_check.csv')
-SortedOTUs
-Taxa2Merge
-plot_bar(Merged_taxa, fill = "Superkingdom")
-plot_bar(Merged_taxa, fill = "Order")
-plot_bar(Merged_taxa, fill = "Family")
-p = plot_bar(Merged_taxa, fill="Family")
-p$data$Sample <- factor(p$data$Sample, levels = c("1_Shelduck_Disease","2_Shelduck_Healthy",
-"3_Turnstone_Temperate_Flu_Negative","4_Turnsone_Temperate_Flu_Positive",
-"5_Temperate_Duck_Flu_positive","6_Temperate_Duck_Flu_negative","7_Outback_Duck",
-"8_Avocet_Tempereate","9_Avocet_Outback"))
-families = levels(p$data$Family)
-families
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-plot_bar(CCMeta_physeq, fill = "Kingdom")
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[22:length(SortedOTUs)] # 22th and beyond will be merged
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-write.csv(otus, file='otus-fam_check.csv')
-otus <- otu_table(Merged_taxa)
-write.csv(otus, file='otus-fam_check.csv')
-plot_bar(Merged_taxa, fill = "Superkingdom")
-plot_bar(Merged_taxa, fill = "Order")
-plot_bar(Merged_taxa, fill = "Family")
-plot_bar(Merged_taxa, fill = "Order")
-p = plot_bar(Merged_taxa, fill="Family")
-p$data$Sample <- factor(p$data$Sample, levels = c("1_Shelduck_Disease","2_Shelduck_Healthy",
-"3_Turnstone_Temperate_Flu_Negative","4_Turnsone_Temperate_Flu_Positive",
-"5_Temperate_Duck_Flu_positive","6_Temperate_Duck_Flu_negative","7_Outback_Duck",
-"8_Avocet_Tempereate","9_Avocet_Outback"))
-families = levels(p$data$Family)
-families
-getPalette = colorRampPalette(c("#ffffcc", "#800026"))( 12 )
-getPalette
-family_Palette <- c("#08519c","#2171b5","#4292c6","#6baed6","#9ecae1",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#ae017e","#999999")
-c
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594",#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#ae017e","#999999")
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594",#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#ae017e","#999999")
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#ae017e","#999999")
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#4a1486","#999999")
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#999999")
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#252525")
-p + scale_fill_manual(values=family_Palette, na.value="grey") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-p + scale_fill_manual(values=family_Palette, na.value="black") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#bdbdbd")
-p + scale_fill_manual(values=family_Palette, na.value="black") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-families = levels(p$data$Family)
-families
-p$data$Family <- factor(p$data$Family, levels = c("Bacillaceae","Bacteroidaceae",
-"Campylobacteraceae","Enterobacteriaceae","Fusobacteriaceae","Pseudomonadaceae","Streptococcaceae",
-"Apiosporaceae","Bulleribasidiaceae","Cladosporiaceae","Cystofilobasidiaceae",
-"Debaryomycetaceae","Dipodascaceae","Metschnikowiaceae","Mucoraceae","Nectriaceae",
-"Rhizopodaceae","Saccotheciaceae","Sporidiobolaceae","Trichomonadidae","NA", "Unclassified"))
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#bdbdbd")
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#bdbdbd")
-p + scale_fill_manual(values=family_Palette, na.value="black") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-families = levels(p$data$Family)
-families
-p$data$Family <- factor(p$data$Family, levels = c("Bacillaceae","Bacteroidaceae",
-"Campylobacteraceae","Enterobacteriaceae","Fusobacteriaceae","Pseudomonadaceae","Streptococcaceae",
-"Apiosporaceae","Bulleribasidiaceae","Cladosporiaceae","Cystofilobasidiaceae",
-"Debaryomycetaceae","Dipodascaceae","Metschnikowiaceae","Mucoraceae","Nectriaceae",
-"Rhizopodaceae","Saccotheciaceae","Sporidiobolaceae","Trichomonadidae","Unclassified"))
-# need 12 colors for fungi
-getPalette = colorRampPalette(c("#ffffcc", "#800026"))( 12 )
-getPalette
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#bdbdbd")
-publication_fig <- p + scale_fill_manual(values=family_Palette, na.value="black") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-publication_fig
-rm(list=ls(all=TRUE))
-#### Instal and load packages + settings ####
-# to install phyloseq:
-#source('http://bioconductor.org/biocLite.R')
-#biocLite('phyloseq')
-library("phyloseq"); packageVersion("phyloseq")
-#[1] ‘1.24.2’
-library("ggplot2"); packageVersion("ggplot2")
-#[1] ‘3.1.0’
-library (RColorBrewer); packageVersion("RColorBrewer")
-# [1] ‘1.1.2’
-# set ggplot colour theme to white
-theme_set(theme_bw())
-raw_CCMetagen_data <-read.csv("Bird_family_table.csv",check.names=FALSE)
-# separate species and taxonomy columns
-taxa_raw <- as.matrix(raw_CCMetagen_data[,10:ncol(raw_CCMetagen_data)])
-rownames(taxa_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-families_raw <- as.matrix(raw_CCMetagen_data[,1:9])
-rownames(families_raw) <- raw_CCMetagen_data[,ncol(raw_CCMetagen_data)]
-# convert to phyloseq object
-tax = tax_table(taxa_raw)
-families_depth = otu_table(families_raw, taxa_are_rows = TRUE)
-families_depth
-CCMeta_physeq = phyloseq(families_depth, tax)
-CCMeta_physeq
-#### Analyse and Plot ####
-# plot bar graph
-# takes a while if the file is large
-#plot_bar(CCMeta_physeq, fill = "Kingdom")
-#plot_bar(CCMeta_physeq, fill = "Family")
-# There are too many taxa to represent, so let's keep the names of the 20 most abundant ones (+Unclassified and 'Others')
-SortedOTUs <- names(sort(taxa_sums(CCMeta_physeq), TRUE))
-Taxa2Merge <- SortedOTUs[22:length(SortedOTUs)] # 22th and beyond will be merged
-Merged_taxa <- merge_taxa(CCMeta_physeq,Taxa2Merge)
-# by default, this will merge all to the taxa with highest abundance, but in the plots, this is called 'NA'
-plot_bar(Merged_taxa, fill = "Family") + guides(fill=guide_legend(ncol=2))
-# cross-checks:
-# otus <- otu_table(Merged_taxa)
-# write.csv(otus, file='otus-fam_check.csv') (note that NA is 'Orthomyxoviridae' in this table)
-########################
-# organise samples and families, add custom colours:
-p = plot_bar(Merged_taxa, fill="Family")
-p$data$Sample <- factor(p$data$Sample, levels = c("1_Shelduck_Disease","2_Shelduck_Healthy",
-"3_Turnstone_Temperate_Flu_Negative","4_Turnsone_Temperate_Flu_Positive",
-"5_Temperate_Duck_Flu_positive","6_Temperate_Duck_Flu_negative","7_Outback_Duck",
-"8_Avocet_Tempereate","9_Avocet_Outback"))
-# organise order of Families, keeping Fungi and Bacteria Apart:
-families = levels(p$data$Family)
-families
-p$data$Family <- factor(p$data$Family, levels = c("Bacillaceae","Bacteroidaceae",
-"Campylobacteraceae","Enterobacteriaceae","Fusobacteriaceae","Pseudomonadaceae","Streptococcaceae",
-"Apiosporaceae","Bulleribasidiaceae","Cladosporiaceae","Cystofilobasidiaceae",
-"Debaryomycetaceae","Dipodascaceae","Metschnikowiaceae","Mucoraceae","Nectriaceae",
-"Rhizopodaceae","Saccotheciaceae","Sporidiobolaceae","Trichomonadidae","Unclassified"))
-# set colors
-# need 12 colors for fungi
-getPalette = colorRampPalette(c("#ffffcc", "#800026"))( 12 )
-getPalette
-# then paste after the blues, and before the pink and grey
-# I got the blues from Color Brewer
-family_Palette <- c("#084594","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#eff3ff",
-"#FFFFCC","#F3E7BC","#E7D0AD","#DCB99E","#D0A28F","#C58B80","#B97371","#AE5C62",
-"#A24553","#972E44","#8B1735","#800026", "#7a0177","#bdbdbd")
-publication_fig <- p + scale_fill_manual(values=family_Palette, na.value="black") +
-geom_bar(aes(fill=Family), stat="identity", position="stack") +
-guides(fill=guide_legend(ncol=2))
-publication_fig
-publication_fig
-publication_fig


Test patch:
diff --git a/tests/test_ccmetagen.py b/tests/test_ccmetagen.py
new file mode 100644
index 0000000..987ec02
--- /dev/null
+++ b/tests/test_ccmetagen.py
@@ -0,0 +1,31 @@
+import unittest
+import shutil
+from pathlib import Path
+
+
+class TestCCMetagen(unittest.TestCase):
+    def test_import_ccmetagen_module(self):
+        try:
+            import ccmetagen
+        except ImportError:
+            self.fail("Failed to import ccmetagen module")
+
+    def test_ccmetagen_scripts_on_path(self):
+        scripts = ["CCMetagen.py", "CCMetagen_extract_seqs.py", "CCMetagen_merge.py"]
+        for script_name in scripts:
+            try:
+                # Check if the script is on the PATH
+                script_path = shutil.which(script_name)
+                self.assertIsNotNone(script_path, f"{script_name} not found on PATH")
+
+                # Check if the script can be executed
+                script_path = Path(script_path)
+                self.assertTrue(script_path.is_file(), f"{script_name} is not a file")
+                self.assertTrue(script_path.exists(), f"{script_name} does not exist")
+
+            except Exception as e:
+                self.fail(f"Unexpected error: {e}")
+
+
+if __name__ == "__main__":
+    unittest.main()