Repository: pachterlab/gget (language: python)

Issue:
Create gget_mutate.py



Discussion / hints:

Fix patch:
diff --git a/docs/src/en/cosmic.md b/docs/src/en/cosmic.md
index a1ea512e..bf1277cd 100644
--- a/docs/src/en/cosmic.md
+++ b/docs/src/en/cosmic.md
@@ -4,6 +4,8 @@ Search for genes, mutations, and other factors associated with cancer using the
 Return format: JSON (command-line) or data frame/CSV (Python).  
 This module was written by [@AubakirovArman](https://github.com/AubakirovArman).
 
+NOTE: License fees apply for the commercial use of COSMIC. You can read more about licensing COSMIC data [here](https://cancer.sanger.ac.uk/cosmic/license).
+
 **Positional argument**  
 `searchterm`   
 Search term, which can be a mutation, gene name (or Ensembl ID), cancer type, tumor site, study ID, PubMed ID, or sample ID, as defined using the `entity` argument. Example: 'EGFR'  
diff --git a/docs/src/en/ref.md b/docs/src/en/ref.md
index 76efc29e..9d95d49c 100644
--- a/docs/src/en/ref.md
+++ b/docs/src/en/ref.md
@@ -8,7 +8,7 @@ Return format: dictionary/JSON.
 Species for which the FTPs will be fetched in the format genus_species, e.g. homo_sapiens.  
 Supports all available vertebrate and invertebrate (plants, fungi, protists, and invertebrate metazoa) genomes from Ensembl, except bacteria.  
 Note: Not required when using flags `--list_species` or `--list_iv_species`.   
-Supported shortcuts: 'human', 'mouse'
+Supported shortcuts: 'human', 'mouse', 'human_grch37' (accesses the GRCh37 genome assembly)
 
 **Optional arguments**  
 `-w` `--which`  
diff --git a/docs/src/es/cosmic.md b/docs/src/es/cosmic.md
index 8877bcac..6e9f7634 100644
--- a/docs/src/es/cosmic.md
+++ b/docs/src/es/cosmic.md
@@ -4,6 +4,8 @@ Busque genes, mutaciones, etc. asociados con cánceres utilizando la base de dat
 Produce: Resultados en formato JSON (Terminal) o Dataframe/CSV (Python).  
 `gget cosmic` fue escrito por [@AubakirovArman](https://github.com/AubakirovArman).
 
+Se aplican tarifas de licencia para el uso comercial de COSMIC. Puede leer más sobre la concesión de licencias de datos COSMIC [aquí](https://cancer.sanger.ac.uk/cosmic/license).
+
 **Parámetro posicional**  
 `searchterm`   
 Término de búsqueda. Puede ser una mutación, un nombre de gen (o ID de Ensembl), tipo de cáncer, sitio del tumor, ID de estudio, ID de PubMed o ID de muestra, tal como se define con el argumento `entity`. Ejemplo: 'EGFR'  
diff --git a/docs/src/es/ref.md b/docs/src/es/ref.md
index a286871d..cdd6a0d5 100644
--- a/docs/src/es/ref.md
+++ b/docs/src/es/ref.md
@@ -7,7 +7,7 @@ Regresa: Resultados en formato JSON.
 `species`  
 La especie por la cual que se buscará los FTP en el formato género_especies, p. ej. homo_sapiens.  
 Nota: No se requiere cuando se llama a la bandera `--list_species`.    
-Accesos directos: 'human', 'mouse'  
+Accesos directos: 'human', 'mouse', 'human_grch37' (accede al ensamblaje del genoma GRCh37)  
 
 **Parámetros optionales**  
 `-w` `--which`  
diff --git a/gget/__init__.py b/gget/__init__.py
index 70554088..3223afbc 100644
--- a/gget/__init__.py
+++ b/gget/__init__.py
@@ -15,13 +15,9 @@
 from .gget_elm import elm
 from .gget_diamond import diamond
 from .gget_cosmic import cosmic
+from .gget_mutate import mutate
 
 import logging
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
 # Mute numexpr threads info
 logging.getLogger("numexpr").setLevel(logging.WARNING)
 
diff --git a/gget/compile.py b/gget/compile.py
index 2647ef66..26983fd3 100644
--- a/gget/compile.py
+++ b/gget/compile.py
@@ -2,19 +2,10 @@
 import subprocess
 import sys
 import platform
-import logging
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
-# Constants
 from .constants import MUSCLE_GITHUB_LINK
+from .utils import set_up_logger
+logger = set_up_logger()
 
 # Get absolute package path
 PACKAGE_PATH = os.path.abspath(os.path.dirname(__file__))
@@ -35,7 +26,7 @@ def compile_muscle():
             f"Muscle compiler currently only supports Linux and Darwin, not {platform.system()}.\n"
         )
 
-    logging.info("Compiling MUSCLE binary from source... ")
+    logger.info("Compiling MUSCLE binary from source... ")
 
     # Record current working directory
     cwd = os.getcwd()
@@ -64,14 +55,14 @@ def compile_muscle():
 
     # Run make command
     if platform.system() == "Linux":
-        logging.warning(
+        logger.warning(
             "Compiling MUSCLE requires that g++, make, sed and git are installed."
         )
     if platform.system() == "Darwin":
-        logging.warning(
+        logger.warning(
             "Compiling MUSCLE requires that gcc v11, make, sed and git are installed."
         )
-        logging.warning(
+        logger.warning(
             "Please run 'brew install gcc' to install gcc v11 if the compile fails."
         )
 
@@ -86,7 +77,7 @@ def compile_muscle():
     if process_2.wait() != 0:
         sys.exit(f"'{command2}' command returned with error {process_2.wait()}.")
 
-    logging.info("MUSCLE compiled.")
+    logger.info("MUSCLE compiled.")
 
     # Change path back to cwd
     os.chdir(cwd)
diff --git a/gget/constants.py b/gget/constants.py
index a752cd4c..cf8a8010 100644
--- a/gget/constants.py
+++ b/gget/constants.py
@@ -3,6 +3,7 @@
 # Ensembl REST API server for gget seq and info
 ENSEMBL_REST_API = "http://rest.ensembl.org/"
 ENSEMBL_FTP_URL = "http://ftp.ensembl.org/pub/"
+ENSEMBL_FTP_URL_GRCH37 = "http://ftp.ensembl.org/pub/grch37/"
 # Non-vertebrate server
 ENSEMBL_FTP_URL_NV = "http://ftp.ensemblgenomes.org/pub/"
 
@@ -38,8 +39,12 @@
 EXPRESSION_URL = "https://maayanlab.cloud/archs4/search/loadExpressionTissue.php?"
 
 # Download links for ELM database
-ELM_INSTANCES_FASTA_DOWNLOAD = "http://elm.eu.org/instances.fasta?q=*&taxon=&instance_logic="
-ELM_INSTANCES_TSV_DOWNLOAD = "http://elm.eu.org/instances.tsv?q=*&taxon=&instance_logic="
+ELM_INSTANCES_FASTA_DOWNLOAD = (
+    "http://elm.eu.org/instances.fasta?q=*&taxon=&instance_logic="
+)
+ELM_INSTANCES_TSV_DOWNLOAD = (
+    "http://elm.eu.org/instances.tsv?q=*&taxon=&instance_logic="
+)
 ELM_CLASSES_TSV_DOWNLOAD = "http://elm.eu.org/elms/elms_index.tsv"
 ELM_INTDOMAINS_TSV_DOWNLOAD = "http://elm.eu.org/interactiondomains.tsv"
 
diff --git a/gget/gget_alphafold.py b/gget/gget_alphafold.py
index 7c598a58..d9f07e67 100644
--- a/gget/gget_alphafold.py
+++ b/gget/gget_alphafold.py
@@ -31,15 +31,8 @@
 from ipywidgets import GridspecLayout
 from ipywidgets import Output
 
-import logging
-
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import set_up_logger
+logger = set_up_logger()
 
 TQDM_BAR_FORMAT = (
     "{l_bar}{bar}| {n_fmt}/{total_fmt} [elapsed: {elapsed} remaining: {remaining}]"
@@ -237,7 +230,7 @@ def alphafold(
     """
 
     if platform.system() == "Windows":
-        logging.warning(
+        logger.warning(
             "gget setup alphafold and gget alphafold are not supported on Windows OS."
         )
 
@@ -261,7 +254,7 @@ def alphafold(
     try:
         import alphafold as AlphaFold
     except ImportError:
-        logging.error(
+        logger.error(
             """
             Some third-party dependencies are missing. Please run the following command: 
             >>> gget.setup('alphafold') or $ gget setup alphafold
@@ -275,7 +268,7 @@ def alphafold(
     pdb_out, err = process.communicate()
 
     if pdb_out.decode() == "":
-        logging.error(
+        logger.error(
             """
             Some third-party dependencies are missing. Please run the following command: 
             >>> gget.setup('alphafold') or $ gget setup alphafold
@@ -285,7 +278,7 @@ def alphafold(
 
     ## Check if model parameters were downloaded
     if not os.path.exists(os.path.join(PARAMS_DIR, "params/")):
-        logging.error(
+        logger.error(
             """
             The AlphaFold model parameters are missing. Please run the following command: 
             >>> gget.setup('alphafold') or $ gget setup alphafold
@@ -294,7 +287,7 @@ def alphafold(
         return
 
     if len(os.listdir(os.path.join(PARAMS_DIR, "params/"))) < 12:
-        logging.error(
+        logger.error(
             """
             The AlphaFold model parameters are missing. Please run the following command: 
             >>> gget.setup('alphafold') or $ gget setup alphafold
@@ -348,7 +341,7 @@ def alphafold(
                 )
 
     ## Move stereo_chemical_props.txt from gget bins to Alphafold package so it can be found
-    # logging.info("Locate files containing stereochemical properties.")
+    # logger.info("Locate files containing stereochemical properties.")
     ALPHAFOLD_PATH = os.path.abspath(os.path.dirname(AlphaFold.__file__))
     os.makedirs(os.path.join(ALPHAFOLD_PATH, "common/"), exist_ok=True)
     shutil.copyfile(
@@ -358,7 +351,7 @@ def alphafold(
 
     ## Validate input sequence(s)
     if verbose:
-        logging.info(f"Validating input sequence(s).")
+        logger.info(f"Validating input sequence(s).")
 
     # Handle command line passing path to FASTA as a list
     if isinstance(sequence, list) and len(sequence) == 1:
@@ -426,17 +419,17 @@ class ModelType(enum.Enum):
     if len(seqs) == 1:
         if multimer_for_monomer:
             if verbose:
-                logging.info(
+                logger.info(
                     "Using the multimer model for a single chain, as requested."
                 )
             model_type_to_use = ModelType.MULTIMER
         else:
             if verbose:
-                logging.info("Using the single-chain (monomer) model.")
+                logger.info("Using the single-chain (monomer) model.")
             model_type_to_use = ModelType.MONOMER
     else:
         if verbose:
-            logging.info(f"Using the multimer model with {len(seqs)} sequences.")
+            logger.info(f"Using the multimer model with {len(seqs)} sequences.")
         model_type_to_use = ModelType.MULTIMER
 
     # Check whether total length exceeds limit
@@ -457,13 +450,13 @@ class ModelType(enum.Enum):
             )
 
     if total_sequence_length > MAX_VALIDATED_LENGTH:
-        logging.warning(
+        logger.warning(
             f"The accuracy of this algorithm has not been fully validated above 3000 residues, and you may experience long running times or run out of memory. Total sequence length is {total_sequence_length} residues."
         )
 
     ## Find the closest source
     if verbose:
-        logging.info(f"Finding closest source for reference database.")
+        logger.info(f"Finding closest source for reference database.")
 
     ex = futures.ThreadPoolExecutor(3)
     fs = [ex.submit(fetch, source) for source in ["", "-europe", "-asia"]]
@@ -526,7 +519,7 @@ class ModelType(enum.Enum):
     features_for_chain = {}
     raw_msa_results_for_sequence = {}
     for sequence_index, sequence in enumerate(sequences, start=1):
-        # logging.info(f"Getting MSA for sequence {sequence_index}.")
+        # logger.info(f"Getting MSA for sequence {sequence_index}.")
 
         ## Manage permissions to jackhmmer binary
         command = f"chmod 755 {JACKHMMER_BINARY_PATH}"
@@ -537,7 +530,7 @@ class ModelType(enum.Enum):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error("Giving chmod 755 permissions to jackhmmer binary failed.")
+            logger.error("Giving chmod 755 permissions to jackhmmer binary failed.")
             return
 
         # Save the target sequence in a fasta file
@@ -568,7 +561,7 @@ class ModelType(enum.Enum):
                 single_chain_msas.append(merged_msa)
                 msa_size = len(set(merged_msa.sequences))
                 if verbose:
-                    logging.info(
+                    logger.info(
                         f"{msa_size} unique sequences found in {db_name} for sequence {sequence_index}."
                     )
             elif merged_msa.sequences and db_name == "uniprot":
@@ -721,7 +714,7 @@ class ModelType(enum.Enum):
                 prot=unrelaxed_proteins[best_model_name]
             )
         else:
-            logging.warning(
+            logger.warning(
                 "\nRunning model without relaxation stage. Use flag [--relax] ('relax=True') to include AMBER relaxation."
             )
             relaxed_pdb = protein.to_pdb(unrelaxed_proteins[best_model_name])
@@ -760,7 +753,7 @@ class ModelType(enum.Enum):
     ## Plotting
     if plot:
         if verbose:
-            logging.info("Plotting prediction results.")
+            logger.info("Plotting prediction results.")
         import py3Dmol
 
         # Construct multiclass b-factors to indicate confidence bands
diff --git a/gget/gget_archs4.py b/gget/gget_archs4.py
index 0aaa49ec..aac3a85b 100644
--- a/gget/gget_archs4.py
+++ b/gget/gget_archs4.py
@@ -2,16 +2,9 @@
 import pandas as pd
 import json as json_package
 import io
-import logging
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import set_up_logger
+logger = set_up_logger()
 
 # Custom functions
 from .gget_info import info
@@ -78,7 +71,7 @@ def archs4(
 
         # Check if Ensembl ID was found
         if isinstance(info_df, type(None)):
-            logging.error(
+            logger.error(
                 f"ID '{gene}' not found. Please double-check spelling/arguments and try again."
             )
             return
@@ -96,7 +89,7 @@ def archs4(
 
     if which == "correlation":
         if verbose:
-            logging.info(
+            logger.info(
                 f"Fetching the {gene_count} most correlated genes to {gene} from ARCHS4."
             )
 
@@ -120,13 +113,13 @@ def archs4(
         # Check if the request returned an error (e.g. gene not found)
         if "error" in corr_data.keys():
             if corr_data["error"] == f"{gene} not in colids":
-                logging.error(
+                logger.error(
                     f"Gene '{gene}' did not return any gene correlation results. \n"
                     "If the gene is an Ensembl ID, please set argument 'ensembl=True' (for terminal, add flag: [--ensembl])."
                 )
                 return
             else:
-                logging.error(
+                logger.error(
                     f"Gene correlation request for search term '{gene}' returned error: {corr_data['error']}"
                 )
                 return
@@ -157,7 +150,7 @@ def archs4(
 
     if which == "tissue":
         if verbose:
-            logging.info(
+            logger.info(
                 f"Fetching the tissue expression atlas of {gene} from {species} ARCHS4 data."
             )
 
@@ -182,7 +175,7 @@ def archs4(
         tissue_exp_df = pd.read_csv(io.StringIO(r.content.decode("utf-8")))
         # Check if any results were returned
         if len(tissue_exp_df) < 2:
-            logging.error(
+            logger.error(
                 f"Gene '{gene}' did not return any tissue expression results. \n"
                 "If the gene is an Ensembl ID, please set argument 'ensembl=True' (for terminal, add flag: [--ensembl])."
             )
diff --git a/gget/gget_blast.py b/gget/gget_blast.py
index 86a306c0..19347886 100644
--- a/gget/gget_blast.py
+++ b/gget/gget_blast.py
@@ -2,16 +2,6 @@
 import json as json_package
 import time
 from bs4 import BeautifulSoup
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
 
 # Using urllib instead of requests here because requests does not
 # support long queries (queries very long here due to input sequence)
@@ -19,7 +9,8 @@
 from urllib.parse import urlencode
 
 # Custom functions
-from .utils import parse_blast_ref_page, wrap_cols_func
+from .utils import parse_blast_ref_page, wrap_cols_func, read_fasta, set_up_logger
+logger = set_up_logger()
 
 # Constants
 from .constants import (
@@ -107,28 +98,8 @@ def blast(
                         seqs.append(line.strip())
 
         elif ".fa" in sequence:
-            # Read the FASTA
-            titles = []
-            seqs = []
-            with open(sequence) as fasta_file:
-                for i, line in enumerate(fasta_file):
-                    # Each second line will be a title line
-                    if i % 2 == 0:
-                        if line[0] != ">":
-                            raise ValueError(
-                                "Expected FASTA to start with a '>' character. "
-                            )
-                        else:
-                            # Append title line to titles list
-                            titles.append(line.strip())
-                    else:
-                        if line[0] == ">":
-                            raise ValueError(
-                                "FASTA contains two lines starting with '>' in a row -> missing sequence line. "
-                            )
-                        # Append sequences line to seqs list
-                        else:
-                            seqs.append(line.strip())
+            titles, seqs = read_fasta(sequence)
+
         else:
             raise ValueError(
                 "File format not recognized. gget BLAST currently only supports '.txt' or '.fa' files. "
@@ -137,7 +108,7 @@ def blast(
         # Set the first sequence from the fasta file as 'sequence'
         sequence = seqs[0]
         if len(seqs) > 1:
-            logging.warning(
+            logger.warning(
                 "File contains more than one sequence. Only the first sequence will be submitted to BLAST."
             )
 
@@ -168,8 +139,8 @@ def blast(
             if database == "default":
                 database = "nt"
                 if verbose:
-                    logging.info("Sequence recognized as nucleotide sequence.")
-                    logging.info("BLAST will use program 'blastn' with database 'nt'.")
+                    logger.info("Sequence recognized as nucleotide sequence.")
+                    logger.info("BLAST will use program 'blastn' with database 'nt'.")
             else:
                 # Check if the user specified database is valid
                 if database not in dbs:
@@ -179,8 +150,8 @@ def blast(
 
                 else:
                     if verbose:
-                        logging.info("Sequence recognized as nucleotide sequence.")
-                        logging.info(
+                        logger.info("Sequence recognized as nucleotide sequence.")
+                        logger.info(
                             "BLAST will use program 'blastn' with user-specified database."
                         )
         # If sequence is an amino acid sequence, set program to blastp
@@ -191,8 +162,8 @@ def blast(
             if database == "default":
                 database = "nr"
                 if verbose:
-                    logging.info("Sequence recognized as amino acid sequence.")
-                    logging.info("BLAST will use program 'blastp' with database 'nr'.")
+                    logger.info("Sequence recognized as amino acid sequence.")
+                    logger.info("BLAST will use program 'blastp' with database 'nr'.")
             else:
                 # Check if the user specified database is valid
                 if database not in dbs:
@@ -202,8 +173,8 @@ def blast(
 
                 else:
                     if verbose:
-                        logging.info("Sequence recognized as amino acid sequence.")
-                        logging.info(
+                        logger.info("Sequence recognized as amino acid sequence.")
+                        logger.info(
                             "BLAST will use program 'blastp' with user-specified database."
                         )
         else:
@@ -285,12 +256,12 @@ def blast(
     if RTOE < 11:
         # Communicate RTOE
         if verbose:
-            logging.info(f"BLAST initiated. Estimated time to completion: 11 seconds.")
+            logger.info(f"BLAST initiated. Estimated time to completion: 11 seconds.")
         time.sleep(11)
     else:
         # Communicate RTOE
         if verbose:
-            logging.info(
+            logger.info(
                 f"BLAST initiated with search ID {RID}. Estimated time to completion: {RTOE} seconds."
             )
         time.sleep(int(RTOE))
@@ -329,23 +300,23 @@ def blast(
 
         if status == "WAITING":
             if verbose:
-                logging.info("BLASTING...")
+                logger.info("BLASTING...")
             i += 1
             continue
 
         elif status == "FAILED":
-            logging.error(
+            logger.error(
                 f"Search {RID} failed; please try again and/or report to blast-help@ncbi.nlm.nih.gov."
             )
             return
 
         elif status == "UNKNOWN":
-            logging.error(f"NCBI status {status}. Search {RID} expired.")
+            logger.error(f"NCBI status {status}. Search {RID} expired.")
             return
 
         elif status == "READY":
             if verbose:
-                logging.info("Retrieving results...")
+                logger.info("Retrieving results...")
             # Stop search
             searching = False
 
@@ -360,7 +331,7 @@ def blast(
             )
 
             if dsc_table is None:
-                logging.error(
+                logger.error(
                     f"No significant similarity found for search {RID}. If your sequence is very short, try increasing the 'expect' argument."
                 )
                 return
@@ -389,7 +360,7 @@ def blast(
                 return results_df
 
         else:
-            logging.error(
+            logger.error(
                 f"Something unexpected happened. Search {RID} possibly failed; please try again and/or report to blast-help@ncbi.nlm.nih.gov"
             )
             return
diff --git a/gget/gget_blat.py b/gget/gget_blat.py
index ccfe430a..f857ca25 100644
--- a/gget/gget_blat.py
+++ b/gget/gget_blat.py
@@ -1,19 +1,11 @@
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 import json as json_package
 from json.decoder import JSONDecodeError
 import pandas as pd
 from urllib.request import urlopen
 
+from .utils import set_up_logger
+logger = set_up_logger()
+
 
 def blat(sequence, seqtype="default", assembly="human", json=False, save=False, verbose=True,):
     """
@@ -82,14 +74,14 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
         sequence = seqs[0]
         if len(seqs) > 1:
             if verbose:
-                logging.info(
+                logger.info(
                     "File contains more than one sequence. Only the first sequence will be submitted to BLAT."
                 )
 
     # Shorten sequence to length limit if necessary
     if len(sequence) > 8000:
         if verbose:
-            logging.info(
+            logger.info(
                 "Length of sequence is > 8000. Only the fist 8000 characters will be submitted to BLAT."
             )
         sequence = sequence[:8000]
@@ -112,7 +104,7 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
         if set(sequence) <= nucleotides:
             seqtype = "DNA"
             if verbose:
-                logging.info(
+                logger.info(
                     f"Sequence recognized as nucleotide sequence. 'seqtype' will be set as {seqtype}."
                 )
 
@@ -120,7 +112,7 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
         elif set(sequence) <= amino_acids:
             seqtype = "protein"
             if verbose:
-                logging.info(
+                logger.info(
                     f"Sequence recognized as amino acid sequence. 'seqtype' will be set as {seqtype}."
                 )
 
@@ -166,7 +158,7 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
         # Read json results into a dictionary
         results = json_package.load(r)
     except JSONDecodeError:
-        logging.error(
+        logger.error(
             f"""
             BLAT of seqtype '{seqtype}' using assembly '{database}' was unsuccesful. 
             Possible causes: 
@@ -178,7 +170,7 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
 
     if len(results["blat"]) == 0:
         if verbose:
-            logging.info(
+            logger.info(
                 f"No {seqtype} BLAT matches were found for this sequence in genome {results['genome']}."
             )
         return
@@ -186,7 +178,7 @@ def blat(sequence, seqtype="default", assembly="human", json=False, save=False,
     # Let user know if assembly was not found
     # If this is the case, BLAT automatically defaults to human (hg38)
     if results["genome"] != database:
-        logging.warning(
+        logger.warning(
             f"Assembly {database} not recognized. Defaulted to {results['genome']} instead."
         )
 
diff --git a/gget/gget_cellxgene.py b/gget/gget_cellxgene.py
index 388f4e74..8e827a2f 100644
--- a/gget/gget_cellxgene.py
+++ b/gget/gget_cellxgene.py
@@ -1,14 +1,5 @@
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
+from .utils import set_up_logger
+logger = set_up_logger()
 
 def convert_to_list(lst):
     """
@@ -117,7 +108,7 @@ def cellxgene(
     try:
         import cellxgene_census
     except ImportError:
-        logging.error(
+        logger.error(
             """
             Some third-party dependencies are missing. Please run the following command: 
             >>> gget.setup('cellxgene') or $ gget setup cellxgene
@@ -151,7 +142,7 @@ def cellxgene(
     ]
 
     if all(el is None for el in args):
-        logging.warning(
+        logger.warning(
             """
             You are attempting to query the entire Census dataset which requires a large amount of RAM (100's of GBs) and high network bandwidth. 
             Use the cell metadata arguments (e.g. 'tissue', 'cell_type', 'disease', etc...) to define the (sub)dataset of interest.
@@ -181,7 +172,7 @@ def cellxgene(
     # Fetch AnnData object
     if not meta_only:
         if verbose:
-            logging.info(
+            logger.info(
                 "Fetching AnnData object from CZ CELLxGENE Discover. This might take a few minutes..."
             )
 
@@ -207,7 +198,7 @@ def cellxgene(
     # Fetch metadata
     else:
         if verbose:
-            logging.info("Fetching metadata from CZ CELLxGENE Discover...")
+            logger.info("Fetching metadata from CZ CELLxGENE Discover...")
         with cellxgene_census.open_soma(census_version=census_version) as census:
             # Reads SOMADataFrame as a slice
             cell_metadata = census["census_data"][species].obs.read(
diff --git a/gget/gget_cosmic.py b/gget/gget_cosmic.py
index c01fe6ce..5a6d41b3 100644
--- a/gget/gget_cosmic.py
+++ b/gget/gget_cosmic.py
@@ -1,244 +1,558 @@
 import requests
 import pandas as pd
+import subprocess
+import os
+import re
 import json as json_package
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+import base64
+import shutil
+import tarfile
+import gzip
+import getpass
 
 # Constants
 from .constants import COSMIC_GET_URL
+from .utils import set_up_logger
+
+logger = set_up_logger()
+
+
+def is_valid_email(email):
+    """
+    Check if an e-mail address is valid.
+    """
+    email_pattern = re.compile(r"(^[a-zA-Z0-9_.+-]+@[a-zA-Z0-9-]+\.[a-zA-Z0-9-.]+$)")
+
+    return re.match(email_pattern, email) is not None
+
+
+def download_reference(download_link, tar_folder_path, file_path, verbose):
+    email = input("Please enter your COSMIC email: ")
+    if not is_valid_email(email):
+        raise ValueError("The email address is not valid.")
+    password = getpass.getpass("Please enter your COSMIC password: ")
+
+    # Concatenate the email and password with a colon
+    input_string = f"{email}:{password}\n"
+
+    encoded_bytes = base64.b64encode(input_string.encode("utf-8"))
+    encoded_string = encoded_bytes.decode("utf-8")
+    curl_command = [
+        "curl",
+        "-H",
+        f"Authorization: Basic {encoded_string}",
+        download_link,
+    ]
+
+    if verbose:
+        logger.info("Downloading data...")
+
+    result = subprocess.run(curl_command, capture_output=True, text=True)
+
+    try:
+        response_data = json_package.loads(result.stdout)
+    except JSONDecodeError:
+        raise RuntimeError(
+            "Failed to download file. Please double-check arguments (especially cosmic_version) and try again."
+        )
+    try:
+        true_download_url = response_data.get("url")
+    except AttributeError:
+        raise AttributeError("Invalid username or password.")
+
+    curl_command2 = ["curl", true_download_url, "--output", f"{tar_folder_path}.tar"]
+    result2 = subprocess.run(curl_command2, capture_output=True, text=True)
+
+    if result2.returncode != 0:
+        raise RuntimeError(
+            f"Failed to download file. Return code: {result.returncode}\n{result.stderr}"
+        )
+
+    with tarfile.open(f"{tar_folder_path}.tar", "r") as tar:
+        tar.extractall(path=tar_folder_path)
+        if verbose:
+            logger.info(f"Extracted tar file to {tar_folder_path}")
+
+    with gzip.open(f"{file_path}.gz", "rb") as f_in:
+        with open(file_path, "wb") as f_out:
+            shutil.copyfileobj(f_in, f_out)
+        if verbose:
+            logger.info(f"Unzipped file to {file_path}")
+
+
+def select_reference(
+    mutation_class, reference_dir, grch_version, cosmic_version, verbose
+):
+    # if mutation_class == "transcriptome":
+    #     download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=grch{grch_version}/cosmic/v{cosmic_version}/Cosmic_Genes_Fasta_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+    #     tarred_folder = f"Cosmic_Genes_Fasta_v{cosmic_version}_GRCh{grch_version}"
+    #     contained_file = f"Cosmic_Genes_v{cosmic_version}_GRCh{grch_version}.fasta"
+
+    if mutation_class == "cancer":
+        # assert (
+        #     grch_version == 37
+        # ), "CancerMutationCensus data is only available for GRCh37 (set grch_version=37)"
+        download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=GRCh{grch_version}/cmc/v{cosmic_version}/CancerMutationCensus_AllData_Tsv_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+        tarred_folder = (
+            f"CancerMutationCensus_AllData_Tsv_v{cosmic_version}_GRCh{grch_version}"
+        )
+        contained_file = (
+            f"CancerMutationCensus_AllData_v{cosmic_version}_GRCh{grch_version}.tsv"
+        )
+
+    elif mutation_class == "cell_line":
+        download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=grch{grch_version}/cell_lines/v{cosmic_version}/CellLinesProject_GenomeScreensMutant_Tsv_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+        tarred_folder = f"CellLinesProject_GenomeScreensMutant_Tsv_v{cosmic_version}_GRCh{grch_version}"
+        contained_file = f"CellLinesProject_GenomeScreensMutant_v{cosmic_version}_GRCh{grch_version}.tsv"
+
+    elif mutation_class == "census":
+        download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=grch{grch_version}/cosmic/v{cosmic_version}/Cosmic_MutantCensus_Tsv_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+        tarred_folder = f"Cosmic_MutantCensus_Tsv_v{cosmic_version}_GRCh{grch_version}"
+        contained_file = f"Cosmic_MutantCensus_v{cosmic_version}_GRCh{grch_version}.tsv"
+
+    elif mutation_class == "resistance":
+        download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=grch{grch_version}/cosmic/v{cosmic_version}/Cosmic_ResistanceMutations_Tsv_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+        tarred_folder = (
+            f"Cosmic_ResistanceMutations_Tsv_v{cosmic_version}_GRCh{grch_version}"
+        )
+        contained_file = (
+            f"Cosmic_ResistanceMutations_v{cosmic_version}_GRCh{grch_version}.tsv"
+        )
+
+    elif mutation_class == "screen":
+        download_link = f"https://cancer.sanger.ac.uk/api/mono/products/v1/downloads/scripted?path=grch{grch_version}/cosmic/v{cosmic_version}/Cosmic_GenomeScreensMutant_Tsv_v{cosmic_version}_GRCh{grch_version}.tar&bucket=downloads"
+        tarred_folder = (
+            f"Cosmic_GenomeScreensMutant_Tsv_v{cosmic_version}_GRCh{grch_version}"
+        )
+        contained_file = (
+            f"Cosmic_GenomeScreensMutant_v{cosmic_version}_GRCh{grch_version}.tsv"
+        )
+
+    elif mutation_class == "cancer_example":
+        download_link = f"https://cog.sanger.ac.uk/cosmic-downloads-production/taster/example_grch{grch_version}.tar"
+        tarred_folder = f"example_GRCh{grch_version}"
+        contained_file = (
+            f"CancerMutationCensus_AllData_v{cosmic_version}_GRCh{grch_version}.tsv"
+        )
+
+    tar_folder_path = os.path.join(reference_dir, tarred_folder)
+    file_path = os.path.join(tar_folder_path, contained_file)
+
+    overwrite = True
+    if os.path.exists(file_path):
+        proceed = (
+            input(
+                "The requested COSMIC database already exists at the destination. Would you like to overwrite the existing files (y/n)? "
+            )
+            .strip()
+            .lower()
+        )
+        if proceed in ["yes", "y"]:
+            overwrite = True
+        else:
+            overwrite = False
+
+    if overwrite:
+        # Only the example database can be downloaded directly (without an account)
+        if mutation_class == "cancer_example":
+            curl_command = [
+                "curl",
+                "-L",
+                "--output",
+                f"{tar_folder_path}.tar",
+                download_link,
+            ]
+            result = subprocess.run(curl_command, capture_output=True, text=True)
+
+            with tarfile.open(f"{tar_folder_path}.tar", "r") as tar:
+                tar.extractall(path=tar_folder_path)
+                if verbose:
+                    logger.info(f"Extracted tar file to {tar_folder_path}")
+
+        # Download full databases
+        else:
+            proceed = (
+                input(
+                    "Downloading complete databases from COSMIC requires an account (https://cancer.sanger.ac.uk/cosmic/register; free for academic use, license for commercial use).\nWould you like to proceed (y/n)? "
+                )
+                .strip()
+                .lower()
+            )
+            if proceed in ["yes", "y"]:
+                download_reference(download_link, tar_folder_path, file_path, verbose)
+            else:
+                raise KeyboardInterrupt(
+                    f"Database download canceled. Learn more about COSMIC at https://cancer.sanger.ac.uk/cosmic/download/cosmic."
+                )
+
+    return file_path, overwrite
 
 
 def cosmic(
-    searchterm, entity="mutations", limit=100, save=False, verbose=True, json=False
+    searchterm,
+    entity="mutations",
+    limit=100,
+    json=False,
+    download_cosmic=False,
+    mutation_class="cancer",
+    cosmic_version=99,
+    grch_version=37,
+    gget_mutate=True,
+    out=None,
+    verbose=True,
 ):
     """
     Search for genes, mutations, etc associated with cancers using the COSMIC
     (Catalogue Of Somatic Mutations In Cancer) database
     (https://cancer.sanger.ac.uk/cosmic).
+    NOTE: Licence fees apply for the commercial use of COSMIC.
 
-    Args:
-    - searchterm    (str) Search term, which can be a mutation, or gene name (or Ensembl ID), or sample, etc.
-                    as defined using the 'entity' argument. Example: 'EGFR'.
-    - entity        (str) Defines the type of the supplied search term. One of the following:
-                    'mutations' (default), 'genes', 'cancer', 'tumour site', 'studies', 'pubmed', or 'samples'.
-    - limit         (int) Number of hits to return. Default: 100
-    - json          (True/False) If True, returns results in json format instead of data frame. Default: False
-    - save          (True/False) whether to save the results in the local directory. Default: False
-    - verbose       (True/False) whether to print progress information. Default: True
+    Args for querying information about specific cancers/genes/etc:
+    - searchterm      (str) Search term, which can be a mutation, or gene name (or Ensembl ID), or sample, etc.
+                      as defined using the 'entity' argument. Example: 'EGFR'.
+                      Set to None when downloading COSMIC databases with download_cosmic=True.
+    - entity          (str) Defines the type of the supplied search term. One of the following:
+                      'mutations' (default), 'genes', 'cancer', 'tumour_site', 'studies', 'pubmed', or 'samples'.
+    - limit           (int) Number of hits to return. Default: 100
+    - json            (True/False) If True, returns results in json format instead of data frame. Default: False
 
     Returns a data frame with the requested results.
+
+    Args for downloading COSMIC databases:
+    NOTE: Downloading complete databases from COSMIC requires an account (https://cancer.sanger.ac.uk/cosmic/register; free for academic use, license for commercial use)
+    - download_cosmic (True/False) whether to switch into database download mode. Default: False
+    - mutation_class  (str) Type of COSMIC database to download. One of the following:
+                      'cancer' (default), 'cell_line', 'census', 'resistance', 'screen', 'cancer_example'
+    - cosmic_version  (int) Version of the COSMIC database. Default: 99
+    - grch_version    (int) Version of the human GRCh reference genome the COSMIC database was based on (37 or 38). Default: 37
+    - gget_mutate     (True/False) whether to create a modified version of the database for use with gget mutate. Default: True
+
+    General args:
+    - out             (str) Path to folder the database will be downloaded into.
+                      Default: None
+                      -> When download_cosmic=False: Results will be returned to standard out
+                      -> When download_cosmic=True: Database will be downloaded into current working directory
+    - verbose         (True/False) whether to print progress information. Default: True
+
+    Saves the requested database into the specified folder (or current working directory if out=None).
     """
 
     if verbose:
-        logging.info("NOTE: Licence fees apply for the commercial use of COSMIC.")
-
-    # Check if 'entity' argument is valid
-    sps = [
-        "mutations",
-        "pubmed",
-        "genes",
-        "studies",
-        "samples",
-        "cancer",
-        "tumour site",
-    ]
-    if entity not in sps:
-        raise ValueError(
-            f"'entity' argument specified as {entity}. Expected one of: {', '.join(sps)}"
-        )
+        logger.info("NOTE: Licence fees apply for the commercial use of COSMIC.")
 
-    # Translate categories to match COSMIC data table IDs
-    if entity == "cancer":
-        entity = "disease"
+    ## Database download
+    if download_cosmic:
+        mut_class_allowed = [
+            "cancer",
+            "cell_line",
+            "census",
+            "resistance",
+            "screen",
+            "cancer_example",
+        ]
+        if mutation_class not in mut_class_allowed:
+            raise ValueError(
+                f"Parameter 'mutation_class' must be one of the following: {', '.join(mut_class_allowed)}.\n"
+            )
 
-    if entity == "tumour site":
-        entity = "tumour"
+        grch_allowed = [37, 38]
+        if grch_version not in grch_allowed:
+            raise ValueError(
+                f"Parameter 'grch_version' must be one of the following: {', '.join(grch_allowed)}.\n"
+            )
 
-    if verbose:
-        logging.info(
-            f"Fetching the {limit} most correlated to {searchterm} from COSMIC."
-        )
+        if not out:
+            out = os.path.dirname(os.getcwd())
 
-    r = requests.get(url=COSMIC_GET_URL + entity + "?q=" + searchterm + "&export=json")
+        if not os.path.exists(out):
+            os.makedirs(out)
 
-    # Check if the request returned an error (e.g. gene not found)
-    if not r.ok:
-        raise RuntimeError(
-            f"COSMIC API request returned error {r.status_code}. "
-            "Please double-check the arguments and try again.\n"
+        ## Download requested database
+        mutation_tsv_file, overwrite = select_reference(
+            mutation_class, out, grch_version, cosmic_version, verbose
         )
 
-    if r.text == "\n":
-        logging.warning(
-            f"searchterm = '{searchterm}' did not return any results with entity = '{entity}'. "
-            "Please double-check the arguments and try again.\n"
-        )
-        return None
-
-    data = r.text.split("\n")
-    dicts = {}
-    counter = 1
-    if entity == "mutations":
-        dicts = {"Gene": [], "Syntax": [], "Alternate IDs": [], "Canonical": []}
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Gene"].append(parsing_mutations[0])
-                dicts["Syntax"].append(parsing_mutations[1])
-                dicts["Alternate IDs"].append(
-                    parsing_mutations[2].replace('" ', "").replace('"', "")
-                )
-                dicts["Canonical"].append(parsing_mutations[3])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "pubmed":
-        dicts = {"Pubmed": [], "Paper title": [], "Author": []}
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Pubmed"].append(parsing_mutations[0])
-                dicts["Paper title"].append(
-                    parsing_mutations[1].replace('" ', "").replace('"', "").capitalize()
-                )
-                dicts["Author"].append(
-                    parsing_mutations[2].replace('" ', "").replace('"', "")
+        if gget_mutate and overwrite:
+            ## Create copy of results formatted for further use by gget mutate
+            if verbose:
+                logger.info(
+                    "Creating modified mutations file for use with gget mutate..."
                 )
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "genes":
-        dicts = {
-            "Gene": [],
-            "Alternate IDs": [],
-            "Tested samples": [],
-            "Simple Mutations": [],
-            "Fusions": [],
-            "Coding Mutations": [],
-        }
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Gene"].append(parsing_mutations[0])
-                dicts["Alternate IDs"].append(
-                    parsing_mutations[1].replace('" ', "").replace('"', "")
-                )
-                dicts["Tested samples"].append(parsing_mutations[2])
-                dicts["Simple Mutations"].append(parsing_mutations[3])
-                dicts["Fusions"].append(parsing_mutations[4])
-                dicts["Coding Mutations"].append(parsing_mutations[5])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "samples":
-        dicts = {
-            "Sample Name": [],
-            "Sites & Histologies": [],
-            "Analysed Genes": [],
-            "Mutations": [],
-            "Fusions": [],
-            "Structual variants": [],
-        }
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Sample Name"].append(parsing_mutations[0])
-                dicts["Sites & Histologies"].append(
-                    parsing_mutations[1].replace(":", ", ")
+
+            if mutation_class == "cancer" or mutation_class == "cancer_example":
+                relevant_cols = [
+                    "ACCESSION_NUMBER",
+                    "GENOMIC_MUTATION_ID",
+                    "MUTATION_URL",
+                    "Mutation CDS",
+                ]
+            else:
+                relevant_cols = [
+                    "TRANSCRIPT_ACCESSION",
+                    "GENOMIC_MUTATION_ID",
+                    "MUTATION_ID",
+                    "MUTATION_CDS",
+                ]
+
+            df = pd.read_csv(mutation_tsv_file, usecols=relevant_cols, sep="\t")
+
+            # Get seq_ID and mutation columns
+            if mutation_class == "cancer" or mutation_class == "cancer_example":
+                df["MUTATION_URL"] = df["MUTATION_URL"].str.extract(r"id=(\d+)")
+                df = df.rename(
+                    columns={
+                        "ACCESSION_NUMBER": "seq_ID",
+                        "Mutation CDS": "mutation",
+                        "MUTATION_URL": "MUTATION_ID",
+                    }
                 )
-                dicts["Analysed Genes"].append(parsing_mutations[2])
-                dicts["Mutations"].append(parsing_mutations[3])
-                dicts["Fusions"].append(parsing_mutations[4])
-                dicts["Structual variants"].append(parsing_mutations[5])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "studies":
-        dicts = {
-            "Study Id": [],
-            "Project Code": [],
-            "Description": [],
-        }
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Study Id"].append(parsing_mutations[0])
-                dicts["Project Code"].append(parsing_mutations[1])
-                dicts["Description"].append(parsing_mutations[2])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "disease":
-        dicts = {
-            "COSMIC classification": [],
-            "Paper description": [],
-            "Tested samples": [],
-            "Mutations": [],
-        }
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["COSMIC classification"].append(
-                    parsing_mutations[0].replace('" ', "").replace('"', "")
+            else:
+                df = df.rename(
+                    columns={
+                        "TRANSCRIPT_ACCESSION": "seq_ID",
+                        "MUTATION_CDS": "mutation",
+                    }
                 )
-                dicts["Paper description"].append(
-                    parsing_mutations[1].replace('" ', "").replace('"', "")
+
+            # Remove version numbers from Ensembl IDs
+            df["seq_ID"] = df["seq_ID"].str.split(".").str[0]
+
+            # Get mut_ID column (by combining GENOMIC_MUTATION_ID and MUTATION_URL/MUTATION_ID)
+            df["GENOMIC_MUTATION_ID"] = df["GENOMIC_MUTATION_ID"].fillna("NA")
+            df["GENOMIC_MUTATION_ID"] = df["GENOMIC_MUTATION_ID"].astype(str)
+            df["MUTATION_ID"] = df["MUTATION_ID"].astype(str)
+            df["mut_ID"] = df["GENOMIC_MUTATION_ID"] + "_" + df["MUTATION_ID"]
+            df = df.drop(columns=["GENOMIC_MUTATION_ID", "MUTATION_ID"])
+
+            mutate_csv_out = mutation_tsv_file.replace(".tsv", "_gget_mutate.csv")
+            df.to_csv(mutate_csv_out, index=False)
+
+            if verbose:
+                logger.info(
+                    f"Modified mutations file for use with gget mutate created at {mutate_csv_out}"
                 )
-                dicts["Tested samples"].append(parsing_mutations[2])
-                dicts["Mutations"].append(parsing_mutations[3])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    elif entity == "tumour":
-        dicts = {
-            "Primary Site": [],
-            "Tested sample": [],
-            "Analyzed genes": [],
-            "Mutations": [],
-            "Fusions": [],
-            "Structural variants": [],
-        }
-        for i in data:
-            if len(i) > 2:
-                parsing_mutations = i.split("\t")
-                dicts["Primary Site"].append(parsing_mutations[0])
-                dicts["Tested sample"].append(parsing_mutations[1])
-                dicts["Analyzed genes"].append(parsing_mutations[2])
-                dicts["Mutations"].append(parsing_mutations[3])
-                dicts["Fusions"].append(parsing_mutations[4])
-                dicts["Structural variants"].append(parsing_mutations[5])
-                counter = counter + 1
-                if limit < counter:
-                    break
-
-    corr_df = pd.DataFrame(dicts)
-
-    if json:
-        results_dict = json_package.loads(corr_df.to_json(orient="records"))
-        if save:
-            with open(
-                f"gget_cosmic_{entity}_{searchterm}.json", "w", encoding="utf-8"
-            ) as f:
-                json_package.dump(results_dict, f, ensure_ascii=False, indent=4)
-
-        return results_dict
 
     else:
-        if save:
-            corr_df.to_csv(f"gget_cosmic_{entity}_{searchterm}.csv", index=False)
+        # Check if 'entity' argument is valid
+        sps = [
+            "mutations",
+            "pubmed",
+            "genes",
+            "studies",
+            "samples",
+            "cancer",
+            "tumour_site",
+        ]
+        if entity not in sps:
+            raise ValueError(
+                f"'entity' argument specified as {entity}. Expected one of: {', '.join(sps)}"
+            )
+
+        # Translate categories to match COSMIC data table IDs
+        if entity == "cancer":
+            entity = "disease"
+
+        if entity == "tumour_site":
+            entity = "tumour"
+
+        if verbose:
+            logger.info(
+                f"Fetching the {limit} most correlated to {searchterm} from COSMIC."
+            )
+
+        r = requests.get(
+            url=COSMIC_GET_URL + entity + "?q=" + searchterm + "&export=json"
+        )
+
+        # Check if the request returned an error (e.g. gene not found)
+        if not r.ok:
+            raise RuntimeError(
+                f"COSMIC API request returned error {r.status_code}. "
+                "Please double-check the arguments and try again.\n"
+            )
+
+        if r.text == "\n":
+            logger.warning(
+                f"searchterm = '{searchterm}' did not return any results with entity = '{entity}'. "
+                "Please double-check the arguments and try again.\n"
+            )
+            return None
+
+        data = r.text.split("\n")
+        dicts = {}
+        counter = 1
+        if entity == "mutations":
+            dicts = {"Gene": [], "Syntax": [], "Alternate IDs": [], "Canonical": []}
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Gene"].append(parsing_mutations[0])
+                    dicts["Syntax"].append(parsing_mutations[1])
+                    dicts["Alternate IDs"].append(
+                        parsing_mutations[2].replace('" ', "").replace('"', "")
+                    )
+                    dicts["Canonical"].append(parsing_mutations[3])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "pubmed":
+            dicts = {"Pubmed": [], "Paper title": [], "Author": []}
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Pubmed"].append(parsing_mutations[0])
+                    dicts["Paper title"].append(
+                        parsing_mutations[1]
+                        .replace('" ', "")
+                        .replace('"', "")
+                        .capitalize()
+                    )
+                    dicts["Author"].append(
+                        parsing_mutations[2].replace('" ', "").replace('"', "")
+                    )
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "genes":
+            dicts = {
+                "Gene": [],
+                "Alternate IDs": [],
+                "Tested samples": [],
+                "Simple Mutations": [],
+                "Fusions": [],
+                "Coding Mutations": [],
+            }
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Gene"].append(parsing_mutations[0])
+                    dicts["Alternate IDs"].append(
+                        parsing_mutations[1].replace('" ', "").replace('"', "")
+                    )
+                    dicts["Tested samples"].append(parsing_mutations[2])
+                    dicts["Simple Mutations"].append(parsing_mutations[3])
+                    dicts["Fusions"].append(parsing_mutations[4])
+                    dicts["Coding Mutations"].append(parsing_mutations[5])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "samples":
+            dicts = {
+                "Sample Name": [],
+                "Sites & Histologies": [],
+                "Analysed Genes": [],
+                "Mutations": [],
+                "Fusions": [],
+                "Structual variants": [],
+            }
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Sample Name"].append(parsing_mutations[0])
+                    dicts["Sites & Histologies"].append(
+                        parsing_mutations[1].replace(":", ", ")
+                    )
+                    dicts["Analysed Genes"].append(parsing_mutations[2])
+                    dicts["Mutations"].append(parsing_mutations[3])
+                    dicts["Fusions"].append(parsing_mutations[4])
+                    dicts["Structual variants"].append(parsing_mutations[5])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "studies":
+            dicts = {
+                "Study Id": [],
+                "Project Code": [],
+                "Description": [],
+            }
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Study Id"].append(parsing_mutations[0])
+                    dicts["Project Code"].append(parsing_mutations[1])
+                    dicts["Description"].append(parsing_mutations[2])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "disease":
+            dicts = {
+                "COSMIC classification": [],
+                "Paper description": [],
+                "Tested samples": [],
+                "Mutations": [],
+            }
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["COSMIC classification"].append(
+                        parsing_mutations[0].replace('" ', "").replace('"', "")
+                    )
+                    dicts["Paper description"].append(
+                        parsing_mutations[1].replace('" ', "").replace('"', "")
+                    )
+                    dicts["Tested samples"].append(parsing_mutations[2])
+                    dicts["Mutations"].append(parsing_mutations[3])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        elif entity == "tumour":
+            dicts = {
+                "Primary Site": [],
+                "Tested sample": [],
+                "Analyzed genes": [],
+                "Mutations": [],
+                "Fusions": [],
+                "Structural variants": [],
+            }
+            for i in data:
+                if len(i) > 2:
+                    parsing_mutations = i.split("\t")
+                    dicts["Primary Site"].append(parsing_mutations[0])
+                    dicts["Tested sample"].append(parsing_mutations[1])
+                    dicts["Analyzed genes"].append(parsing_mutations[2])
+                    dicts["Mutations"].append(parsing_mutations[3])
+                    dicts["Fusions"].append(parsing_mutations[4])
+                    dicts["Structural variants"].append(parsing_mutations[5])
+                    counter = counter + 1
+                    if limit < counter:
+                        break
+
+        corr_df = pd.DataFrame(dicts)
+
+        json_out = os.path.join(
+            out, f"Cosmic_v{cosmic_version}_GRCh{grch_version}_clean.csv"
+        )
+
+        if json:
+            results_dict = json_package.loads(corr_df.to_json(orient="records"))
+            if out:
+                # Create saving directory
+                directory = "/".join(out.split("/")[:-1])
+                if directory != "":
+                    os.makedirs(directory, exist_ok=True)
+
+                json_out = os.path.join(out, f"gget_cosmic_{entity}_{searchterm}.json")
+                with open(json_out, "w", encoding="utf-8") as f:
+                    json_package.dump(results_dict, f, ensure_ascii=False, indent=4)
+
+            else:
+                return results_dict
+
+        else:
+            if out:
+                # Create saving directory
+                directory = "/".join(out.split("/")[:-1])
+                if directory != "":
+                    os.makedirs(directory, exist_ok=True)
+
+                df_out = os.path.join(out, f"gget_cosmic_{entity}_{searchterm}.csv")
+                corr_df.to_csv(df_out, index=False)
 
-        return corr_df
+            else:
+                return corr_df
diff --git a/gget/gget_diamond.py b/gget/gget_diamond.py
index 24f760a1..196a2ac3 100644
--- a/gget/gget_diamond.py
+++ b/gget/gget_diamond.py
@@ -1,4 +1,3 @@
-import logging
 import subprocess
 import sys
 import platform
@@ -7,17 +6,9 @@
 import uuid
 import json as json_package
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 from .compile import PACKAGE_PATH
-from .utils import tsv_to_df, create_tmp_fasta, remove_temp_files
+from .utils import tsv_to_df, create_tmp_fasta, remove_temp_files, set_up_logger
+logger = set_up_logger()
 
 # Path to precompiled diamond binary
 if platform.system() == "Windows":
@@ -138,7 +129,7 @@ def diamond(
 
     # Run DIAMOND
     if verbose:
-        logging.info(f"Creating DIAMOND database and initiating alignment...")
+        logger.info(f"Creating DIAMOND database and initiating alignment...")
 
     with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process:
         stderr = process.stderr.read().decode("utf-8")
@@ -151,7 +142,7 @@ def diamond(
         raise RuntimeError("DIAMOND alignment failed.")
     else:
         if verbose:
-            logging.info(f"DIAMOND alignment complete.")
+            logger.info(f"DIAMOND alignment complete.")
 
     df_diamond = tsv_to_df(
         output,
diff --git a/gget/gget_elm.py b/gget/gget_elm.py
index 0c6bf665..df5a362a 100644
--- a/gget/gget_elm.py
+++ b/gget/gget_elm.py
@@ -1,20 +1,12 @@
 import pandas as pd
 import numpy as np
 import os
-import logging
 import json as json_package
 import re
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import get_uniprot_seqs, tsv_to_df, set_up_logger
+logger = set_up_logger()
 
-from .utils import get_uniprot_seqs, tsv_to_df
 from .constants import UNIPROT_REST_API
 from .gget_diamond import diamond
 from .gget_setup import (
@@ -114,7 +106,7 @@ def seq_workflow(
     Returns: data frame consisting of ELM instances, class information, start, end in query, and if motif overlaps with subject sequence
     """
     if verbose:
-        logging.info(
+        logger.info(
             f"ORTHO Performing pairwise sequence alignment against ELM database using DIAMOND for {len(sequences)} sequence(s)..."
         )
 
@@ -132,7 +124,7 @@ def seq_workflow(
         )
 
         if len(df_diamond) == 0:
-            logging.warning(
+            logger.warning(
                 f"ORTHO Sequence {seq_number}/{len(sequences)}: No orthologous proteins found in ELM database."
             )
 
@@ -144,7 +136,7 @@ def seq_workflow(
                     str(id).split("|")[1]
                     for id in df_diamond["subject_accession"].values
                 ]
-                logging.info(
+                logger.info(
                     f"ORTHO Sequence {seq_number}/{len(sequences)}: DIAMOND found the following orthologous proteins: {', '.join(uniprot_ids)}. Retrieving ELMs for each UniProt Acc..."
                 )
 
@@ -291,11 +283,11 @@ def elm(
     with open(ELM_CLASSES_TSV) as input_file:
         head = [next(input_file) for _ in range(lines_number)]
     if verbose:
-        logging.info(", ".join(head).replace("#", "").replace("\n", ""))
+        logger.info(", ".join(head).replace("#", "").replace("\n", ""))
     with open(ELM_INSTANCES_TSV) as input_file:
         head = [next(input_file) for _ in range(lines_number)]
     if verbose:
-        logging.info(", ".join(head).replace("#", "").replace("\n", ""))
+        logger.info(", ".join(head).replace("#", "").replace("\n", ""))
 
     # Check validity of amino acid seq
     if not uniprot:
@@ -305,19 +297,19 @@ def elm(
 
         # If sequence is not a valid amino sequence, raise error
         if not set(sequence) <= amino_acids:
-            logging.warning(
+            logger.warning(
                 f"Input amino acid sequence contains invalid characters. If the input is a UniProt Acc, please use flag --uniprot (Python: uniprot=True)."
             )
 
     # Build ortholog dataframe
     if verbose:
-        logging.info(f"ORTHO Compiling ortholog information...")
+        logger.info(f"ORTHO Compiling ortholog information...")
     ortho_df = pd.DataFrame()
     if uniprot:
         ortho_df = get_elm_instances(sequence)
 
         if len(ortho_df) == 0:
-            logging.warning(
+            logger.warning(
                 "ORTHO The provided UniProt Accession does not match UniProt Accessions in the ELM database. Fetching amino acid sequence from UniProt..."
             )
             df_uniprot = get_uniprot_seqs(server=UNIPROT_REST_API, ensembl_ids=sequence)
@@ -357,7 +349,7 @@ def elm(
         )
 
         if len(ortho_df) == 0:
-            logging.warning(
+            logger.warning(
                 "ORTHO No ELM database orthologs found for input sequence or UniProt Acc."
             )
 
@@ -408,7 +400,7 @@ def elm(
 
     # Build data frame containing regex motif matches
     if verbose:
-        logging.info(f"REGEX Finding regex motif matches...")
+        logger.info(f"REGEX Finding regex motif matches...")
     fetch_aa_failed = False
     if uniprot:
         # use amino acid sequence associated with UniProt Acc to do regex match
@@ -424,13 +416,13 @@ def elm(
             ].values
 
             if len(sequences) == 0:
-                logging.warning(
+                logger.warning(
                     f"REGEX No amino acid sequences found for UniProt Acc {sequence} from the UniProt server."
                 )
                 fetch_aa_failed = True
             else:
                 if len(sequences) > 1:
-                    logging.warning(
+                    logger.warning(
                         f"REGEX More than one amino acid sequence found for UniProt Acc {sequence}. Using best match to find regex motifs."
                     )
                 sequence = sequences[0]
@@ -440,7 +432,7 @@ def elm(
         df_regex_matches = regex_match(sequence)
 
     if len(df_regex_matches) == 0:
-        logging.warning(
+        logger.warning(
             "REGEX No regex matches found for input sequence or UniProt Acc."
         )
 
diff --git a/gget/gget_enrichr.py b/gget/gget_enrichr.py
index 0eea3db3..69562a56 100644
--- a/gget/gget_enrichr.py
+++ b/gget/gget_enrichr.py
@@ -2,16 +2,6 @@
 import pandas as pd
 import json as json_package
 import numpy as np
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
 
 # Plotting packages
 import matplotlib.pyplot as plt
@@ -27,6 +17,9 @@
 from .compile import PACKAGE_PATH
 from .gget_info import info
 
+from .utils import set_up_logger
+logger = set_up_logger()
+
 def ensembl_to_gene_names(ensembl_ids):
     """
     Function to fetch gene names from a list of Ensembl IDs using gget info.
@@ -41,7 +34,7 @@ def ensembl_to_gene_names(ensembl_ids):
 
         # Check if Ensembl ID was found
         if isinstance(info_df, type(None)):
-            logging.warning(
+            logger.warning(
                 f"ID '{gene_id}' not found. Please double-check spelling/arguments."
             )
             continue
@@ -128,61 +121,61 @@ def enrichr(
     if database == "pathway":
         database = "KEGG_2021_Human"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     elif database == "transcription":
         database = "ChEA_2016"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     elif database == "ontology":
         database = "GO_Biological_Process_2021"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     elif database == "diseases_drugs":
         database = "GWAS_Catalog_2019"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     elif database == "celltypes":
         database = "PanglaoDB_Augmented_2021"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     elif database == "kinase_interactions":
         database = "KEA_2015"
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using database {database}. " + db_message
             )
 
     else:
         database = database
         if verbose:
-            logging.info(f"Performing Enichr analysis using database {database}.")
+            logger.info(f"Performing Enichr analysis using database {database}.")
 
     # To generate a KEGG pathway image, confirm that the database is a KEGG database and pykegg is installed
     if kegg_out:
         if not database.startswith("KEGG"):
-            logging.error(
+            logger.error(
                 "Please specify a KEGG database when generating a KEGG pathway image."
             )
             return
         try:
             import pykegg
         except ImportError:
-            logging.error(
+            logger.error(
                 "Please install `pykegg` to generate a KEGG pathway image. Pykegg can be installed using pip: 'pip install pykegg'"
             )
             return
@@ -194,7 +187,7 @@ def enrichr(
     ## Transform Ensembl IDs to gene symbols
     if ensembl:
         if verbose:
-            logging.info("Getting gene symbols from Ensembl IDs.")
+            logger.info("Getting gene symbols from Ensembl IDs.")
 
         genes_v2 = ensembl_to_gene_names(genes)
 
@@ -202,7 +195,7 @@ def enrichr(
         genes_v2 = genes
 
     if len(genes_v2) == 0 and ensembl:
-        logging.error("No gene symbols found for given Ensembl IDs.")
+        logger.error("No gene symbols found for given Ensembl IDs.")
         return
 
     # Transform Ensembl IDs to gene symbols for background genes
@@ -211,14 +204,14 @@ def enrichr(
 
     if not isinstance(background_list, type(None)):
         if len(background_list) == 0 and ensembl:
-            logging.error("No background gene symbols found for given Ensembl IDs.")
+            logger.error("No background gene symbols found for given Ensembl IDs.")
             return
 
     genes_clean = clean_genes_list(genes_v2)
 
     if ensembl:
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis on the following gene symbols: {', '.join(genes_clean)}"
             )
 
@@ -253,19 +246,19 @@ def enrichr(
     # If user gives a background list, use the user input instead of the default
     if background_list:
         if verbose:
-            logging.info(
+            logger.info(
                 f"Performing Enichr analysis using user-defined background gene list."
             )
 
         if background:
-            logging.warning(
+            logger.warning(
                 "Since you provided a list of background genes, the 'background==True' argument to use the default background gene list is being ignored."
             )
         background_final = "\n".join(background_list)
 
     elif background:
         if verbose:
-            logging.info(
+            logger.info(
                 "Background genes set to > 20,000 default background genes listed here: https://github.com/pachterlab/gget/blob/main/gget/constants/enrichr_bkg_genes.txt."
             )
         with open(f"{PACKAGE_PATH}/constants/enrichr_bkg_genes.txt") as f:
@@ -342,7 +335,7 @@ def enrichr(
         df = pd.DataFrame(enrichr_results[database], columns=columns)
 
     except KeyError:
-        logging.error(
+        logger.error(
             f"""
             Database {database} not found. Go to https://maayanlab.cloud/Enrichr/#libraries 
             for a full list of supported databases.
@@ -357,7 +350,7 @@ def enrichr(
     df["database"] = database
 
     if len(df) == 0:
-        logging.error(
+        logger.error(
             f"""
             No Enrichr results were found for genes {genes_clean} and database {database}. \n
             If the genes are Ensembl IDs, please set argument 'ensembl=True' (for terminal, add flag: [--ensembl]).
diff --git a/gget/gget_gpt.py b/gget/gget_gpt.py
index d1d4b5f3..1a7e7d5d 100644
--- a/gget/gget_gpt.py
+++ b/gget/gget_gpt.py
@@ -1,13 +1,5 @@
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import set_up_logger
+logger = set_up_logger()
 
 
 def gpt(
@@ -61,7 +53,7 @@ def gpt(
     try:
         import openai
     except ImportError:
-        logging.error(
+        logger.error(
             """
             Some third-party dependencies are missing. Please run the following command: 
             >>> gget.setup('gpt') or $ gget setup gpt
@@ -107,7 +99,7 @@ def gpt(
         )
 
     if verbose:
-        logging.info(
+        logger.info(
             f"Total tokens used for API call to model '{model}': {response['usage']['total_tokens']}"
         )
 
diff --git a/gget/gget_info.py b/gget/gget_info.py
index 045365a3..629bd8a4 100644
--- a/gget/gget_info.py
+++ b/gget/gget_info.py
@@ -4,20 +4,9 @@
 import requests
 from bs4 import BeautifulSoup
 
-# import json
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 # Custom functions
-from .utils import rest_query, get_uniprot_info, wrap_cols_func, get_pdb_ids
+from .utils import rest_query, get_uniprot_info, wrap_cols_func, get_pdb_ids, set_up_logger
+logger = set_up_logger()
 
 # Constants
 from .constants import ENSEMBL_REST_API, UNIPROT_REST_API, NCBI_URL
@@ -59,12 +48,12 @@ def info(
     # Handle deprecated arguments
     if expand:
         if verbose:
-            logging.warning(
+            logger.warning(
                 "'expand' argument deprecated! gget info now always returns all of the available information."
             )
     if ensembl_only:
         if verbose:
-            logging.warning(
+            logger.warning(
                 "'ensembl_only' argument deprecated! Please use arguments 'ncbi=False' and 'uniprot=False'."
             )
 
@@ -98,7 +87,7 @@ def info(
 
             if "." in ensembl_ID and temp == 0:
                 if verbose is True:
-                    logging.info(
+                    logger.info(
                         "We noticed that you passed a version number with your Ensembl ID.\n"
                         "Please note that gget info will always return information linked to the latest Ensembl ID version (see 'ensembl_id')."
                     )
@@ -154,7 +143,7 @@ def info(
             # Log error if this also did not work
             except RuntimeError:
                 if verbose:
-                    logging.warning(
+                    logger.warning(
                         f"ID '{ensembl_ID}' not found. Please double-check spelling/arguments and try again."
                     )
 
@@ -197,7 +186,7 @@ def info(
 
                 except Exception as e:
                     if verbose:
-                        logging.warning(
+                        logger.warning(
                             f"UniProt server request for ID '{ens_id}' return following error:\n{e}"
                         )
                     continue
@@ -210,7 +199,7 @@ def info(
                         df_uniprot = df_uniprot.iloc[[0]]
 
                         if verbose:
-                            logging.warning(
+                            logger.warning(
                                 f"More than one UniProt match was found for ID {ens_id}. Only the first match and its associated information will be returned."
                             )
 
@@ -226,7 +215,7 @@ def info(
 
                 else:
                     if verbose:
-                        logging.warning(f"No UniProt entry was found for ID {ens_id}.")
+                        logger.warning(f"No UniProt entry was found for ID {ens_id}.")
 
             if fetch_ncbi is True:
                 ## Get NCBI gene ID and description (for genes only)
@@ -236,7 +225,7 @@ def info(
                     html = requests.get(url)
                     # Raise error if status code not "OK" Response
                     if html.status_code != 200:
-                        logging.error(
+                        logger.error(
                             f"NCBI server request for {ens_id} returned error status code:\n{html.status_code}.\nPlease double-check arguments or try again later."
                         )
 
@@ -253,7 +242,7 @@ def info(
                             "li", class_="error icon"
                         ).text.strip()
 
-                        logging.error(
+                        logger.error(
                             f"The NCBI server request for Ensembl ID '{ens_id}' returned the following error:\n{error_message}"
                         )
 
@@ -295,7 +284,7 @@ def info(
                         ncbi_synonyms = None
 
                 except Exception as e:
-                    logging.error(
+                    logger.error(
                         f"The NCBI server request for Ensembl ID '{ens_id}' returned the following error:\n{e}"
                     )
 
@@ -323,7 +312,7 @@ def info(
 
                 except Exception as e:
                     if verbose:
-                        logging.warning(
+                        logger.warning(
                             f"The PDBe server request for Ensembl ID '{ens_id}' returned the following error:\n{e}"
                         )
                     continue
diff --git a/gget/gget_muscle.py b/gget/gget_muscle.py
index b0dd0014..61553635 100644
--- a/gget/gget_muscle.py
+++ b/gget/gget_muscle.py
@@ -4,21 +4,12 @@
 import itertools
 import sys
 import time
-import logging
 import uuid
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 # Custom functions
 from .compile import compile_muscle, MUSCLE_PATH, PACKAGE_PATH
-from .utils import aa_colors, n_colors, create_tmp_fasta
+from .utils import aa_colors, n_colors, create_tmp_fasta, set_up_logger
+logger = set_up_logger()
 
 # Path to precompiled muscle binary
 if platform.system() == "Windows":
@@ -78,7 +69,7 @@ def muscle(fasta, super5=False, out=None, verbose=True):
 
     else:
         if verbose:
-            logging.info("MUSCLE compiled. ")
+            logger.info("MUSCLE compiled. ")
         muscle_path = PRECOMPILED_MUSCLE_PATH
 
     # Replace slashes in path for Windows compatibility
@@ -118,7 +109,7 @@ def muscle(fasta, super5=False, out=None, verbose=True):
     # Record MUSCLE align start
     start_time = time.time()
     if verbose:
-        logging.info("MUSCLE aligning... ")
+        logger.info("MUSCLE aligning... ")
 
     # Run muscle command and write command output
     with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process_2:
@@ -131,7 +122,7 @@ def muscle(fasta, super5=False, out=None, verbose=True):
         return
     else:
         if verbose:
-            logging.info(
+            logger.info(
                 f"MUSCLE alignment complete. Alignment time: {round(time.time() - start_time, 2)} seconds"
             )
 
diff --git a/gget/gget_mutate.py b/gget/gget_mutate.py
new file mode 100644
index 00000000..8871c232
--- /dev/null
+++ b/gget/gget_mutate.py
@@ -0,0 +1,757 @@
+import pandas as pd
+import re
+from tqdm import tqdm
+
+tqdm.pandas()
+
+from .utils import read_fasta, set_up_logger
+
+logger = set_up_logger()
+
+# Define global variables to count occurences of weird mutations
+intronic_mutations = 0
+posttranslational_region_mutations = 0
+unknown_mutations = 0
+uncertain_mutations = 0
+ambiguous_position_mutations = 0
+cosmic_incorrect_wt_base = 0
+mut_idx_outside_seq = 0
+
+mutation_pattern = r"c\.([0-9_\-\+\*]+)([a-zA-Z>]+)"  # more complex: r'c\.([0-9_\-\+\*\(\)\?]+)([a-zA-Z>\(\)0-9]+)'
+
+
+def extract_mutation_type(mutation):
+    match = re.match(mutation_pattern, mutation)
+
+    if match:
+        letters = match.group(2)  # The letter sequence
+        mutation_type = re.findall(r"[a-z>-]+", letters)
+        if mutation_type[0] == ">":  # substitution - eg c.1252C>T
+            return "substitution"
+
+        elif mutation_type[0] == "del":  # deletion - eg c.2126_2128del OR c.1627del
+            return "deletion"
+
+        elif (
+            mutation_type[0] == "delins"
+        ):  # insertion-deletion - eg c.2239_2253delinsAAT OR c.646delins
+            return "delins"
+
+        elif (
+            mutation_type[0] == "ins"
+        ):  # insertion - eg c.2239_2240insAAT   (interval is always 1)
+            return "insertion"
+
+        elif mutation_type[0] == "dup":  # duplication - eg c.2239_2253dup OR c.646dup
+            return "duplication"
+
+        elif (
+            mutation_type[0] == "inv"
+        ):  # inversion - eg c.2239_2253inv    (always has an interval)
+            return "inversion"
+
+        else:
+            logger.debug(
+                f"mutation {mutation} matches re but is not a known mutation type"
+            )
+            # raise_pytest_error()
+            return "unknown"
+    else:
+        logger.debug(f"mutation {mutation} does not match re")
+        # raise_pytest_error()
+        return "unknown"
+
+
+def create_mutant_sequence(
+    row, mutation_function, kmer_flanking_length, mut_column, seq_id_column
+):
+    global intronic_mutations, posttranslational_region_mutations, unknown_mutations, uncertain_mutations, ambiguous_position_mutations
+
+    if "?" in row[mut_column]:
+        logger.debug(
+            f"Uncertain mutation found in mutation {row[mut_column]} - mutation is ignored"
+        )
+        uncertain_mutations += 1
+        return ""
+
+    if "(" in row[mut_column]:
+        logger.debug(
+            f"Ambiguous mutational position found in mutation {row[mut_column]} - mutation is ignored"
+        )
+        ambiguous_position_mutations += 1
+        return ""
+
+    match = re.match(mutation_pattern, row[mut_column])
+    if match:
+        nucleotide_position = match.group(1)  # The number sequence
+        letters = match.group(2)  # The letter sequence
+
+        try:
+            if "-" in nucleotide_position or "+" in nucleotide_position:
+                logger.debug(
+                    f"Intronic nucleotide position found in mutation {row[mut_column]} - {nucleotide_position}{letters} - mutation is ignored"
+                )
+                intronic_mutations += 1
+                return ""
+
+            elif "*" in nucleotide_position:
+                logger.debug(
+                    f"Posttranslational region nucleotide position found in mutation {row[mut_column]} - {nucleotide_position}{letters} - mutation is ignored"
+                )
+                posttranslational_region_mutations += 1
+                return ""
+
+            else:
+                if "_" in nucleotide_position:
+                    starting_nucleotide_position_index_0 = (
+                        int(nucleotide_position.split("_")[0]) - 1
+                    )
+                    ending_nucleotide_position_index_0 = (
+                        int(nucleotide_position.split("_")[1]) - 1
+                    )
+                else:
+                    starting_nucleotide_position_index_0 = int(nucleotide_position) - 1
+                    ending_nucleotide_position_index_0 = (
+                        starting_nucleotide_position_index_0
+                    )
+        except:
+            logger.debug(f"Error with mutation {row[mut_column]}]")
+            unknown_mutations += 1
+            # raise_pytest_error()
+            return ""
+
+        mutant_sequence, adjusted_end_position = mutation_function(
+            row,
+            letters,
+            starting_nucleotide_position_index_0,
+            ending_nucleotide_position_index_0,
+            mut_column,
+            seq_id_column,
+        )
+
+        if not mutant_sequence:
+            return ""
+
+        kmer_start = max(0, starting_nucleotide_position_index_0 - kmer_flanking_length)
+        kmer_end = min(
+            len(mutant_sequence), adjusted_end_position + kmer_flanking_length + 1
+        )
+        return str(mutant_sequence[kmer_start:kmer_end])
+    else:
+        logger.debug(f"Error with mutation {row[mut_column]}")
+        unknown_mutations += 1
+        # raise_pytest_error()
+        return ""
+
+
+def substitution_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    # assert letters[0] == row['full_sequence'][starting_nucleotide_position_index_0], f"Transcript has {row['full_sequence'][starting_nucleotide_position_index_0]} at position {starting_nucleotide_position_index_0} but mutation is {letters[-1]} at position {starting_nucleotide_position_index_0} in {row[mut_column]}"
+    global cosmic_incorrect_wt_base
+    global mut_idx_outside_seq
+    try:
+        if letters[0] != row["full_sequence"][starting_nucleotide_position_index_0]:
+            logger.debug(
+                f"Sequence {row[seq_id_column]} has nucleotide '{row['full_sequence'][starting_nucleotide_position_index_0]}' at position {starting_nucleotide_position_index_0}, but mutation {row[mut_column]} expected '{letters[0]}'."
+            )
+            cosmic_incorrect_wt_base += 1
+        mutant_sequence = (
+            row["full_sequence"][:starting_nucleotide_position_index_0]
+            + letters[-1]
+            + row["full_sequence"][ending_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = ending_nucleotide_position_index_0
+
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def deletion_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    global mut_idx_outside_seq
+    try:
+        # Accessing the sequence to trigger IndexError if out of bounds
+        _ = row["full_sequence"][starting_nucleotide_position_index_0]
+
+        mutant_sequence = (
+            row["full_sequence"][:starting_nucleotide_position_index_0]
+            + row["full_sequence"][ending_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = starting_nucleotide_position_index_0 - 1
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def delins_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    global mut_idx_outside_seq
+    try:
+        # Accessing the sequence to trigger IndexError if out of bounds
+        _ = row["full_sequence"][starting_nucleotide_position_index_0]
+
+        # yields 60+(length of insertion)-mer
+        insertion_string = "".join(re.findall(r"[A-Z]+", letters))
+        mutant_sequence = (
+            row["full_sequence"][:starting_nucleotide_position_index_0]
+            + insertion_string
+            + row["full_sequence"][ending_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = (
+            starting_nucleotide_position_index_0 + len(insertion_string) - 1
+        )
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def insertion_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    global mut_idx_outside_seq
+    try:
+        # Accessing the sequence to trigger IndexError if out of bounds
+        _ = row["full_sequence"][starting_nucleotide_position_index_0]
+
+        # yields 61+(length of insertion)-mer  - k before, length of insertion, k-1 after (k before and k-1 after is a little unconventional, but if I want to make it k-1 before then I have to return an adjusted start as well as an adjusted end)
+        insertion_string = "".join(re.findall(r"[A-Z]+", letters))
+        mutant_sequence = (
+            row["full_sequence"][: starting_nucleotide_position_index_0 + 1]
+            + insertion_string
+            + row["full_sequence"][starting_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = (
+            ending_nucleotide_position_index_0 + len(insertion_string) - 1
+        )
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def duplication_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    global mut_idx_outside_seq
+    try:
+        # Accessing the sequence to trigger IndexError if out of bounds
+        _ = row["full_sequence"][starting_nucleotide_position_index_0]
+        insertion_string = row["full_sequence"][
+            starting_nucleotide_position_index_0 : ending_nucleotide_position_index_0
+            + 1
+        ]
+        mutant_sequence = (
+            row["full_sequence"][:starting_nucleotide_position_index_0]
+            + insertion_string * 2
+            + row["full_sequence"][ending_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = ending_nucleotide_position_index_0 + len(
+            insertion_string
+        )
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def inversion_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    global mut_idx_outside_seq
+    try:
+        # Accessing the sequence to trigger IndexError if out of bounds
+        _ = row["full_sequence"][starting_nucleotide_position_index_0]
+        insertion_string = row["full_sequence"][
+            starting_nucleotide_position_index_0 : ending_nucleotide_position_index_0
+            + 1
+        ]
+        mutant_sequence = (
+            row["full_sequence"][:starting_nucleotide_position_index_0]
+            + insertion_string[::-1]
+            + row["full_sequence"][ending_nucleotide_position_index_0 + 1 :]
+        )
+        adjusted_end_position = ending_nucleotide_position_index_0
+        return mutant_sequence, adjusted_end_position
+
+    # Mutation index is outside sequence length
+    except IndexError:
+        mut_idx_outside_seq += 1
+        return "", ""
+
+
+def unknown_mutation(
+    row,
+    letters,
+    starting_nucleotide_position_index_0,
+    ending_nucleotide_position_index_0,
+    mut_column,
+    seq_id_column,
+):
+    return "", ""
+
+
+def mutate(
+    sequences,
+    mutations,
+    k=31,
+    mut_column="mutation",
+    mut_id_column="mut_ID",
+    seq_id_column="seq_ID",
+    out="gget_mutate_out.fa",
+    verbose=True,
+):
+    """
+    Takes in nucleotide sequences and mutations (in standard mutation annotation - see below)
+    and returns mutated versions of the input sequences according to the provided mutations.
+
+    Args:
+    - sequences     (str) Path to the fasta file containing the sequences to be mutated, e.g., 'seqs.fa'.
+                    Sequence identifiers following the '>' character must correspond to the identifiers
+                    in the seq_ID column of 'mutations'.
+                    NOTE: Only string until first space or dot will be used as sequence identifier
+                    - Version numbers of Ensembl IDs will be ignored.
+
+                    Example:
+                    >seq1 (or ENSG00000106443)
+                    ACTGCGATAGACT
+                    >seq2
+                    AGATCGCTAG
+
+                    Alternatively: Input sequence(s) as a string or list, e.g. 'AGCTAGCT' or 'ACTGCTAGCT' 'AGCTAGCT'.
+    - mutations     Path to csv or tsv file (str) (e.g., 'mutations.csv') or data frame (DataFrame object),
+                    containing information about the mutations in the following format:
+
+                    | mutation         | mut_ID | seq_ID |
+                    | c.2C>T           | mut1   | seq1   |
+                    | c.9_13inv        | mut2   | seq2   |
+                    | c.9_13inv        | mut2   | seq4   |
+                    | c.9_13delinsAAT  | mut3   | seq4   |
+                    | ...              | ...    | ...    |
+
+                    'mutation' = Column containing the mutations to be performed written in standard mutation annotation (see below)
+                    'mut_ID' = Column containing an identifier for each mutation
+                    'seq_ID' = Column containing the identifiers of the sequences to be mutated (must correspond to the string following
+                    the > character in the 'sequences' fasta file; do NOT include spaces or dots)
+
+                    Alternatively: Input mutation(s) as a string or list, e.g., 'c.2C>T' or ['c.2C>T', 'c.1A>C'].
+                    If a list is passed, the number of mutations must equal the number of input sequences.
+    - k             (int) Length of sequences flanking the mutation. Default: 31.
+                    If k > total length of the sequence, the entire sequence will be kept.
+    - mut_column    (str) Name of the column containing the mutations to be performed in 'mutations'. Default: 'mutation'.
+    - mut_id_column (str) Name of the column containing the IDs of each mutation in 'mutations'. Default: 'mut_ID'.
+    - seq_id_column (str) Name of the column containing the IDs of the sequences to be mutated in 'mutations'. Default: 'seq_ID'.
+    - out           (str) Path to output fasta file containing the mutated sequences. Default: 'gget_mutate_out.fa'.
+                    Set out=None to return a list of the mutated sequences to standard out.
+                    The identifiers (following the '>') of the mutated sequences in the output fasta will be '>[seq_ID]_[mut_ID]'.
+    - verbose       (True/False) whether to print progress information. Default: True
+
+    For more information on the standard mutation annotation, see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1867422/.
+
+    Saves mutated sequences in fasta format ( or returns a list containing the mutated sequences if out=None).
+    """
+
+    global intronic_mutations, posttranslational_region_mutations, unknown_mutations, uncertain_mutations, ambiguous_position_mutations, cosmic_incorrect_wt_base, mut_idx_outside_seq
+
+    # Load input sequences and their identifiers from fasta file
+    if "." in sequences:
+        titles, seqs = read_fasta(sequences)
+
+    # Handle input sequences passed as a list
+    elif isinstance(sequences, list):
+        titles = [f"seq{i+1}" for i in range(len(sequences))]
+        seqs = sequences
+
+    # Handle a single sequence passed as a string
+    elif isinstance(sequences, str):
+        titles = ["seq1"]
+        seqs = [sequences]
+
+    else:
+        raise ValueError(
+            """
+            Format of the input to the 'sequences' argument not recognized. 
+            'sequences' must be one of the following:
+            - Path to the fasta file containing the sequences to be mutated (e.g. 'seqs.fa')
+            - A list of sequences to be mutated (e.g. ['ACTGCTAGCT', 'AGCTAGCT'])
+            - A single sequence to be mutated passed as a string (e.g. 'AGCTAGCT')
+            """
+        )
+
+    # Read in 'mutations' if passed as filepath to comma-separated csv
+    if isinstance(mutations, str) and ".csv" in mutations:
+        mutations = pd.read_csv(mutations)
+
+    elif isinstance(mutations, str) and ".tsv" in mutations:
+        mutations = pd.read_csv(mutations, sep="\t")
+
+    # Handle mutations passed as a list
+    elif isinstance(mutations, list):
+        if len(mutations) > 1:
+            if len(mutations) != len(seqs):
+                raise ValueError(
+                    "If a list is passed, the number of mutations must equal the number of input sequences."
+                )
+
+            temp = pd.DataFrame()
+            temp["mutation"] = mutations
+            temp["mut_ID"] = [f"mut{i+1}" for i in range(len(mutations))]
+            temp["seq_ID"] = [f"seq{i+1}" for i in range(len(mutations))]
+            mutations = temp
+        else:
+            temp = pd.DataFrame()
+            temp["mutation"] = [mutations[0]] * len(seqs)
+            temp["mut_ID"] = [f"mut{i+1}" for i in range(len(seqs))]
+            temp["seq_ID"] = [f"seq{i+1}" for i in range(len(seqs))]
+            mutations = temp
+
+    # Handle single mutation passed as a string
+    elif isinstance(mutations, str):
+        # This will work for one mutation for one sequence as well as one mutation for multiple sequences
+        temp = pd.DataFrame()
+        temp["mutation"] = [mutations] * len(seqs)
+        temp["mut_ID"] = [f"mut{i+1}" for i in range(len(seqs))]
+        temp["seq_ID"] = [f"seq{i+1}" for i in range(len(seqs))]
+        mutations = temp
+
+    elif isinstance(mutations, pd.DataFrame):
+        pass
+
+    else:
+        raise ValueError(
+            """
+            Format of the input to the 'mutations' argument not recognized. 
+            'mutations' must be one of the following:
+            - Path to comma-separated csv file (e.g. 'mutations.csv')
+            - A pandas DataFrame object
+            - A single mutation to be applied to all input sequences (e.g. 'c.2C>T')
+            - A list of mutations (the number of mutations must equal the number of input sequences) (e.g. ['c.2C>T', 'c.1A>C'])
+            """
+        )
+
+    seq_dict = {}
+    for title, seq in zip(titles, seqs):
+        # Keep text following the > until the first space/dot as the sequence identifier
+        # Dots are removed so Ensembl version numbers are removed
+        seq_dict[title.split(" ")[0].split(".")[0]] = seq
+
+    # Get all mutation types
+    if verbose:
+        tqdm.pandas(desc="Extracting mutation types")
+        mutations["mutation_type"] = mutations[mut_column].progress_apply(
+            extract_mutation_type
+        )
+    else:
+        mutations["mutation_type"] = mutations[mut_column].apply(extract_mutation_type)
+
+    # Link sequences to their mutations using the sequence identifiers
+    mutations["full_sequence"] = mutations[seq_id_column].map(seq_dict)
+
+    # Handle sequences that were not found based on their sequence IDs
+    seqs_not_found = mutations[mutations["full_sequence"].isnull()]
+    if 0 < len(seqs_not_found) < 20:
+        logger.warning(
+            f"""
+            The sequences with the following {len(seqs_not_found)} sequence ID(s) were not found: {", ".join(seqs_not_found["seq_ID"].values)}  
+            These sequences and their corresponding mutations will not be included in the output.  
+            Ensure that the sequence IDs correspond to the string following the > character in the 'sequences' fasta file (do NOT include spaces or dots).
+            """
+        )
+    elif len(seqs_not_found) > 0:
+        logger.warning(
+            f"""
+            The sequences corresponding to {len(seqs_not_found)} sequence IDs were not found.  
+            These sequences and their corresponding mutations will not be included in the output.  
+            Ensure that the sequence IDs correspond to the string following the > character in the 'sequences' fasta file (do NOT include spaces or dots).
+            """
+        )
+
+    # Drop inputs for sequences that were not found
+    mutations = mutations.dropna()
+    if len(mutations) < 1:
+        raise ValueError(
+            """
+            None of the input sequences match the sequence IDs provided in 'mutations'. 
+            Ensure that the sequence IDs correspond to the string following the > character in the 'sequences' fasta file (do NOT include spaces or dots).
+            """
+        )
+
+    # Split data frame by mutation type
+    mutation_types = [
+        "substitution",
+        "deletion",
+        "delins",
+        "insertion",
+        "duplication",
+        "inversion",
+        "unknown",
+    ]
+    mutation_dict = {mutation_type: [] for mutation_type in mutation_types}
+    for mutation_type in mutation_types:
+        df_mutation_type = mutations[mutations["mutation_type"] == mutation_type].copy()
+        df_mutation_type["mutant_sequence_kmer"] = ""
+        # Define header for mutated sequences in output fasta
+        df_mutation_type["header"] = (
+            ">"
+            + df_mutation_type[seq_id_column]
+            + "_"
+            + df_mutation_type[mut_id_column]
+        )
+        mutation_dict[mutation_type] = df_mutation_type
+
+    # Create mutated sequences
+    if verbose:
+        logger.info("Mutating sequences...")
+    if not mutation_dict["substitution"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing substitutions")
+            mutation_dict["substitution"]["mutant_sequence_kmer"] = mutation_dict[
+                "substitution"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(substitution_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["substitution"]["mutant_sequence_kmer"] = mutation_dict[
+                "substitution"
+            ].apply(
+                create_mutant_sequence,
+                args=(substitution_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["deletion"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing deletions")
+            mutation_dict["deletion"]["mutant_sequence_kmer"] = mutation_dict[
+                "deletion"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(deletion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["deletion"]["mutant_sequence_kmer"] = mutation_dict[
+                "deletion"
+            ].apply(
+                create_mutant_sequence,
+                args=(deletion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["delins"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing delins")
+            mutation_dict["delins"]["mutant_sequence_kmer"] = mutation_dict[
+                "delins"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(delins_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["delins"]["mutant_sequence_kmer"] = mutation_dict[
+                "delins"
+            ].apply(
+                create_mutant_sequence,
+                args=(delins_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["insertion"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing insertions")
+            mutation_dict["insertion"]["mutant_sequence_kmer"] = mutation_dict[
+                "insertion"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(insertion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["insertion"]["mutant_sequence_kmer"] = mutation_dict[
+                "insertion"
+            ].apply(
+                create_mutant_sequence,
+                args=(insertion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["duplication"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing duplications")
+            mutation_dict["duplication"]["mutant_sequence_kmer"] = mutation_dict[
+                "duplication"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(duplication_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["duplication"]["mutant_sequence_kmer"] = mutation_dict[
+                "duplication"
+            ].apply(
+                create_mutant_sequence,
+                args=(duplication_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["inversion"].empty:
+        if verbose:
+            tqdm.pandas(desc="Performing inversions")
+            mutation_dict["inversion"]["mutant_sequence_kmer"] = mutation_dict[
+                "inversion"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(inversion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["inversion"]["mutant_sequence_kmer"] = mutation_dict[
+                "inversion"
+            ].apply(
+                create_mutant_sequence,
+                args=(inversion_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+    if not mutation_dict["unknown"].empty:
+        if verbose:
+            tqdm.pandas(desc="Unknown mutations")
+            mutation_dict["unknown"]["mutant_sequence_kmer"] = mutation_dict[
+                "unknown"
+            ].progress_apply(
+                create_mutant_sequence,
+                args=(unknown_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+        else:
+            mutation_dict["unknown"]["mutant_sequence_kmer"] = mutation_dict[
+                "unknown"
+            ].apply(
+                create_mutant_sequence,
+                args=(unknown_mutation, k, mut_column, seq_id_column),
+                axis=1,
+            )
+
+    # Report status of mutations back to user
+    total_mutations = mutations.shape[0]
+    good_mutations = (
+        total_mutations
+        - intronic_mutations
+        - posttranslational_region_mutations
+        - unknown_mutations
+        - uncertain_mutations
+        - ambiguous_position_mutations
+        - cosmic_incorrect_wt_base
+        - mut_idx_outside_seq
+    )
+
+    if good_mutations != total_mutations:
+        logger.warning(
+            f"""
+            {good_mutations} mutations correctly recorded ({good_mutations/total_mutations*100:.2f}%)
+            {intronic_mutations} intronic mutations found ({intronic_mutations/total_mutations*100:.2f}%)
+            {posttranslational_region_mutations} posttranslational region mutations found ({posttranslational_region_mutations/total_mutations*100:.2f}%)
+            {unknown_mutations} unknown mutations found ({unknown_mutations/total_mutations*100:.2f}%)
+            {uncertain_mutations} mutations with uncertain mutation found ({uncertain_mutations/total_mutations*100:.2f}%)
+            {ambiguous_position_mutations} mutations with ambiguous position found ({ambiguous_position_mutations/total_mutations*100:.2f}%)
+            {cosmic_incorrect_wt_base} mutations with incorrect wildtype base found ({cosmic_incorrect_wt_base/total_mutations*100:.2f}%)
+            {mut_idx_outside_seq} mutations with indeces outside of the sequence length found ({mut_idx_outside_seq/total_mutations*100:.2f}%)
+            """
+        )
+    else:
+        if verbose:
+            logger.info(
+                f"""
+                {good_mutations} mutations correctly recorded ({good_mutations/total_mutations*100:.2f}%)
+                {intronic_mutations} intronic mutations found ({intronic_mutations/total_mutations*100:.2f}%)
+                {posttranslational_region_mutations} posttranslational region mutations found ({posttranslational_region_mutations/total_mutations*100:.2f}%)
+                {unknown_mutations} unknown mutations found ({unknown_mutations/total_mutations*100:.2f}%)
+                {uncertain_mutations} mutations with uncertain mutation found ({uncertain_mutations/total_mutations*100:.2f}%)
+                {ambiguous_position_mutations} mutations with ambiguous position found ({ambiguous_position_mutations/total_mutations*100:.2f}%)
+                {cosmic_incorrect_wt_base} mutations with incorrect wildtype base found ({cosmic_incorrect_wt_base/total_mutations*100:.2f}%)
+                {mut_idx_outside_seq} mutations with indeces outside of the sequence length found ({mut_idx_outside_seq/total_mutations*100:.2f}%)
+                """
+            )
+
+    if out:
+        # Save mutated sequences in new fasta file
+        with open(out, "w") as fasta_file:
+            for mutation in mutation_dict:
+                if not mutation_dict[mutation].empty:
+                    df = mutation_dict[mutation]
+                    df["fasta_format"] = (
+                        df["header"] + "\n" + df["mutant_sequence_kmer"] + "\n"
+                    )
+
+                    # Do not include an empty string as a mutated sequence
+                    # (these are introduced when unknown mutations are encountered)
+                    df = df[df["mutant_sequence_kmer"] != ""]
+
+                    fasta_file.write("".join(df["fasta_format"].values))
+
+        # with open(out, "r") as fasta_file:
+        #     lines = [line for line in fasta_file if line.strip()]
+
+        # with open(out, "w") as fasta_file:
+        #     fasta_file.writelines(lines)
+
+        if verbose:
+            logger.info(f"FASTA file containing mutated sequences created at {out}.")
+
+    # When out=None, return list of mutated seqs
+    else:
+        all_mut_seqs = []
+        for mutation in mutation_dict:
+            df = mutation_dict[mutation]
+            all_mut_seqs.extend(df["mutant_sequence_kmer"].values)
+
+        # Remove empty strings from final list of mutated sequences
+        # (these are introduced when unknown mutations are encountered)
+        while "" in all_mut_seqs:
+            all_mut_seqs.remove("")
+
+        if len(all_mut_seqs) > 0:
+            return all_mut_seqs
diff --git a/gget/gget_pdb.py b/gget/gget_pdb.py
index 4b0a5a0d..d93f8a94 100644
--- a/gget/gget_pdb.py
+++ b/gget/gget_pdb.py
@@ -1,16 +1,9 @@
 from urllib.request import urlopen
 from urllib.error import HTTPError
 import json
-import logging
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import set_up_logger
+logger = set_up_logger()
 
 from .constants import RCSB_PDB_API
 
@@ -97,19 +90,19 @@ def pdb(pdb_id, resource="pdb", identifier=None, save=False):
         r = urlopen(url)
     except HTTPError:
         if resource == "assembly":
-            logging.error(
+            logger.error(
                 f"{resource} for {pdb_id} assembly {identifier} was not found. Please double-check arguments and try again."
             )
         elif resource in need_entitiy_id:
-            logging.error(
+            logger.error(
                 f"{resource} for {pdb_id} entity {identifier} was not found. Please double-check arguments and try again."
             )
         elif resource in need_chain_id:
-            logging.error(
+            logger.error(
                 f"{resource} for {pdb_id} chain {identifier} was not found. Please double-check arguments and try again."
             )
         else:
-            logging.error(
+            logger.error(
                 f"{resource} for {pdb_id} was not found. Please double-check arguments and try again."
             )
         return
diff --git a/gget/gget_ref.py b/gget/gget_ref.py
index 82e3c3bd..da4b65a6 100644
--- a/gget/gget_ref.py
+++ b/gget/gget_ref.py
@@ -1,21 +1,18 @@
 from bs4 import BeautifulSoup
 import requests
 import json
-import logging
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
+# Custom functions
+from .utils import (
+    ref_species_options,
+    find_latest_ens_rel,
+    find_nv_kingdom,
+    set_up_logger,
 )
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
 
-# Custom functions
-from .utils import ref_species_options, find_latest_ens_rel, find_nv_kingdom
+logger = set_up_logger()
 
-from .constants import ENSEMBL_FTP_URL, ENSEMBL_FTP_URL_NV
+from .constants import ENSEMBL_FTP_URL, ENSEMBL_FTP_URL_NV, ENSEMBL_FTP_URL_GRCH37
 
 
 def find_FTP_link(url, link_substring):
@@ -71,6 +68,7 @@ def ref(
     Args:
     - species         Defines the species for which the reference should be fetched in the format "<genus>_<species>",
                       e.g. species = "homo_sapiens".
+                      Supported shortcuts: "human", "mouse", "human_grch37" (accesses the GRCh37 genome assembly)
     - which           Defines which results to return.
                       Default: 'all' -> Returns all available results.
                       Possible entries are one or a combination (as a list of strings) of the following:
@@ -97,12 +95,12 @@ def ref(
     if list_species:
         if release is None:
             if verbose:
-                logging.info(
+                logger.info(
                     f"Fetching available vertebrate genomes (GTF and FASTA available) from Ensembl release {find_latest_ens_rel()} (latest)."
                 )
         else:
             if verbose:
-                logging.info(
+                logger.info(
                     f"Fetching available vertebrate genomes (GTF and FASTA available) from Ensembl release {release}."
                 )
 
@@ -116,8 +114,8 @@ def ref(
         species_list = list(set(species_list_gtf) & set(species_list_dna))
 
         if save:
-            with open("ensembl_species.txt", 'w') as tfile:
-                tfile.write('\n'.join(species_list))
+            with open("ensembl_species.txt", "w") as tfile:
+                tfile.write("\n".join(species_list))
 
         return sorted(species_list)
 
@@ -125,12 +123,12 @@ def ref(
     elif list_iv_species:
         if release is None:
             if verbose:
-                logging.info(
+                logger.info(
                     f"Fetching available invertebrate genomes (GTF and FASTA present) from Ensembl release {find_latest_ens_rel(database=ENSEMBL_FTP_URL_NV)} (latest)."
                 )
         else:
             if verbose:
-                logging.info(
+                logger.info(
                     f"Fetching available invertebrate genomes (GTF and FASTA present) from Ensembl release {release}."
                 )
 
@@ -148,8 +146,8 @@ def ref(
         species_list = list(set(species_list_gtf) & set(species_list_dna))
 
         if save:
-            with open("ensembl_iv_species.txt", 'w') as tfile:
-                tfile.write('\n'.join(species_list))
+            with open("ensembl_iv_species.txt", "w") as tfile:
+                tfile.write("\n".join(species_list))
 
         return sorted(species_list)
 
@@ -171,19 +169,30 @@ def ref(
         )
 
     # Species shortcuts
+    grch37 = False
     if species == "human":
         species = "homo_sapiens"
     if species == "mouse":
         species = "mus_musculus"
+    if species == "human_grch37":
+        species = "homo_sapiens"
+        grch37 = True
 
     # In case species was passed with upper case letters
     species = species.lower()
 
-    ## For non-vertebrates, switch to non-vertebrate databases
-    if species in ref_species_options("dna", database=ENSEMBL_FTP_URL, release=release):
+    # GRCh37 database (releases same as standard database)
+    if grch37:
+        database = ENSEMBL_FTP_URL_GRCH37
+        ENS_rel = find_latest_ens_rel(ENSEMBL_FTP_URL)
+    # Standard database
+    elif species in ref_species_options(
+        "dna", database=ENSEMBL_FTP_URL, release=release
+    ):
         database = ENSEMBL_FTP_URL
         # Find latest vertebrate Ensembl release
-        ENS_rel = find_latest_ens_rel(ENSEMBL_FTP_URL)
+        ENS_rel = find_latest_ens_rel(database)
+    # For non-vertebrates, switch to non-vertebrate databases
     else:
         database = ENSEMBL_FTP_URL_NV
         # Find latest NV Ensembl release
@@ -193,26 +202,31 @@ def ref(
     if release != None:
         # Warn user when release is higher than the latest release
         if release > ENS_rel:
-            logging.warning(
+            logger.warning(
                 f"Provided Ensembl release number {release} is greater than the latest release ({ENS_rel})."
             )
         ENS_rel = release
 
-    ## Raise error if species not found (both FASTA and GTF have to be available)
-    # Find all available species for genome FASTAs for this Ensembl release
-    species_list_dna = ref_species_options("dna", database=database, release=ENS_rel)
-    # Find all available species for GTFs for this Ensembl release
-    species_list_gtf = ref_species_options("gtf", database=database, release=ENS_rel)
-    # Find intersection of the two lists
-    # (Only species which have GTF and FASTAs available can continue)
-    species_list = list(set(species_list_gtf) & set(species_list_dna))
-
-    if species not in species_list:
-        raise ValueError(
-            f"Species does not match any available species for Ensembl release {ENS_rel}. Please double-check spelling.\n"
-            "'gget ref --list_species' -> lists out all available species (Python: 'gget.ref(None, list_species=True)').\n"
-            "Combine with `release` argument to define specific Ensembl release (default: latest).\n"
+    if not grch37:
+        ## Raise error if species not found (both FASTA and GTF have to be available)
+        # Find all available species for genome FASTAs for this Ensembl release
+        species_list_dna = ref_species_options(
+            "dna", database=database, release=ENS_rel
         )
+        # Find all available species for GTFs for this Ensembl release
+        species_list_gtf = ref_species_options(
+            "gtf", database=database, release=ENS_rel
+        )
+        # Find intersection of the two lists
+        # (Only species which have GTF and FASTAs available can continue)
+        species_list = list(set(species_list_gtf) & set(species_list_dna))
+
+        if species not in species_list:
+            raise ValueError(
+                f"Species does not match any available species for Ensembl release {ENS_rel}. Please double-check spelling.\n"
+                "'gget ref --list_species' -> lists out all available species (Python: 'gget.ref(None, list_species=True)').\n"
+                "Combine with `release` argument to define specific Ensembl release (default: latest).\n"
+            )
 
     ## Find kingdom for non-vertebrate species
     if database == ENSEMBL_FTP_URL_NV:
@@ -226,9 +240,14 @@ def ref(
         else:
             gtf_search_url = database + f"release-{ENS_rel}/gtf/{species}/"
 
+        if grch37:
+            link_substring = "GRCh37.87.gtf.gz"
+        else:
+            link_substring = f"{ENS_rel}.gtf.gz"
+
         # Get link, release date and dataset size
         gtf_str, gtf_date, gtf_size = find_FTP_link(
-            url=gtf_search_url, link_substring=f"{ENS_rel}.gtf.gz"
+            url=gtf_search_url, link_substring=link_substring
         )
         # Build the final download link
         if not isinstance(gtf_str, type(None)):
@@ -493,7 +512,7 @@ def ref(
             with open("gget_ref_results.json", "w", encoding="utf-8") as file:
                 json.dump(ref_dict, file, ensure_ascii=False, indent=4)
         if verbose:
-            logging.info(
+            logger.info(
                 f"Fetching reference information for {species} from Ensembl release: {ENS_rel}."
             )
         return ref_dict
@@ -501,7 +520,7 @@ def ref(
     # If FTP==True, return only the specified URLs as a list
     if ftp:
         if verbose:
-            logging.info(
+            logger.info(
                 f"Fetching reference information for {species} from Ensembl release: {ENS_rel}."
             )
         results = []
@@ -531,7 +550,7 @@ def ref(
                 )
 
         if save:
-            with open('gget_ref_results.txt', 'w') as tfile:
-                tfile.write('\n'.join(results))
+            with open("gget_ref_results.txt", "w") as tfile:
+                tfile.write("\n".join(results))
 
         return results
diff --git a/gget/gget_search.py b/gget/gget_search.py
index 25fb8fc4..dc50f06d 100644
--- a/gget/gget_search.py
+++ b/gget/gget_search.py
@@ -1,31 +1,21 @@
 import numpy as np
+import pandas as pd
 import json as json_package
 import mysql.connector as sql
 import time
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 import warnings
 
 warnings.simplefilter(action="ignore", category=UserWarning)
-import pandas as pd
-
 
 # Custom functions
-from gget.utils import (
+from .utils import (
     search_species_options,
     find_latest_ens_rel,
     wrap_cols_func,
     find_nv_kingdom,
+    set_up_logger,
 )
+logger = set_up_logger()
 
 from gget.constants import ENSEMBL_FTP_URL, ENSEMBL_FTP_URL_NV
 
@@ -101,7 +91,7 @@ def search(
     """
     # Handle deprecated arguments
     if seqtype:
-        logging.error(
+        logger.error(
             "'seqtype' argument deprecated! Please use argument 'id_type' instead."
         )
         return
@@ -142,7 +132,7 @@ def search(
     if "core" in species:
         db = species
         if release:
-            logging.warning(
+            logger.warning(
                 "Specified release overwritten because database name was provided."
             )
     else:
@@ -167,7 +157,7 @@ def search(
         # the standard core database will be used
         if len(db) > 1 and "mus_musculus" in species:
             db = f"mus_musculus_core_{ens_rel}_39"
-            logging.warning(
+            logger.warning(
                 f"Defaulting to mus musculus core database: {db}.\n"
                 "All available vertebrate databases can be found here:\n"
                 f"http://ftp.ensembl.org/pub/release-{ens_rel}/mysql/ \n"
@@ -175,7 +165,7 @@ def search(
 
         # Check for ambigious species matches in species other than mouse
         elif len(db) > 1 and "mus_musculus" not in species:
-            logging.warning(
+            logger.warning(
                 f"Species matches more than one database. Defaulting to first database: {db[0]}.\n"
                 "All available databases can be found here:\n"
                 f"Vertebrates: http://ftp.ensembl.org/pub/release-{ens_rel}/mysql/ \n"
@@ -196,7 +186,7 @@ def search(
             db = db[0]
 
     if verbose:
-        logging.info(f"Fetching results from database: {db}")
+        logger.info(f"Fetching results from database: {db}")
 
     ## Connect to data base
     try:
@@ -333,18 +323,18 @@ def search(
     if limit != None:
         # Print number of genes/transcripts found versus fetched
         if verbose:
-            logging.info(f"Returning {limit} matches of {len(df)} total matches found.")
+            logger.info(f"Returning {limit} matches of {len(df)} total matches found.")
         # Remove all but limit rows
         df = df.head(limit)
 
     else:
         # Print number of genes/transcripts fetched
         if verbose:
-            logging.info(f"Total matches found: {len(df)}.")
+            logger.info(f"Total matches found: {len(df)}.")
 
     # Print query time
     if verbose:
-        logging.info(f"Query time: {round(time.time() - start_time, 2)} seconds.")
+        logger.info(f"Query time: {round(time.time() - start_time, 2)} seconds.")
 
     # Remove database numbers to retain only species name
     clean_db = "_".join(db.split("_")[:3]).replace("_core", "")
diff --git a/gget/gget_seq.py b/gget/gget_seq.py
index 8aa56139..a1145956 100644
--- a/gget/gget_seq.py
+++ b/gget/gget_seq.py
@@ -1,18 +1,8 @@
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 import numpy as np
 
 # Custom functions
-from .utils import rest_query, get_uniprot_seqs
+from .utils import rest_query, get_uniprot_seqs, set_up_logger
+logger = set_up_logger()
 from .gget_info import info
 
 # Constants
@@ -50,7 +40,7 @@ def seq(
     """
     # Handle deprecated arguments
     if seqtype:
-        logging.error(
+        logger.error(
             "'seqtype' argument deprecated! Please use True/False argument 'translate' instead."
         )
         return
@@ -72,7 +62,7 @@ def seq(
 
             if "." in ensembl_ID and temp == 0:
                 if verbose:
-                    logging.info(
+                    logger.info(
                         "We noticed that you may have passed a version number with your Ensembl ID.\n"
                         "Please note that gget seq will return information linked to the latest Ensembl ID version."
                     )
@@ -119,12 +109,12 @@ def seq(
                     results_dict[ensembl_ID].update({"seq": df_temp})
 
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting nucleotide sequence of {ensembl_ID} from Ensembl."
                         )
 
                 except RuntimeError:
-                    logging.error(
+                    logger.error(
                         f"ID {ensembl_ID} not found. Please double-check spelling/arguments and try again."
                     )
 
@@ -135,7 +125,7 @@ def seq(
 
                 # Check if Ensembl ID was found
                 if isinstance(info_df, type(None)):
-                    logging.warning(
+                    logger.warning(
                         f"ID '{ensembl_ID}' not found. Please double-check spelling/arguments and try again."
                     )
                     continue
@@ -145,7 +135,7 @@ def seq(
                 # If the ID is a gene, get the IDs of all its transcripts
                 if ens_ID_type == "Gene":
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting nucleotide sequences of all transcripts of {ensembl_ID} from Ensembl."
                         )
 
@@ -173,7 +163,7 @@ def seq(
                             )
 
                         except RuntimeError:
-                            logging.error(
+                            logger.error(
                                 f"ID {transcipt_id} not found. "
                                 "Please double-check spelling/arguments and try again."
                             )
@@ -197,13 +187,13 @@ def seq(
                         # Add results to main dict
                         results_dict[ensembl_ID].update({"seq": df_temp})
 
-                        logging.info(
+                        logger.info(
                             f"Requesting nucleotide sequence of {ensembl_ID} from Ensembl."
                         )
-                        logging.warning("The isoform option only applies to gene IDs.")
+                        logger.warning("The isoform option only applies to gene IDs.")
 
                     except RuntimeError:
-                        logging.error(
+                        logger.error(
                             f"ID {ensembl_ID} not found. "
                             "Please double-check spelling/arguments and try again."
                         )
@@ -238,7 +228,7 @@ def seq(
 
                 # Check that Ensembl ID was found
                 if isinstance(info_df, type(None)):
-                    logging.warning(
+                    logger.warning(
                         f"ID '{ensembl_ID}' not found. Please double-check spelling/arguments."
                     )
                     continue
@@ -268,7 +258,7 @@ def seq(
                         trans_ids.append(temp_trans_id)
 
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting amino acid sequence of the canonical transcript {temp_trans_id} of gene {ensembl_ID} from UniProt."
                         )
 
@@ -281,12 +271,12 @@ def seq(
                     trans_ids.append(ensembl_ID)
 
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting amino acid sequence of {ensembl_ID} from UniProt."
                         )
 
                 else:
-                    logging.warning(
+                    logger.warning(
                         f"{ensembl_ID} not recognized as either a gene or transcript ID. It will not be included in the UniProt query."
                     )
 
@@ -303,7 +293,7 @@ def seq(
 
                 # Check that Ensembl ID was found
                 if isinstance(info_df, type(None)):
-                    logging.warning(
+                    logger.warning(
                         f"ID '{ensembl_ID}' not found. Please double-check spelling/arguments."
                     )
                     continue
@@ -328,7 +318,7 @@ def seq(
                             trans_ids.append(transcipt_id)
 
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting amino acid sequences of all transcripts of gene {ensembl_ID} from UniProt."
                         )
 
@@ -341,13 +331,13 @@ def seq(
                     trans_ids.append(ensembl_ID)
 
                     if verbose:
-                        logging.info(
+                        logger.info(
                             f"Requesting amino acid sequence of {ensembl_ID} from UniProt."
                         )
-                    logging.warning("The isoform option only applies to gene IDs.")
+                    logger.warning("The isoform option only applies to gene IDs.")
 
                 else:
-                    logging.warning(
+                    logger.warning(
                         f"{ensembl_ID} not recognized as either a gene or transcript ID. It will not be included in the UniProt query."
                     )
 
@@ -356,7 +346,7 @@ def seq(
 
         # Check if any results were found
         if len(df_uniprot) < 1:
-            logging.error("No UniProt amino acid sequences were found for these ID(s).")
+            logger.error("No UniProt amino acid sequences were found for these ID(s).")
 
         else:
             # Build FASTA file from UniProt results
diff --git a/gget/gget_setup.py b/gget/gget_setup.py
index a704ee66..4281241d 100644
--- a/gget/gget_setup.py
+++ b/gget/gget_setup.py
@@ -5,15 +5,9 @@
 import platform
 import uuid
 from platform import python_version
-import logging
 
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
+from .utils import set_up_logger
+logger = set_up_logger()
 
 from .compile import PACKAGE_PATH
 from .constants import (
@@ -67,7 +61,7 @@ def setup(module, verbose=True, out=None):
 
     if module == "gpt":
         if verbose:
-            logging.info("Installing openai version <=0.28.1 (requires pip).")
+            logger.info("Installing openai version <=0.28.1 (requires pip).")
         command = "pip install -q -U 'openai<=0.28.1'"
         with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process:
             stderr = process.stderr.read().decode("utf-8")
@@ -76,7 +70,7 @@ def setup(module, verbose=True, out=None):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error(
+            logger.error(
                 "Installation of openai version <=0.28.1 with pip (https://pypi.org/project/openai) failed."
             )
             return
@@ -86,16 +80,16 @@ def setup(module, verbose=True, out=None):
             import openai
 
             if verbose:
-                logging.info(f"openai installed succesfully.")
+                logger.info(f"openai installed succesfully.")
         except ImportError as e:
-            logging.error(
+            logger.error(
                 f"openai installation with pip (https://pypi.org/project/openai) failed. Import error:\n{e}"
             )
             return
 
     if module == "cellxgene":
         if verbose:
-            logging.info("Installing cellxgene-census package (requires pip).")
+            logger.info("Installing cellxgene-census package (requires pip).")
         command = "pip install -q -U cellxgene-census"
         with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process:
             stderr = process.stderr.read().decode("utf-8")
@@ -104,7 +98,7 @@ def setup(module, verbose=True, out=None):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error(
+            logger.error(
                 "cellxgene-census installation with pip (https://pypi.org/project/cellxgene-census) failed."
             )
             return
@@ -114,19 +108,19 @@ def setup(module, verbose=True, out=None):
             import cellxgene_census
 
             if verbose:
-                logging.info(f"cellxgene_census installed succesfully.")
+                logger.info(f"cellxgene_census installed succesfully.")
         except ImportError as e:
-            logging.error(
+            logger.error(
                 f"cellxgene-census installation with pip (https://pypi.org/project/cellxgene-census) failed. Import error:\n{e}"
             )
             return
 
     if module == "elm":
         if verbose:
-            logging.info(
+            logger.info(
                 "ELM data can be downloaded & distributed for non-commercial use according to the following license: http://elm.eu.org/media/Elm_academic_license.pdf"
             )
-            logging.info(
+            logger.info(
                 "Downloading ELM database files (requires curl to be installed)..."
             )
 
@@ -176,37 +170,37 @@ def setup(module, verbose=True, out=None):
 
         # Exit system if the subprocess returned with an error
         if process.wait() != 0:
-            logging.error("ELM database files download failed.")
+            logger.error("ELM database files download failed.")
             return
 
         # Check if files are present
         if os.path.exists(elm_instances_fasta):
             if verbose:
-                logging.info(f"ELM sequences file present.")
+                logger.info(f"ELM sequences file present.")
         else:
-            logging.error("ELM FASTA file missing.")
+            logger.error("ELM FASTA file missing.")
 
         if os.path.exists(elm_classes_tsv):
             if verbose:
-                logging.info("ELM classes file present.")
+                logger.info("ELM classes file present.")
         else:
-            logging.error("ELM classes file missing.")
+            logger.error("ELM classes file missing.")
 
         if os.path.exists(elm_instances_tsv):
             if verbose:
-                logging.info("ELM instances file present.")
+                logger.info("ELM instances file present.")
         else:
-            logging.error("ELM instances file missing.")
+            logger.error("ELM instances file missing.")
 
         if os.path.exists(elm_intdomains_tsv):
             if verbose:
-                logging.info("ELM interactions domains file present.")
+                logger.info("ELM interactions domains file present.")
         else:
-            logging.error("ELM interactions domains file missing.")
+            logger.error("ELM interactions domains file missing.")
 
     if module == "alphafold":
         if platform.system() == "Windows":
-            logging.error(
+            logger.error(
                 "gget setup alphafold and gget alphafold are not supported on Windows OS."
             )
 
@@ -215,7 +209,7 @@ def setup(module, verbose=True, out=None):
             import simtk.openmm as openmm
 
             # Silence openmm logger
-            logging.getLogger("openmm").setLevel(logging.WARNING)
+            logger.getLogger("openmm").setLevel(logger.WARNING)
 
             # Commenting the following out because openmm v7.7.0 does not support __version__
             # # Check if correct version was installed
@@ -223,7 +217,7 @@ def setup(module, verbose=True, out=None):
             #     raise ImportError()
 
             # if verbose:
-            #   logging.info(f"openmm v{openmm.__version__} already installed.")
+            #   logger.info(f"openmm v{openmm.__version__} already installed.")
 
         except ImportError as e:
             raise ImportError(
@@ -239,7 +233,7 @@ def setup(module, verbose=True, out=None):
 
         ## Install py3Dmol
         if verbose:
-            logging.info("Installing py3Dmol (requires pip).")
+            logger.info("Installing py3Dmol (requires pip).")
         command = "pip install -q py3Dmol"
         with subprocess.Popen(command, shell=True, stderr=subprocess.PIPE) as process:
             stderr = process.stderr.read().decode("utf-8")
@@ -248,7 +242,7 @@ def setup(module, verbose=True, out=None):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error(
+            logger.error(
                 "py3Dmol installation with pip (https://pypi.org/project/py3Dmol) failed."
             )
             return
@@ -258,16 +252,16 @@ def setup(module, verbose=True, out=None):
             import py3Dmol
 
             if verbose:
-                logging.info(f"py3Dmol installed succesfully.")
+                logger.info(f"py3Dmol installed succesfully.")
         except ImportError as e:
-            logging.error(
+            logger.error(
                 f"py3Dmol installation with pip (https://pypi.org/project/py3Dmol/) failed. Import error:\n{e}"
             )
             return
 
         ## Install Alphafold if not already installed
         if verbose:
-            logging.info("Installing AlphaFold from source (requires pip and git).")
+            logger.info("Installing AlphaFold from source (requires pip and git).")
 
         ## Install AlphaFold and change jackhmmer directory where database chunks are saved in
         # Define AlphaFold folder name and location
@@ -283,7 +277,7 @@ def setup(module, verbose=True, out=None):
             command = """
                 git clone --branch main -q --branch {} {} {} \
                 && sed -i '' 's/\/tmp\/ramdisk/{}/g' {}/alphafold/data/tools/jackhmmer.py \
-                && sed -i '' 's/from absl import logging/from absl import logging\\\nlogging.set_verbosity(logging.WARNING)/g' {}/alphafold/data/tools/jackhmmer.py \
+                && sed -i '' 's/from absl import logging/from absl import logging\\\logging.set_verbosity(logging.WARNING)/g' {}/alphafold/data/tools/jackhmmer.py \
                 && pip install -q -r {}/requirements.txt \
                 && pip install -q --no-dependencies {}
                 """.format(
@@ -325,7 +319,7 @@ def setup(module, verbose=True, out=None):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error("AlphaFold installation failed.")
+            logger.error("AlphaFold installation failed.")
             return
 
         # Remove cloned directory
@@ -335,9 +329,9 @@ def setup(module, verbose=True, out=None):
             import alphafold as AlphaFold
 
             if verbose:
-                logging.info(f"AlphaFold installed succesfully.")
+                logger.info(f"AlphaFold installed succesfully.")
         except ImportError as e:
-            logging.error(f"AlphaFold installation failed. Import error:\n{e}")
+            logger.error(f"AlphaFold installation failed. Import error:\n{e}")
             return
 
         ## Append AlphaFold to path
@@ -347,7 +341,7 @@ def setup(module, verbose=True, out=None):
 
         ## Install pdbfixer
         if verbose:
-            logging.info("Installing pdbfixer from source (requires pip and git).")
+            logger.info("Installing pdbfixer from source (requires pip and git).")
 
         pdbfixer_folder = os.path.join(
             PACKAGE_PATH, "tmp_pdbfixer_" + str(uuid.uuid4())
@@ -374,7 +368,7 @@ def setup(module, verbose=True, out=None):
             if stderr:
                 # Log the standard error if it is not empty
                 sys.stderr.write(stderr)
-            logging.error("pdbfixer installation failed.")
+            logger.error("pdbfixer installation failed.")
             return
 
         # Remove cloned directory
@@ -386,9 +380,9 @@ def setup(module, verbose=True, out=None):
         pdb_out, err = process.communicate()
 
         if pdb_out.decode() != "":
-            logging.info(f"pdbfixer installed succesfully.")
+            logger.info(f"pdbfixer installed succesfully.")
         else:
-            logging.error("pdbfixer installation failed.")
+            logger.error("pdbfixer installation failed.")
             return
 
         ## Download model parameters
@@ -399,7 +393,7 @@ def setup(module, verbose=True, out=None):
 
         if len(os.listdir(os.path.join(PARAMS_DIR, "params/"))) < 12:
             if verbose:
-                logging.info(
+                logger.info(
                     "Downloading AlphaFold model parameters (requires 4.1 GB of storage). This might take a few minutes."
                 )
             if platform.system() == "Windows":
@@ -425,11 +419,11 @@ def setup(module, verbose=True, out=None):
                     sys.stderr.write(stderr)
             # Exit system if the subprocess returned with an error
             if process.wait() != 0:
-                logging.error("Model parameter download failed.")
+                logger.error("Model parameter download failed.")
                 return
             else:
                 if verbose:
-                    logging.info("Model parameter download complete.")
+                    logger.info("Model parameter download complete.")
         else:
             if verbose:
-                logging.info("AlphaFold model parameters already downloaded.")
+                logger.info("AlphaFold model parameters already downloaded.")
diff --git a/gget/main.py b/gget/main.py
index 8cf5231a..913a9e98 100644
--- a/gget/main.py
+++ b/gget/main.py
@@ -1,16 +1,5 @@
 import argparse
 import sys
-import logging
-
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
-# Mute numexpr threads info
-logging.getLogger("numexpr").setLevel(logging.WARNING)
-
 from datetime import datetime
 
 # Get current date and time for alphafold default foldername
@@ -19,8 +8,13 @@
 import os
 import json
 
-# Custom functions
+from .utils import set_up_logger
+
+logger = set_up_logger()
+
 from .__init__ import __version__
+
+# Module functions
 from .gget_ref import ref
 from .gget_search import search
 from .gget_info import info
@@ -38,6 +32,20 @@
 from .gget_elm import elm
 from .gget_diamond import diamond
 from .gget_cosmic import cosmic
+from .gget_mutate import mutate
+
+
+# Custom formatter for help messages that preserved the text formatting and adds the default value to the end of the help message
+class CustomHelpFormatter(argparse.RawTextHelpFormatter):
+    def _get_help_string(self, action):
+        help_str = action.help if action.help else ""
+        if (
+            "%(default)" not in help_str
+            and action.default is not argparse.SUPPRESS
+            and action.default is not None
+        ):
+            help_str += " (default: %(default)s)"
+        return help_str
 
 
 def main():
@@ -65,7 +73,12 @@ def main():
     ## gget ref subparser
     ref_desc = "Fetch FTPs for reference genomes and annotations by species."
     parser_ref = parent_subparsers.add_parser(
-        "ref", parents=[parent], description=ref_desc, help=ref_desc, add_help=True
+        "ref",
+        parents=[parent],
+        description=ref_desc,
+        help=ref_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # ref parser arguments
     parser_ref.add_argument(
@@ -74,11 +87,8 @@ def main():
         nargs="?",
         default=None,
         help=(
-            """
-            Species or database to be searched. Species should be passed in the format "genus_species", e.g. "homo_sapiens".
-            To pass a specific database, enter the name of the core database and release number, e.g. 'mus_musculus_dba2j_core_105_1'.
-            All available databases for each Ensembl release can be found here: http://ftp.ensembl.org/pub/
-            """
+            "Species or database to be searched. Species should be passed in the format 'genus_species', e.g. 'homo_sapiens'.\n"
+            "Supported shortcuts: 'human', 'mouse', 'human_grch37' (accesses the GRCh37 genome assembly)"
         ),
     )
     parser_ref.add_argument(
@@ -88,10 +98,8 @@ def main():
         action="store_true",
         required=False,
         help=(
-            """
-            List all available vertebrate species from the Ensembl database. 
-            (Combine with `--release` to get the available species from a specific Ensembl release.)
-            """
+            "List all available vertebrate species from the Ensembl database.\n"
+            "(Combine with `--release` to get the available species from a specific Ensembl release.)"
         ),
     )
     parser_ref.add_argument(
@@ -101,10 +109,8 @@ def main():
         action="store_true",
         required=False,
         help=(
-            """
-            List all available invertebrate species from the Ensembl database. 
-            (Combine with `--release` to get the available species from a specific Ensembl release.)
-            """
+            "List all available invertebrate species from the Ensembl database.\n"
+            "(Combine with `--release` to get the available species from a specific Ensembl release.)"
         ),
     )
     parser_ref.add_argument(
@@ -114,18 +120,16 @@ def main():
         type=str,
         required=False,
         help=(
-            """
-        Defines which results to return. \n
-        Default: 'all' -> Returns all available results. \n
-        Possible entries are one or a combination (as a comma-separated list) of the following: \n
-        'gtf' - Returns the annotation (GTF). \n
-        'cdna' - Returns the trancriptome (cDNA). \n
-        'dna' - Returns the genome (DNA). \n
-        'cds - Returns the coding sequences corresponding to Ensembl genes. (Does not contain UTR or intronic sequence.) \n
-        'cdrna' - Returns transcript sequences corresponding to non-coding RNA genes (ncRNA). \n
-        'pep' - Returns the protein translations of Ensembl genes. \n
-        Example: '-w dna,gtf'
-        """
+            "Defines which results to return.\n"
+            "Default: 'all' -> Returns all available results.\n"
+            "Possible entries are one or a combination (as a comma-separated list) of the following:\n"
+            "'gtf' - Returns the annotation (GTF).\n"
+            "'cdna' - Returns the trancriptome (cDNA).\n"
+            "'dna' - Returns the genome (DNA).\n"
+            "'cds - Returns the coding sequences corresponding to Ensembl genes. (Does not contain UTR or intronic sequence.)\n"
+            "'cdrna' - Returns transcript sequences corresponding to non-coding RNA genes (ncRNA).\n"
+            "'pep' - Returns the protein translations of Ensembl genes.\n"
+            "Example: '-w dna,gtf'"
         ),
     )
     parser_ref.add_argument(
@@ -189,6 +193,7 @@ def main():
         description=search_desc,
         help=search_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # Search parser arguments
     parser_gget.add_argument(
@@ -203,14 +208,12 @@ def main():
         type=str,
         required=True,
         help=(
-            """
-            Species or database to be queried, e.g. 'homo_sapiens' or 'arabidopsis_thaliana'.  
-            To pass a specific database, pass the name of the CORE database, e.g. 'mus_musculus_dba2j_core_105_1'.  
-            All available core databases can be found here:  
-            Vertebrates: http://ftp.ensembl.org/pub/current/mysql/  
-            Invertebrates: http://ftp.ensemblgenomes.org/pub/current/ + kingdom + mysql/  
-            Supported shortcuts: 'human', 'mouse'. 
-            """
+            "Species or database to be queried, e.g. 'homo_sapiens' or 'arabidopsis_thaliana'.\n"
+            "To pass a specific database, pass the name of the CORE database, e.g. 'mus_musculus_dba2j_core_105_1'.\n"
+            "All available core databases can be found here:\n"
+            "Vertebrates: http://ftp.ensembl.org/pub/current/mysql/\n"
+            "Invertebrates: http://ftp.ensemblgenomes.org/pub/current/ + kingdom + mysql/\n"
+            "Supported shortcuts: 'human', 'mouse'."
         ),
     )
     parser_gget.add_argument(
@@ -220,12 +223,10 @@ def main():
         type=int,
         required=False,
         help=(
-            """
-            Defines the Ensembl release number from which the files are fetched, e.g. 104.
-            Note: Does not apply to invertebrate species (you can pass a specific core database (which include a release number) to the species argument instead). 
-            This argument is overwritten if a specific database (which includes a release number) is passed to the species argument.
-            Default: None -> latest Ensembl release is used.
-            """
+            "Defines the Ensembl release number from which the files are fetched, e.g. 104.\n"
+            "Note: Does not apply to invertebrate species (you can pass a specific core database (which include a release number) to the species argument instead).\n"
+            "This argument is overwritten if a specific database (which includes a release number) is passed to the species argument.\n"
+            "Default: None -> latest Ensembl release is used."
         ),
     )
     parser_gget.add_argument(
@@ -237,7 +238,7 @@ def main():
         required=False,
         help=(
             "'gene': Returns genes that match the searchwords. (default).\n"
-            "'transcript': Returns transcripts that match the searchwords. \n"
+            "'transcript': Returns transcripts that match the searchwords."
         ),
     )
     parser_gget.add_argument(
@@ -249,7 +250,7 @@ def main():
         required=False,
         help=(
             "'or': Gene descriptions must include at least one of the searchwords (default).\n"
-            "'and': Only return genes whose descriptions include all searchwords.\n"
+            "'and': Only return genes whose descriptions include all searchwords."
         ),
     )
     parser_gget.add_argument(
@@ -313,7 +314,12 @@ def main():
     ## gget elm subparser
     elm_desc = "Locally predicts Eukaryotic Linear Motifs from an amino acid sequence or UniProt Acc using data from the ELM database (http://elm.eu.org/media/Elm_academic_license.pdf)."
     parser_elm = parent_subparsers.add_parser(
-        "elm", parents=[parent], description=elm_desc, help=elm_desc, add_help=True
+        "elm",
+        parents=[parent],
+        description=elm_desc,
+        help=elm_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # elm parser arguments
     parser_elm.add_argument(
@@ -403,6 +409,7 @@ def main():
         description=diamond_desc,
         help=diamond_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     parser_diamond.add_argument(
         "query",
@@ -425,10 +432,8 @@ def main():
         default=None,
         required=False,
         help=(
-            """
-            Path to save DIAMOND database created from reference. 
-            Default: None -> Temporary db file will be deleted after alignment or saved in 'out' if 'out' is provided.
-            """
+            "Path to save DIAMOND database created from reference.\n"
+            "Default: None -> Temporary db file will be deleted after alignment or saved in 'out' if 'out' is provided."
         ),
     )
     parser_diamond.add_argument(
@@ -446,10 +451,8 @@ def main():
         type=str,
         required=False,
         help=(
-            """
-            One of the following:'fast', 'mid-sensitive', 'sensitive', 'more-sensitive', 'very-sensitive' or 'ultra-sensitive'. 
-            Sensitivity of DIAMOND alignment. Default: 'very-sensitive'. 
-            """
+            "One of the following:'fast', 'mid-sensitive', 'sensitive', 'more-sensitive', 'very-sensitive' or 'ultra-sensitive'.\n"
+            "Sensitivity of DIAMOND alignment. Default: 'very-sensitive'."
         ),
     )
     parser_diamond.add_argument(
@@ -467,10 +470,8 @@ def main():
         default=None,
         required=False,
         help=(
-            """
-            Path to DIAMOND binary,  e.g. path/bins/Linux/diamond.
-            Default: None -> Uses DIAMOND binary installed with gget.
-            """
+            "Path to DIAMOND binary,  e.g. path/bins/Linux/diamond.\n"
+            "Default: None -> Uses DIAMOND binary installed with gget."
         ),
     )
     parser_diamond.add_argument(
@@ -495,17 +496,20 @@ def main():
         type=str,
         required=False,
         help=(
-            """
-            Path to folder to save DIAMOND results in, e.g. path/to/directory/results.json. 
-            Default: Standard out, temporary files are deleted.
-            """
+            "Path to folder to save DIAMOND results in, e.g. path/to/directory/results.json.\n"
+            "Default: Standard out, temporary files are deleted."
         ),
     )
 
     ## gget info subparser
     info_desc = "Fetch gene and transcript metadata using Ensembl IDs."
     parser_info = parent_subparsers.add_parser(
-        "info", parents=[parent], description=info_desc, help=info_desc, add_help=True
+        "info",
+        parents=[parent],
+        description=info_desc,
+        help=info_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # info parser arguments
     parser_info.add_argument(
@@ -601,7 +605,12 @@ def main():
     ## gget seq subparser
     seq_desc = "Fetch nucleotide or amino acid sequence (FASTA) of a gene (and all isoforms) or transcript by Ensembl, WormBase or FlyBase ID. "
     parser_seq = parent_subparsers.add_parser(
-        "seq", parents=[parent], description=seq_desc, help=seq_desc, add_help=True
+        "seq",
+        parents=[parent],
+        description=seq_desc,
+        help=seq_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # seq parser arguments
     parser_seq.add_argument(
@@ -678,6 +687,7 @@ def main():
         description=muscle_desc,
         help=muscle_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # muscle parser arguments
     parser_muscle.add_argument(
@@ -730,6 +740,7 @@ def main():
         description=blast_desc,
         help=blast_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # blast parser arguments
     parser_blast.add_argument(
@@ -850,7 +861,12 @@ def main():
         "BLAT a nucleotide or amino acid sequence against any BLAT UCSC assembly."
     )
     parser_blat = parent_subparsers.add_parser(
-        "blat", parents=[parent], description=blat_desc, help=blat_desc, add_help=True
+        "blat",
+        parents=[parent],
+        description=blat_desc,
+        help=blat_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # blat parser arguments
     parser_blat.add_argument(
@@ -933,6 +949,7 @@ def main():
         description=enrichr_desc,
         help=enrichr_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # enrichr parser arguments
     parser_enrichr.add_argument(
@@ -1053,6 +1070,7 @@ def main():
         description=archs4_desc,
         help=archs4_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # archs4 parser arguments
     parser_archs4.add_argument(
@@ -1079,14 +1097,9 @@ def main():
         type=str,
         required=False,
         help=(
-            """
-            'correlation' (default) or 'tissue'.
-            - 'correlation' returns a gene correlation table that contains the
-            100 most correlated genes to the gene of interest. The Pearson
-            correlation is calculated over all samples and tissues in ARCHS4.
-            - 'tissue' returns a tissue expression atlas calculated from
-            human or mouse samples (as defined by 'species') in ARCHS4.
-            """
+            "'correlation' (default) or 'tissue'.\n"
+            "- 'correlation' returns a gene correlation table that contains the 100 most correlated genes to the gene of interest. The Pearson correlation is calculated over all samples and tissues in ARCHS4.\n"
+            "- 'tissue' returns a tissue expression atlas calculated from human or mouse samples (as defined by 'species') in ARCHS4."
         ),
     )
     parser_archs4.add_argument(
@@ -1096,10 +1109,8 @@ def main():
         type=int,
         required=False,
         help=(
-            """
-            Number of correlated genes to return (default: 100).
-            (Only for gene correlation.)
-            """
+            "Number of correlated genes to return (default: 100).\n"
+            "(Only for gene correlation.)"
         ),
     )
     parser_archs4.add_argument(
@@ -1165,6 +1176,7 @@ def main():
         description=setup_desc,
         help=setup_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # setup parser arguments
     parser_setup.add_argument(
@@ -1200,6 +1212,7 @@ def main():
         description=alphafold_desc,
         help=alphafold_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # alphafold parser arguments
     parser_alphafold.add_argument(
@@ -1223,10 +1236,8 @@ def main():
         type=int,
         required=False,
         help=(
-            """
-            The multimer model will continue recycling until the predictions stop changing, up to the limit set here.
-            For higher accuracy, at the potential cost of longer inference times, set this to 20.
-            """
+            "The multimer model will continue recycling until the predictions stop changing, up to the limit set here.\n"
+            "For higher accuracy, at the potential cost of longer inference times, set this to 20."
         ),
     )
     parser_alphafold.add_argument(
@@ -1264,6 +1275,7 @@ def main():
         description=pdb_desc,
         help=pdb_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     # alphafold parser arguments
     parser_pdb.add_argument(
@@ -1291,20 +1303,18 @@ def main():
         ],
         required=False,
         help=(
-            """
-            Defines type of information to be returned.
-            "pdb": Returns the protein structure in PDB format.
-            "entry": Information about PDB structures at the top level of PDB structure hierarchical data organization.
-            "pubmed": Get PubMed annotations (data integrated from PubMed) for a given entry's primary citation.
-            "assembly": Information about PDB structures at the quaternary structure level.
-            "branched_entity": Get branched entity description (define entity ID as "identifier").
-            "nonpolymer_entity": Get non-polymer entity data (define entity ID as "identifier").
-            "polymer_entity": Get polymer entity data (define entity ID as "identifier").
-            "uniprot": Get UniProt annotations for a given macromolecular entity (define entity ID as "identifier").
-            "branched_entity_instance": Get branched entity instance description (define chain ID as "identifier").
-            "polymer_entity_instance": Get polymer entity instance (a.k.a chain) data (define chain ID as "identifier").
-            "nonpolymer_entity_instance": Get non-polymer entity instance description (define chain ID as "identifier").
-            """
+            "Defines type of information to be returned.\n\n"
+            '"pdb": Returns the protein structure in PDB format.\n'
+            '"entry": Information about PDB structures at the top level of PDB structure hierarchical data organization.\n'
+            '"pubmed": Get PubMed annotations (data integrated from PubMed) for a given entry\'s primary citation.\n'
+            '"assembly": Information about PDB structures at the quaternary structure level.\n'
+            '"branched_entity": Get branched entity description (define entity ID as "identifier").\n'
+            '"nonpolymer_entity": Get non-polymer entity data (define entity ID as "identifier").\n'
+            '"polymer_entity": Get polymer entity data (define entity ID as "identifier").\n'
+            '"uniprot": Get UniProt annotations for a given macromolecular entity (define entity ID as "identifier").\n'
+            '"branched_entity_instance": Get branched entity instance description (define chain ID as "identifier").\n'
+            '"polymer_entity_instance": Get polymer entity instance (a.k.a chain) data (define chain ID as "identifier").\n'
+            '"nonpolymer_entity_instance": Get non-polymer entity instance description (define chain ID as "identifier").'
         ),
     )
     parser_pdb.add_argument(
@@ -1314,10 +1324,8 @@ def main():
         type=str,
         required=False,
         help=(
-            """
-            Can be used to define assembly, entity or chain ID if applicable (default: None).
-            Assembly/entity IDs are numbers (e.g. 1), and chain IDs are letters (e.g. A).
-            """
+            "Can be used to define assembly, entity or chain ID if applicable (default: None).\n"
+            "Assembly/entity IDs are numbers (e.g. 1), and chain IDs are letters (e.g. A)."
         ),
     )
     parser_pdb.add_argument(
@@ -1340,6 +1348,7 @@ def main():
         description=gpt_desc,
         help=gpt_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     parser_gpt.add_argument(
         "prompt",
@@ -1442,6 +1451,7 @@ def main():
         description=cellxgene_desc,
         help=cellxgene_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     parser_cellxgene.add_argument(
         "-o",
@@ -1503,11 +1513,8 @@ def main():
             "cell_type",
         ],
         help=(
-            """
-            List of metadata columns to return (stored in .obs).
-            Default: ["dataset_id", "assay", "suspension_type", "sex", "tissue_general", "tissue", "cell_type"]
-            For more options see: https://api.cellxgene.cziscience.com/curation/ui/#/ -> Schemas -> dataset
-            """
+            "List of metadata columns to return (stored in .obs).\n"
+            "For more options see: https://api.cellxgene.cziscience.com/curation/ui/#/ -> Schemas -> dataset"
         ),
     )
     parser_cellxgene.add_argument(
@@ -1701,10 +1708,13 @@ def main():
         description=cosmic_desc,
         help=cosmic_desc,
         add_help=True,
+        formatter_class=CustomHelpFormatter,
     )
     parser_cosmic.add_argument(
         "searchterm",
         type=str,
+        nargs="?",
+        default=None,
         help="Search term, which can be a mutation, or gene (or Ensembl ID), or sample, etc. as defined using the 'entity' argument. Example: 'EGFR'",
     )
     parser_cosmic.add_argument(
@@ -1714,7 +1724,7 @@ def main():
             "mutations",
             "genes",
             "cancer",
-            "tumour site",
+            "tumour_site",
             "studies",
             "pubmed",
             "samples",
@@ -1740,14 +1750,66 @@ def main():
         required=False,
         help="Returns results in csv format instead of json.",
     )
+    parser_cosmic.add_argument(
+        "-d",
+        "--download_cosmic",
+        default=False,
+        action="store_true",
+        required=False,
+        help="Switch into database download mode.",
+    )
+    parser_cosmic.add_argument(
+        "-mc",
+        "--mutation_class",
+        choices=[
+            "cancer",
+            "cell_line",
+            "census",
+            "resistance",
+            "screen",
+            "cancer_example",
+        ],
+        default="cancer",
+        type=str,
+        required=False,
+        help="Type of COSMIC database to download (only for use with --download_cosmic).",
+    )
+    parser_cosmic.add_argument(
+        "-cv",
+        "--cosmic_version",
+        default=99,
+        type=int,
+        required=False,
+        help="Version of the COSMIC database (only for use with --download_cosmic).",
+    )
+    parser_cosmic.add_argument(
+        "-gv",
+        "--grch_version",
+        default=37,
+        type=int,
+        choices=[37, 38],
+        required=False,
+        help="Version of the human GRCh reference genome (only for use with --download_cosmic).",
+    )
+    parser_cosmic.add_argument(
+        "-gm",
+        "--gget_mutate",
+        default=True,
+        action="store_false",
+        required=False,
+        help="Do NOT create a modified version of the database for use with gget mutate (only for use with --download_cosmic).",
+    )
     parser_cosmic.add_argument(
         "-o",
         "--out",
         type=str,
         required=False,
+        default=None,
         help=(
-            "Path to the file the results will be saved in, e.g. path/to/directory/results.json.\n"
-            "Default: Standard out."
+            "Path to the file (or folder when downloading databases with the download_cosmic flag) the results will be saved in, e.g. path/to/results.json.\n"
+            "Default: None\n"
+            "-> When download_cosmic=False: Results will be returned to standard out\n"
+            "-> When download_cosmic=True: Database will be downloaded into current working directory"
         ),
     )
     parser_cosmic.add_argument(
@@ -1759,6 +1821,109 @@ def main():
         help="Do not print progress information.",
     )
 
+    # mutate parser arguments
+    mutate_desc = "Mutate nucleotide sequences according to defined mutations."
+    parser_mutate = parent_subparsers.add_parser(
+        "mutate",
+        parents=[parent],
+        description=mutate_desc,
+        help=mutate_desc,
+        add_help=True,
+        formatter_class=CustomHelpFormatter,
+    )
+    parser_mutate.add_argument(
+        "sequences",
+        type=str,
+        nargs="+",
+        help=(
+            "(str) Path to the fasta file containing the sequences to be mutated, e.g., 'seqs.fa'.\n"
+            "Sequence identifiers following the '>' character must correspond to the identifiers\n"
+            "in the seq_ID column of 'mutations'.\n"
+            "NOTE: Only string until first space or dot will be used as sequence identifier\n"
+            "- Version numbers of Ensembl IDs will be ignored.\n\n"
+            "Example:\n"
+            ">seq1 (or ENSG00000106443)\n"
+            "ACTGCGATAGACT\n"
+            ">seq2\n"
+            "AGATCGCTAG\n\n"
+            "Alternatively: Input sequence(s) as a string or list, e.g. 'AGCTAGCT' or 'ACTGCTAGCT' 'AGCTAGCT'."
+        ),
+    )
+    parser_mutate.add_argument(
+        "-m",
+        "--mutations",
+        type=str,
+        nargs="+",
+        required=True,
+        help=(
+            "Path to csv or tsv file (e.g., 'mutations.csv') containing information about the mutations in the following format:\n\n"
+            "| mutation             | mut_ID | seq_ID |\n"
+            "| c.1252C>T            | mut1   | seq1   |\n"
+            "| c.2239_2253inv       | mut2   | seq2   |\n"
+            "| c.2239_2253inv       | mut2   | seq3   |\n"
+            "| c.2239_2253delinsAAT | mut3   | seq3   |\n"
+            "| ...                  | ...    | ...    |\n\n"
+            "'mutation' = Column containing the mutations to be performed written in standard mutation annotation (see below)\n"
+            "'mut_ID' = Column containing an identifier for each mutation\n"
+            "'seq_ID' = Column containing the identifiers of the sequences to be mutated\n"
+            "(sequence IDs must correspond to the string following the > character in the input fasta; do NOT include spaces or dots)\n\n"
+            "Alternatively: Input mutation(s) as a string or list, e.g. 'c.2C>T' or 'c.2C>T' 'c.1A>C'\n\n"
+            "NOTE: Enclose individual mutation annotations in quotation marks to prevent terminal parsing errors.\n"
+            "If a list is passed, the number of mutations must equal the number of input sequences."
+        ),
+    )
+    parser_mutate.add_argument(
+        "-k",
+        "--k",
+        default=31,
+        type=int,
+        required=False,
+        help="Length of sequences flanking the mutation. If k > total length of the sequence, the entire sequence will be kept.",
+    )
+    parser_mutate.add_argument(
+        "-mc",
+        "--mut_column",
+        default="mutation",
+        type=str,
+        required=False,
+        help="Name of the column containing the mutations to be performed in 'mutations'.",
+    )
+    parser_mutate.add_argument(
+        "-mic",
+        "--mut_id_column",
+        default="mut_ID",
+        type=str,
+        required=False,
+        help="Name of the column containing the IDs of each mutation in 'mutations'.",
+    )
+    parser_mutate.add_argument(
+        "-sic",
+        "--seq_id_column",
+        default="seq_ID",
+        type=str,
+        required=False,
+        help="Name of the column containing the IDs of the sequences to be mutated in 'mutations'.",
+    )
+    parser_mutate.add_argument(
+        "-o",
+        "--out",
+        default="gget_mutate_out.fa",
+        type=str,
+        required=False,
+        help=(
+            "Path to output fasta file containing the mutated sequences. The identifiers (following the '>') of the mutated sequences in the output fasta will be '>[seq_ID]_[mut_ID]'."
+            "Set to None to return a list of the mutated sequences to standard out."
+        ),
+    )
+    parser_mutate.add_argument(
+        "-q",
+        "--quiet",
+        default=True,
+        action="store_false",
+        required=False,
+        help="Do not print progress information.",
+    )
+
     ### Define return values
     args = parent_parser.parse_args()
 
@@ -1806,6 +1971,7 @@ def main():
         "elm": parser_elm,
         "diamond": parser_diamond,
         "cosmic": parser_cosmic,
+        "mutate": parser_mutate,
     }
 
     if len(sys.argv) == 2:
@@ -1872,12 +2038,12 @@ def main():
     if args.command == "blat":
         # Handle deprecated flags for backwards compatibility
         if args.seq_deprecated and args.sequence:
-            logging.warning(
+            logger.warning(
                 "The [-seq][--sequence] argument is deprecated, using positional argument [sequence] instead."
             )
         if args.seq_deprecated and not args.sequence:
             args.sequence = args.seq_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-seq][--sequence] argument is deprecated, please use positional argument [sequence] instead."
             )
         if not args.seq_deprecated and not args.sequence:
@@ -1922,12 +2088,12 @@ def main():
     if args.command == "blast":
         # Handle deprecated flags for backwards compatibility
         if args.seq_deprecated and args.sequence:
-            logging.warning(
+            logger.warning(
                 "The [-seq][--sequence] argument is deprecated, using positional argument [sequence] instead."
             )
         if args.seq_deprecated and not args.sequence:
             args.sequence = args.seq_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-seq][--sequence] argument is deprecated, please use positional argument [sequence] instead."
             )
         if not args.seq_deprecated and not args.sequence:
@@ -1972,6 +2138,35 @@ def main():
             if not args.out and args.csv:
                 print(json.dumps(blast_results, ensure_ascii=False, indent=4))
 
+    ## mutate return
+    if args.command == "mutate":
+        if isinstance(args.sequences, list) and len(args.sequences) == 1:
+            seqs = args.sequences[0]
+        else:
+            seqs = args.sequences
+
+        if isinstance(args.mutations, list) and len(args.mutations) == 1:
+            muts = args.mutations[0]
+        else:
+            muts = args.mutations
+
+        # Run mutate function (automatically saves output)
+        mutate_results = mutate(
+            sequences=seqs,
+            mutations=muts,
+            k=args.k,
+            mut_column=args.mut_column,
+            mut_id_column=args.mut_id_column,
+            seq_id_column=args.seq_id_column,
+            out=args.out,
+            verbose=args.quiet,
+        )
+
+        # Print list of mutated sequences if any are returned (this should only happen when out=None)
+        if mutate_results:
+            for mut_seq in mutate_results:
+                print(mut_seq)
+
     ## cosmic return
     if args.command == "cosmic":
         # Run gget cosmic function
@@ -1979,46 +2174,52 @@ def main():
             searchterm=args.searchterm,
             entity=args.entity,
             limit=args.limit,
-            verbose=args.quiet,
             json=args.csv,
+            download_cosmic=args.download_cosmic,
+            mutation_class=args.mutation_class,
+            cosmic_version=args.cosmic_version,
+            grch_version=args.grch_version,
+            gget_mutate=args.gget_mutate,
+            out=args.out,
+            verbose=args.quiet,
         )
 
-        # Check if the function returned something
-        if not isinstance(cosmic_results, type(None)):
-            # Save blast results if args.out specified
-            if args.out and not args.csv:
-                # Create saving directory
-                directory = "/".join(args.out.split("/")[:-1])
-                if directory != "":
-                    os.makedirs(directory, exist_ok=True)
-                # Save to csv
-                cosmic_results.to_csv(args.out, index=False)
-
-            if args.out and args.csv:
-                # Create saving directory
-                directory = "/".join(args.out.split("/")[:-1])
-                if directory != "":
-                    os.makedirs(directory, exist_ok=True)
-                # Save json
-                with open(args.out, "w", encoding="utf-8") as f:
-                    json.dump(cosmic_results, f, ensure_ascii=False, indent=4)
-
-            # Print results if no directory specified
-            if not args.out and not args.csv:
-                cosmic_results.to_csv(sys.stdout, index=False)
-            if not args.out and args.csv:
-                print(json.dumps(cosmic_results, ensure_ascii=False, indent=4))
+        # # Check if the function returned something
+        # if not isinstance(cosmic_results, type(None)):
+        #     # Save cosmic results if args.out specified
+        #     if args.out and not args.csv:
+        #         # Create saving directory
+        #         directory = "/".join(args.out.split("/")[:-1])
+        #         if directory != "":
+        #             os.makedirs(directory, exist_ok=True)
+        #         # Save to csv
+        #         cosmic_results.to_csv(args.out, index=False)
+
+        #     if args.out and args.csv:
+        #         # Create saving directory
+        #         directory = "/".join(args.out.split("/")[:-1])
+        #         if directory != "":
+        #             os.makedirs(directory, exist_ok=True)
+        #         # Save json
+        #         with open(args.out, "w", encoding="utf-8") as f:
+        #             json.dump(cosmic_results, f, ensure_ascii=False, indent=4)
+
+        #     # Print results if no directory specified
+        #     if not args.out and not args.csv:
+        #         cosmic_results.to_csv(sys.stdout, index=False)
+        #     if not args.out and args.csv:
+        #         print(json.dumps(cosmic_results, ensure_ascii=False, indent=4))
 
     ## archs4 return
     if args.command == "archs4":
         # Handle deprecated flags for backwards compatibility
         if args.gene_deprecated and args.gene:
-            logging.warning(
+            logger.warning(
                 "The [-g][--gene] argument is deprecated, using positional argument [gene] instead."
             )
         if args.gene_deprecated and not args.gene:
             args.gene = args.gene_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-g][--gene] argument is deprecated, please use positional argument [gene] instead."
             )
         if not args.gene_deprecated and not args.gene:
@@ -2065,12 +2266,12 @@ def main():
     if args.command == "muscle":
         # Handle deprecated flags for backwards compatibility
         if args.fasta_deprecated and args.fasta:
-            logging.warning(
+            logger.warning(
                 "The [-fa][--fasta] argument is deprecated, using positional argument [fasta] instead."
             )
         if args.fasta_deprecated and not args.fasta:
             args.fasta = args.fasta_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-fa][--fasta] argument is deprecated, please use positional argument [fasta] instead."
             )
         if not args.fasta_deprecated and not args.fasta:
@@ -2131,8 +2332,8 @@ def main():
                 directory = "/".join(args.out.split("/")[:-1])
                 if directory != "":
                     os.makedirs(directory, exist_ok=True)
-                with open(args.out, 'w') as tfile:
-                    tfile.write('\n'.join(species_list))
+                with open(args.out, "w") as tfile:
+                    tfile.write("\n".join(species_list))
             else:
                 for species in species_list:
                     print(species)
@@ -2147,20 +2348,20 @@ def main():
                 directory = "/".join(args.out.split("/")[:-1])
                 if directory != "":
                     os.makedirs(directory, exist_ok=True)
-                with open(args.out, 'w') as tfile:
-                    tfile.write('\n'.join(species_list))
+                with open(args.out, "w") as tfile:
+                    tfile.write("\n".join(species_list))
             else:
                 for species in species_list:
                     print(species)
 
         # Handle deprecated flags for backwards compatibility
         if args.species_deprecated and args.species:
-            logging.warning(
+            logger.warning(
                 "The [-s][--species] argument is deprecated, using positional argument [species] instead."
             )
         if args.species_deprecated and not args.species:
             args.species = args.species_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-s][--species] argument is deprecated, please use positional argument [species] instead."
             )
 
@@ -2198,8 +2399,8 @@ def main():
                     directory = "/".join(args.out.split("/")[:-1])
                     if directory != "":
                         os.makedirs(directory, exist_ok=True)
-                    with open(args.out, 'w') as tfile:
-                        tfile.write('\n'.join(ref_results))
+                    with open(args.out, "w") as tfile:
+                        tfile.write("\n".join(ref_results))
 
                     if args.download == True:
                         # Download list of URLs
@@ -2208,7 +2409,7 @@ def main():
                             command = "curl -O " + link
                             os.system(command)
                 #                     else:
-                #                         logging.info(
+                #                         logger.info(
                 #                             "To download the FTPs to the current directory, add flag [-d]."
                 #                         )
 
@@ -2225,7 +2426,7 @@ def main():
                             command = "curl -O " + link
                             os.system(command)
             #                     else:
-            #                         logging.info(
+            #                         logger.info(
             #                             "To download the FTPs to the current directory, add flag [-d]."
             #                         )
 
@@ -2249,7 +2450,7 @@ def main():
                                     command = "curl -O " + link
                                     os.system(command)
                 #                     else:
-                #                         logging.info(
+                #                         logger.info(
                 #                             "To download the FTPs to the current directory, add flag [-d]."
                 #                         )
 
@@ -2267,7 +2468,7 @@ def main():
                                     command = "curl -O " + link
                                     os.system(command)
     #                     else:
-    #                         logging.info(
+    #                         logger.info(
     #                             "To download the FTPs to the current directory, add flag [-d]."
     #                         )
 
@@ -2275,12 +2476,12 @@ def main():
     if args.command == "search":
         # Handle deprecated flags for backwards compatibility
         if args.sw_deprecated and args.searchwords:
-            logging.warning(
+            logger.warning(
                 "The [-sw][--searchwords] argument is deprecated, using positional argument [searchwords] instead."
             )
         if args.sw_deprecated and not args.searchwords:
             args.searchwords = args.sw_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-sw][--searchwords] argument is deprecated, please use positional argument [searchwords] instead."
             )
         if not args.sw_deprecated and not args.searchwords:
@@ -2338,12 +2539,12 @@ def main():
     if args.command == "enrichr":
         # Handle deprecated flags for backwards compatibility
         if args.genes_deprecated and args.genes:
-            logging.warning(
+            logger.warning(
                 "The [-g][--genes] argument is deprecated, using positional argument [genes] instead."
             )
         if args.genes_deprecated and not args.genes:
             args.genes = args.genes_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-g][--genes] argument is deprecated, please use positional argument [genes] instead."
             )
         if not args.genes_deprecated and not args.genes:
@@ -2419,16 +2620,16 @@ def main():
     if args.command == "info":
         # Handle deprecated flags for backwards compatibility
         if args.id_deprecated and args.ens_ids:
-            logging.warning(
+            logger.warning(
                 "The [-id][--ens_ids] argument is deprecated, using positional argument [ens_ids] instead."
             )
         if args.id_deprecated and not args.ens_ids:
             args.ens_ids = args.id_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-id][--genes] argument is deprecated, please use arguments [ens_ids] instead."
             )
         if args.ensembl_only:
-            logging.warning(
+            logger.warning(
                 "The [-eo][--ensembl_only] argument is deprecated, please use arguments [ncbi] and [uniprot] instead."
             )
         if not args.id_deprecated and not args.ens_ids:
@@ -2487,12 +2688,12 @@ def main():
     if args.command == "seq":
         # Handle deprecated flags for backwards compatibility
         if args.id_deprecated and args.ens_ids:
-            logging.warning(
+            logger.warning(
                 "The [-id][--ens_ids] argument is deprecated, using positional argument [ens_ids] instead."
             )
         if args.id_deprecated and not args.ens_ids:
             args.ens_ids = args.id_deprecated
-            logging.warning(
+            logger.warning(
                 "The [-id][--ens_ids] argument is deprecated, please use positional argument [ens_ids] instead."
             )
         if not args.id_deprecated and not args.ens_ids:
diff --git a/gget/utils.py b/gget/utils.py
index 4e11bce1..7cdf78ea 100644
--- a/gget/utils.py
+++ b/gget/utils.py
@@ -10,24 +10,91 @@
 import numpy as np
 from IPython.display import display, HTML
 import logging
+from datetime import datetime
 
-# Add and format time stamp in logging messages
-logging.basicConfig(
-    format="%(asctime)s %(levelname)s %(message)s",
-    level=logging.INFO,
-    datefmt="%c",
-)
 # Mute numexpr threads info
 logging.getLogger("numexpr").setLevel(logging.WARNING)
 
 from .constants import ENSEMBL_FTP_URL, ENSEMBL_FTP_URL_NV, ENS_TO_PDB_API
 
+def set_up_logger():
+    logging_level = logging.INFO
+    logger = logging.getLogger(__name__)
+    logger.setLevel(logging_level)
+
+    if not logger.hasHandlers():
+        formatter = logging.Formatter("%(asctime)s - %(levelname)s - %(message)s", "%H:%M:%S")
+
+        console_handler = logging.StreamHandler()
+        console_handler.setFormatter(formatter)
+        logger.addHandler(console_handler)
+
+        # log_dir = "logs"
+        # if not os.path.exists(log_dir):
+        #     os.makedirs(log_dir)
+
+        # log_file = os.path.join(
+        #     log_dir, f"log_file_{datetime.now().strftime('%Y%m%d_%H%M%S')}.txt"
+        # )
+
+        # file_handler = logging.FileHandler(log_file)
+        # file_handler.setFormatter(formatter)
+        # logger.addHandler(file_handler)
+
+    return logger
+
+logger = set_up_logger()
+
 def flatten(xss):
     """
     Function to flatten a list of lists.
     """
     return [x for xs in xss for x in xs]
 
+
+def read_fasta(fasta):
+    """
+    Args:
+    - fasta     (str) Path to fasta file.
+
+    Returns titles and seqs from fasta file as two list objects.
+    """
+    titles = []
+    seqs = []
+    title_last = False
+    new_seq = False
+    with open(fasta) as fasta_file:
+        for i, line in enumerate(fasta_file):
+            if i == 0 and line[0] != ">":
+                raise ValueError("Expected FASTA file to start with a '>' character. ")
+            elif line[0] == ">":
+                if title_last:
+                    raise ValueError(
+                        "FASTA file contains two lines starting with '>' in a row -> missing sequence line. "
+                    )
+
+                if new_seq:
+                    # Append last recorded sequence
+                    seqs.append(new_seq)
+
+                # Append title line to titles list
+                titles.append(line.strip().replace(">", ""))
+                title_last = True
+
+            else:
+                if title_last:
+                    # Start recording new sequence if title was last
+                    new_seq = line.strip()
+                else:
+                    new_seq = new_seq + line.strip()
+                title_last = False
+
+        # Append last sequence
+        seqs.append(new_seq)
+
+    return titles, seqs
+
+
 def n_colors(nucleotide):
     """
     Returns a string format to print the nucleotide
@@ -165,7 +232,7 @@ def get_uniprot_seqs(server, ensembl_ids):
         # Submit server request
         r = requests.get(server + id_ + "+AND+reviewed:true")
         if not r.ok:
-            logging.error(
+            logger.error(
                 f"UniProt server request returned with error status code: {r.status_code}. Please double-check arguments or try again later."
             )
         # Convert to json
@@ -176,7 +243,7 @@ def get_uniprot_seqs(server, ensembl_ids):
             # Submit server request
             r = requests.get(server + id_)
             if not r.ok:
-                logging.error(
+                logger.error(
                     f"UniProt server request returned with error status code: {r.status_code}. Please double-check arguments or try again later."
                 )
             # Convert to json
@@ -184,7 +251,7 @@ def get_uniprot_seqs(server, ensembl_ids):
 
             # Warn user if unreviewed results were found
             if len(json["results"]) > 0:
-                logging.warning(
+                logger.warning(
                     f"No reviewed UniProt results were found for ID {id_}. Returning all unreviewed results."
                 )
 
@@ -227,7 +294,7 @@ def get_uniprot_seqs(server, ensembl_ids):
 
         else:
             # If no results were found, warn user and do nothing -> returns empty df
-            logging.warning(f"No UniProt sequences were found for ID {id_}.")
+            logger.warning(f"No UniProt sequences were found for ID {id_}.")
 
     return master_df
 
@@ -247,7 +314,7 @@ def get_uniprot_info(server, ensembl_id, verbose=True):
     # Submit server request for reviewed entries
     r = requests.get(server + ensembl_id + "+AND+reviewed:true")
     if not r.ok:
-        logging.error(
+        logger.error(
             f"UniProt server request returned with error status code: {r.status_code}. Please double-check arguments or try again later."
         )
     # Convert to json
@@ -258,7 +325,7 @@ def get_uniprot_info(server, ensembl_id, verbose=True):
         # Submit server request
         r = requests.get(server + ensembl_id)
         if not r.ok:
-            logging.error(
+            logger.error(
                 f"UniProt server request returned with error status code: {r.status_code}. Please double-check arguments or try again later."
             )
         # Convert to json
@@ -267,7 +334,7 @@ def get_uniprot_info(server, ensembl_id, verbose=True):
         # Warn user if unreviewed results were found
         if len(json["results"]) > 0:
             if verbose is True:
-                logging.warning(
+                logger.warning(
                     f"No reviewed UniProt results were found for ID {ensembl_id}. Returning all unreviewed results."
                 )
 
@@ -423,7 +490,7 @@ def get_uniprot_info(server, ensembl_id, verbose=True):
 #         try:
 #             response.raise_for_status()
 #         except requests.HTTPError:
-#             logging.error(
+#             logger.error(
 #                 f"UniProt ID mapping to fetch PDB IDs returned HTTP Error:\n{response.json()}"
 #             )
 #             return
@@ -454,7 +521,7 @@ def get_uniprot_info(server, ensembl_id, verbose=True):
 #             if "jobStatus" in j:
 #                 # Sleep for POLLING_INTERVAL seconds if job status if "RUNNING"
 #                 if j["jobStatus"] == "RUNNING":
-#                     logging.info("Checking if PDB IDs are available...")
+#                     logger.info("Checking if PDB IDs are available...")
 #                     time.sleep(POLLING_INTERVAL)
 #                 else:
 #                     # Raise error if job status is other than "RUNNING"
@@ -611,7 +678,7 @@ def search_species_options(database=ENSEMBL_FTP_URL, release=None):
     # If release != None, use user-defined Ensembl release
     if release != None:
         if release > ENS_rel:
-            logging.warning(
+            logger.warning(
                 f"Provided Ensembl release number {release} is greater than the latest release ({ENS_rel})."
             )
         ENS_rel = release
@@ -658,6 +725,7 @@ def search_species_options(database=ENSEMBL_FTP_URL, release=None):
 
     return databases
 
+
 def find_nv_kingdom(species, release):
     kds = ["plants", "protists", "metazoa", "fungi"]
     for kingdom in kds:
@@ -681,6 +749,7 @@ def find_nv_kingdom(species, release):
         if species in sps[5:]:
             return kingdom
 
+
 def ref_species_options(which, database=ENSEMBL_FTP_URL, release=None):
     """
     Function to find all available species for gget ref.
@@ -700,7 +769,7 @@ def ref_species_options(which, database=ENSEMBL_FTP_URL, release=None):
     if release != None:
         # Warn user if user-defined release is higher than the latest release
         if release > ENS_rel:
-            logging.warning(
+            logger.warning(
                 f"Provided Ensembl release number {release} is greater than the latest release ({ENS_rel})."
             )
         ENS_rel = release


Test patch:
diff --git a/tests/fixtures.py b/tests/fixtures.py
index deb04e94..de7ab515 100644
--- a/tests/fixtures.py
+++ b/tests/fixtures.py
@@ -1,5 +1,5 @@
 # Latest Ensembl release for unittests
-LATEST_ENS_RELEASE = 111
+LATEST_ENS_RELEASE = 112
 
 # gget search species options for Ensembl release 106
 SPECIES_OPTIONS = [
diff --git a/tests/fixtures/muscle_nt_test.fa b/tests/fixtures/muscle_nt_test.fa
index 2cb40372..e4f11ae3 100644
--- a/tests/fixtures/muscle_nt_test.fa
+++ b/tests/fixtures/muscle_nt_test.fa
@@ -1,4 +1,6 @@
 >ENSTGUT00000006367 primary_assembly:bTaeGut1_v1.p:1:107513786:107528106:-1
-TAGTTGTGATACTGGTCAGTGTCACTGGGGCTGTTTCACCCTGTCATTTGAAGACTTCCATTT
+TAGTTGTGATACTGGTCAGTGTCACTGGGGCT
+GTTTCACCCTGTCATTTGAAGACTTCCATTT
 >ENSTGUT00000027003 primary_assembly:bTaeGut1_v1.p:1:107513786:107526965:-1
 CTTTGTTGGAGGCTTTTTGCTGGGGCTGGCGTCGTAAGGACCCGAACATCTCTTCCAGATGGGAT
+
diff --git a/tests/test_mutate.py b/tests/test_mutate.py
new file mode 100644
index 00000000..6eb153bf
--- /dev/null
+++ b/tests/test_mutate.py
@@ -0,0 +1,115 @@
+import unittest
+
+from gget.gget_mutate import (
+    create_mutant_sequence,
+    substitution_mutation,
+    deletion_mutation,
+    insertion_mutation,
+    delins_mutation,
+    duplication_mutation,
+    inversion_mutation,
+)
+
+
+def alphabet_sequence():
+    return "ABCDEFG"
+
+
+def alphabet_sequence_long():
+    return (
+        "ABCDEFGHIJKLMNOPQRSTUVWXYZABCDEFGHIJKLMNOPQRSTUVWXYZABCDEFGHIJKLMNOPQRSTUVWXYZ"
+    )
+
+
+def create_test_row(mutation, sequence):
+    return {"mutation": mutation, "full_sequence": sequence}
+
+
+class TestMutate(unittest.TestCase):
+    def test_substitution(alphabet_sequence):
+        test_row = create_test_row("c.3C>X", alphabet_sequence)
+        result = create_mutant_sequence(test_row, substitution_mutation)
+        assert result == "ABXDEFG", f"Expected ABXDEFG, got {result}"
+
+    def test_substitution_long(alphabet_sequence_long):
+        test_row = create_test_row("c.35I>X", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, substitution_mutation)
+        assert (
+            result == "EFGHIJKLMNOPQRSTUVWXYZABCDEFGHXJKLMNOPQRSTUVWXYZABCDEFGHIJKLM"
+        ), f"Expected EFGHIJKLMNOPQRSTUVWXYZABCDEFGHXJKLMNOPQRSTUVWXYZABCDEFGHIJKLM, got {result}"
+
+    def test_multi_deletion(alphabet_sequence):
+        test_row = create_test_row("c.3_6del", alphabet_sequence)
+        result = create_mutant_sequence(test_row, deletion_mutation)
+        assert result == "ABG", f"Expected ABG, got {result}"
+
+    def test_multi_deletion_long(alphabet_sequence_long):
+        test_row = create_test_row("c.35_37del", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, deletion_mutation)
+        assert (
+            result == "EFGHIJKLMNOPQRSTUVWXYZABCDEFGHLMNOPQRSTUVWXYZABCDEFGHIJKLMNO"
+        ), f"Expected EFGHIJKLMNOPQRSTUVWXYZABCDEFGHLMNOPQRSTUVWXYZABCDEFGHIJKLMNO, got {result}"
+
+    def test_single_deletion(alphabet_sequence):
+        test_row = create_test_row("c.3del", alphabet_sequence)
+        result = create_mutant_sequence(test_row, deletion_mutation)
+        assert result == "ABDEFG", f"Expected ABDEFG, got {result}"
+
+    def test_multi_delins(alphabet_sequence):
+        test_row = create_test_row("c.3_6delinsXYZ", alphabet_sequence)
+        result = create_mutant_sequence(test_row, delins_mutation)
+        assert result == "ABXYZG", f"Expected ABXYZG, got {result}"
+
+    def test_multi_delins_long(alphabet_sequence_long):
+        test_row = create_test_row("c.35_37delinsWXYZ", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, delins_mutation)
+        assert (
+            result == "EFGHIJKLMNOPQRSTUVWXYZABCDEFGHWXYZLMNOPQRSTUVWXYZABCDEFGHIJKLMNO"
+        ), f"Expected EFGHIJKLMNOPQRSTUVWXYZABCDEFGHWXYZLMNOPQRSTUVWXYZABCDEFGHIJKLMNO, got {result}"
+
+    def test_single_delins(alphabet_sequence):
+        test_row = create_test_row("c.3delinsXYZ", alphabet_sequence)
+        result = create_mutant_sequence(test_row, delins_mutation)
+        assert result == "ABXYZDEFG", f"Expected ABXYZDEFG, got {result}"
+
+    def test_ins(alphabet_sequence):
+        test_row = create_test_row("c.3_4insXYZ", alphabet_sequence)
+        result = create_mutant_sequence(test_row, insertion_mutation)
+        assert result == "ABCXYZDEFG", f"Expected ABCXYZDEFG, got {result}"
+
+    def test_ins_long(alphabet_sequence_long):
+        test_row = create_test_row("c.35_36insPPP", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, insertion_mutation)
+        assert (
+            result == "EFGHIJKLMNOPQRSTUVWXYZABCDEFGHIPPPJKLMNOPQRSTUVWXYZABCDEFGHIJKLM"
+        ), f"Expected EFGHIJKLMNOPQRSTUVWXYZABCDEFGHIPPPJKLMNOPQRSTUVWXYZABCDEFGHIJKLM, got {result}"
+
+    def test_multi_dup(alphabet_sequence):
+        test_row = create_test_row("c.3_5dup", alphabet_sequence)
+        result = create_mutant_sequence(test_row, duplication_mutation)
+        assert result == "ABCDECDEFG", f"Expected ABCDECDEFG, got {result}"
+
+    def test_multi_dup_long(alphabet_sequence_long):
+        test_row = create_test_row("c.39_42dup", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, duplication_mutation)
+        assert (
+            result
+            == "IJKLMNOPQRSTUVWXYZABCDEFGHIJKLMNOPMNOPQRSTUVWXYZABCDEFGHIJKLMNOPQRST"
+        ), f"Expected IJKLMNOPQRSTUVWXYZABCDEFGHIJKLMNOPMNOPQRSTUVWXYZABCDEFGHIJKLMNOPQRST, got {result}"
+
+    def test_single_dup(alphabet_sequence):
+        test_row = create_test_row("c.5dup", alphabet_sequence)
+        result = create_mutant_sequence(test_row, duplication_mutation)
+        assert result == "ABCDEEFG", f"Expected ABCDEEFG, got {result}"
+
+    def test_inv(alphabet_sequence):
+        test_row = create_test_row("c.3_4inv", alphabet_sequence)
+        result = create_mutant_sequence(test_row, inversion_mutation)
+        assert result == "ABDCEFG", f"Expected ABDCEFG, got {result}"
+
+    def test_inv_long(alphabet_sequence_long):
+        test_row = create_test_row("c.41_42inv", alphabet_sequence_long)
+        result = create_mutant_sequence(test_row, inversion_mutation)
+        assert (
+            result == "KLMNOPQRSTUVWXYZABCDEFGHIJKLMNPOQRSTUVWXYZABCDEFGHIJKLMNOPQRST"
+        ), f"Expected KLMNOPQRSTUVWXYZABCDEFGHIJKLMNPOQRSTUVWXYZABCDEFGHIJKLMNOPQRST, got {result}"
diff --git a/tests/test_mutation_reference_creation.py b/tests/test_mutation_reference_creation.py
new file mode 100644
index 00000000..f8f0c0f8
--- /dev/null
+++ b/tests/test_mutation_reference_creation.py
@@ -0,0 +1,38 @@
+import pytest
+import os
+import src.cosmic_cancer_mutation_fasta_creation
+import sys
+from pdb import set_trace as st
+
+# # uncomment these if I move global variables back into the main function in src.cosmic_cancer_mutation_fasta_creation
+# intronic_mutations = 0
+# posttranslational_region_mutations = 0
+# unknown_mutations = 0
+# uncertain_mutations = 0
+# ambiguous_position_mutations = 0
+# cosmic_incorrect_wt_base = 0
+
+# List of argument combinations to test
+argument_combinations = [
+    ('cancer', '37'),
+    ('cell_line', '38'),
+    ('census', '38'),
+    ('resistance', '38'),
+    ('screen', '38'),
+]
+
+@pytest.mark.parametrize("mutation_class_name, grch_number", argument_combinations)
+def test_mutation_reference_creation_main(mutation_class_name, grch_number, monkeypatch):
+    # Set the TEST_MODE environment variable
+    monkeypatch.setenv('TEST_MODE', '1')
+    
+    # Mock command-line arguments
+    monkeypatch.setattr(sys, 'argv', ['my_script.py', '-m', mutation_class_name, '-g', grch_number])
+    
+    src.cosmic_cancer_mutation_fasta_creation.main()
+
+    # access one of the global variables with referencing eg src.cosmic_cancer_mutation_fasta_creation.intronic_mutations
+
+    # st()
+    
+# run with pytest tests/test_mutation_reference_creation.py -v -s
\ No newline at end of file
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 05603945..1db54ae2 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -9,6 +9,7 @@
     find_latest_ens_rel,
     search_species_options,
     ref_species_options,
+    read_fasta,
 )
 
 from gget.constants import UNIPROT_REST_API, ENSEMBL_REST_API, ENSEMBL_FTP_URL_NV
@@ -23,6 +24,24 @@
 
 
 class TestUtils(unittest.TestCase):
+    def test_read_fasta(self):
+        result_to_test1, result_to_test2 = read_fasta(
+            "tests/fixtures/muscle_nt_test.fa"
+        )
+        result_to_test = result_to_test1 + result_to_test2
+
+        expected_result1 = [
+            ">ENSTGUT00000006367 primary_assembly:bTaeGut1_v1.p:1:107513786:107528106:-1",
+            ">ENSTGUT00000027003 primary_assembly:bTaeGut1_v1.p:1:107513786:107526965:-1",
+        ]
+        expected_result2 = [
+            "TAGTTGTGATACTGGTCAGTGTCACTGGGGCTGTTTCACCCTGTCATTTGAAGACTTCCATTT",
+            "CTTTGTTGGAGGCTTTTTGCTGGGGCTGGCGTCGTAAGGACCCGAACATCTCTTCCAGATGGGAT",
+        ]
+        expected_result = expected_result1 + expected_result2
+
+        self.assertEqual(result_to_test, expected_result)
+
     def test_n_colors(self):
         result_to_test = n_colors("A")
         expected_result = "\x1b[38;5;15m\x1b[48;5;9mA\x1b[0;0m"