Codetta is a Python program for predicting the genetic code (codon table) of an organism from nucleotide sequence data.
The analysis consists of three steps:
- Aligning the input nucleotide sequence to a database of profile Hidden Markov models (HMMs) of proteins (such as the Pfam database)
- Collating the resulting alignments into a single output file
- Inferring the genetic code from the alignment output file
A detailed explanation of the underlying probability model can be found in Shulgina & Eddy (2021). If you are looking to reproduce results from Shulgina & Eddy (2021), please follow the README for Codetta v1.0.
If you encounter any problems in the installation or usage of Codetta, please leave a Github issue or email me at [email protected]
Codetta was developed for Python versions 3.5+ on Linux and MacOS. Required
Python packages are numpy and scipy.
Once inside the Codetta directory, run the setup script
./setup.sh
This script will check Codetta requirements and set up a local version of
HMMER. If you are missing wget or gzip, on MacOS, you can use
brew install. (For Linux, you'll have to use your system's package
management tool.) Note that on MacOS, you'll also need XCode installed so
that HMMER can compile.
Codetta aligns the input sequence to a database of protein profile HMMs.
By default, Codetta will use the Pfam database. You can download a version of
Pfam 35.0 specially built for Codetta from our website. From the codetta
directory, use these commands to download and uncompress it (note: this will
take about 3 Gb of disk space)
cd resources/
wget http://eddylab.org/publications/Shulgina21/Pfam-A_enone.tar.gz
tar xf Pfam-A_enone.tar.gz
rm Pfam-A_enone.tar.gz
cd ..
If you want to build your own custom profile HMM database, this is be described in a later section.
The program codetta.py rolls the three main analysis steps into a single
script. Usage for codetta.py is
./codetta.py [input sequence file] [optional arguments]
The three steps can also be run separately using the programs:
codetta_align.py: Align profile HMMs to the input nucleotide sequence.codetta_summary.py: Summarize profile HMM alignments into an alignment output file.codetta_infer.py: Infer the genetic code from the alignment output file.
General usage for these programs is
./[program name] [input file or prefix] [optional arguments]
For any of these programs, type ./[program name] --help for complete usage
details.
The simplest way to run Codetta is by using the codetta.py program.
This program performs the three analysis steps in order. All you have to do is specify a nucleotide sequence input file. This file should contain nucleotide sequences from a single organism in FASTA format. This can be a genome, transcriptome, collection of genes, etc.
Make sure you have first set up a profile HMM database (see above)!
In the examples/ directory, you will find the file GCA_014211875.1.fna
which contains the genome of the bacterium Nasuia deltocephalinicola.
We can predict the genetic code of this bacterium simply by running
./codetta.py examples/GCA_014211875.1.fna
You will see outputs written to the terminal indicating that each of the three Codetta steps is executing. Codetta will create five files in the directory containing the input sequence file.
- Processed sequence file:
examples/GCA_014211875.1.fna.sequence_pieces.fna - Preliminary translation file + ssi index:
examples/GCA_014211875.1.fna.preliminary_translation.faa - Alignment output file:
examples/GCA_014211875.1.fna.Pfam-A_enone.hmm.alignment_output.txt - Inference output file:
examples/GCA_014211875.1.fna.Pfam-A_enone.hmm.1e-10_0.9999_0.01_excl-mtvuy.genetic_code.out
At the end, the inferred genetic code is printed to the terminal
Genetic code: FFLLSSSSYY??CCWWL?LLPPPPHHQQ????I?IMTTTTNNKKSSRRVVVVAAAADDEEGGGG
This corresponds to the inferred translation of each of the 64 codons, in order of 'UUU, UUC, UUA, UUG, UCU, UCC, ..., GGA, GGG' (iterating 3rd, 2nd, then 1st base through UCAG).
The ?s correspond to codons that had no inferred amino acid meaning-- this means that there was insufficient or ambiguous information about the codon to make a confident inference. This is also the expected inference for stop codons (since Codetta does not explicitly predict stop codons).
A detailed overview of the analysis can be found at
examples/GCA_014211875.1.fna.Pfam-A_enone.hmm.1e-10_0.9999_0.01_excl-mtvuy.genetic_code.out.
The long file extension specifies the inference parameters. You can specify
your own (simpler) name for the inference output file using the
--inference_output argument.
This file contains a detailed overview of the genetic code inference results:
# Analysis arguments
alignment_prefix examples/GCA_014211875.1.fna
profile_database Pfam-A_enone.hmm
results_summary None
evalue_threshold 1e-10
prob_threshold 0.9999
max_fraction 0.01
excluded_pfams mtvuy
#
# Codon inferences Consensus columns
# codon inference std code diff? N aligned N used
TTT F F N 546 546
TTC F F N 48 48
TTA L L N 703 703
TTG L L N 82 82
TCT S S N 307 307
...
GGA G G N 248 248
GGG G G N 73 72
#
# Log decoding probabilities
# codon logP(A) logP(C) logP(D) logP(E) logP(F) logP(G) ...
TTT -1250.3530 -1453.8418 -2107.4121 -1850.9243 0.0000 -2034.8423 ...
...
#
# Final genetic code inference
FFLLSSSSYY??CCWWL?LLPPPPHHQQ????I?IMTTTTNNKKSSRRVVVVAAAADDEEGGGG
You might choose to change some parameters of the analysis. Some commonly used options are
- Use the
-pargument to specify a different profile HMM database. Note that the database must be located in theresources/directory. - The
excluded_pfamsline in the above output file refers to which groups of problematic Pfam domains are excluded from the analysis (mmitochondrial,ttransposon and other mobile genetic element,vviral,uselenocysteine-containing, andypyrrolysine-containing). By default, all of these groups are excluded. Use-m -t -v -u -yto include these groups. For instance, if you're analyzing a mitochondrial genome you will want to use-mto include Pfams commonly found in mitochondrial genomes. Likewise, use the-vand-tflags when analyzing a viral genome. - You can specify your own name for the inference output file with the
--inference_outputargument. - Use
-eto change the profile HMM hit e-value threshold (default is 1e-10). Use-rto change the probability threshold to call an amino acid meaning for a codon (default is 0.9999). Use-fto change the maximum fraction of observations for a codon coming from a single profile HMM position (default is 0.01). - If you're running several Codetta analyses, you can use
--results_summaryto specify a file to which a one-line summary of the results will be appended to.
See the full list of options with ./codetta.py --help
The basic usage of Codetta is to use codetta.py, which rolls the three main
steps of the analysis into a single program. See the "Starting out..."
section above for an example and description of the output.
If, for any reason, you want to run the three Codetta analysis steps
separately, you can use the programs codetta_align.py, codetta_summary.py
, and codetta_infer.py.
The example in the previous section
./codetta.py examples/GCA_014211875.1.fna
can be also run by breaking the steps up as
./codetta_align.py examples/GCA_014211875.1.fna
./codetta_summary.py examples/GCA_014211875.1.fna
./codetta_infer.py examples/GCA_014211875.1.fna
You might want to do this if:
- You want to infer the genetic code multiple times using different
parameters without having to re-align the profile HMMs to your input
sequence. In this case, you would run
codetta_align.pyandcodetta_summary.pyonce and thencodetta_infer.pyseveral times with different parameters. - You want to parallelize the computationally intensive alignment step on a computing cluster. See the next section.
Codetta is default set up to run locally. However, the codetta_align step
can start to take a long time for longer genomes. For instance, a 4 Mb
bacterial genome takes about 10 minutes on a MacBook Pro (3.5 GHz
processor). If you're analyzing a large genome or many sequences, we
recommended parallelizing the codetta_align step on a computing cluster.
This is simple for clusters using a SLURM job scheduler. Install Codetta
following the usual steps. Then, manually edit the
resources/template_jobarray.sh file to have the right parameters for your
computing cluster. The template_jobarray.sh file looks like this:
#!/bin/bash
#SBATCH -p [SPECIFY PARTITION NAME] # partition
#SBATCH --time=8:30:00 # wall-clock time (mins:secs)
#SBATCH -c 1 # requesting 1 core
#SBATCH --mem=4000M
#SBATCH -o /dev/null # file to which STDOUT + STDERR will be written
You'll probably only need to specify a partition name.
Then, when you run the codetta_align step, add the argument
--parallelize_hmmscan 's'. This will cause the hmmscan jobs to be sent as
a SLURM job array to the specified partition instead of being executed
locally. Make sure to wait for the SLURM jobs to finish before proceeding to
codetta_summary. That's it!
If your cluster uses a different job scheduler, you will have to manually
modify the job array template file and the code in codetta.py
(processing_genome function) that writes and sends the job array file.
Pfam domains are expected to align to about ~50% of coding sequence. If you're focused on a specific clade, you could increase the proportion of aligned sequence by using a customized profile HMM database. Starting from a set of multiple sequence alignments of expected proteins in the organism of interest, you can build a custom profile HMM database.
In this example, we will build a profile HMM database tailored for metazoan
mitochondria which have only 13 protein-coding genes. In the
examples/mito-db/ directory, you can find multiple sequence alignment files
in Stockholm format for each of the 13 genes. *
The first step is to use hmmbuild to create profile HMMs from each of the
alignment files.
cd examples/mito-db/
ls metazoan_mito*.msa | xargs -I {} ../../hmmer-3.3.2/bin/hmmbuild --enone {}.hmm {}
Then we concatenate all of these .hmm files into a single database and
finally run hmmpress to finish
cat metazoan_mito*.hmm > metazoan_mito_proteins.hmm
../../hmmer-3.3.2/bin/hmmpress metazoan_mito_proteins.hmm
mv metazoan_mito_proteins.hmm* ../../resources/.
At the end, the profile HMM database is moved into the resources directory,
where Codetta expects to find all profile HMM databases. When running a
Codetta analysis, you will have to specify the custom profile HMM database
with -p [name of profile HMM db].
*These alignments were created by searching the human mitochondrial protein
sequences against SwissProt with jackhmmer
We have also provided a simple program for a downloading FASTA file from GenBank by specifying either a genome assembly accession or a nucleotide accession.
codetta_download: Download a genome assembly or nucleotide sequence from GenBank
We can download the Nasuia deltocephalinicola genome using the Genbank genome assembly accession by
./codetta_download.py GCA_014211875.1 a --sequence_file examples/GCA_014211875.1.fna
This will download a FASTA file containing the GCA_014211875.1 sequence into
examples/GCA_014211875.1.fna. The argument a specifies that this is an
assembly database accession and not a nucleotide accession (which would be
c).
We can download the mitochondrial genome of the green algae Pycnococcus provasolii, which is under NCBI nucleotide accession NC_013935.1 by
./codetta_download.py NC_013935.1 c --sequence_file examples/NC_013935.1.fna
Notice the argument c specifying that this is a nucleotide database
accession.
Now let's pull it all together by predicting the genetic code of the P. provasolii mitochondrial genome:
./codetta_download.py NC_013935.1 c --sequence_file examples/NC_013935.1.fna
./codetta.py examples/NC_013935.1.fna -m
Notice how we specified the -m argument to indicate that we do not want to
exclude Pfam domains associated with mitochondrial genomes, which is expected to have two codon reassignments.
Alternatively, we could also run the same analysis as
./codetta_download.py NC_013935.1 c --sequence_file examples/NC_013935.1.fna
./codetta_align.py examples/NC_013935.1.fna --align_output examples/Pprovasolii_mito
./codetta_summary.py examples/Pprovasolii_mito
./codetta_infer.py examples/Pprovasolii_mito -m --inference_output examples/Pprovasolii_mito_Pfam_genetic_code.out
Here we are also showing how to use the --align_output argument to write
the alignment output files with the more informative prefix Pprovasolii_mito
and the --inference_output argument to write the inference output file
examples/Pprovasolii_mito_Pfam_genetic_code.out.
The output genetic code is:
FF??SSSSYY??CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSR?V?VVAAAADDEEGGGG
Comparing to the standard genetic code (below), you can see that two codons have alternative meanings: the stop codon UGA is now tryptophan codon and the isoleucine codon AUA is now a methionine codon. Some codons are uninferred (?) due to few aligned Pfam consensus columns (look at the inference output file for more detail).
P mt code : FF??SSSSYY??CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSR?V?VVAAAADDEEGGGG
std code : FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
^ ^
This alternative genetic code in Pycnococcus mitochondria is summarized in Noutahi et al (2019).
Occasionally, Codetta will make an incorrect codon inference. This usually occurs when an assumption made by the method is broken.
Common situations to look out for:
- Stop codons erroneously aligning to profile HMMs. A common cause of this are pseudogenes that maintain enough homology for profile HMM alignment, but have acquired in-frame stop codons, causing alignment between stop codons and consensus columns of a profile HMM. This may lead to stop codons erroneously inferred to code for an amino acid. This could also occur due to profile HMMs aligning over intron boundaries, and stop codons that are readthrough / frameshifted over (often seen in viruses).
- The input sequence contaminated by another organism. This can lead to incorrect genetic code prediction if the contaminating organism uses a different genetic code. A common source of such contamination are Mycoplasma bacteria.
- The input nucleotide sequence not being representative of the mature mRNA coding sequence. In some mitochondrial genomes (plants, trypanosomes, etc), there is pervasive mRNA editing such that some codons in the genome are systematically converted to another codon in the mature mRNA. This may lead to incorrect codon inference by Codetta.
- The composition of the input sequence significantly deviating from the profile HMMs used. If there is a significant difference in amino acid composition, it may appear that a codon has systematically substituted into sites conserved for another amino acid in the profile HMMs. To Codetta, this resembles the signature of codon reassignment, leading to incorrect inferences. A way around this is to use a custom profile HMM database matched for amino acid composition / GC content.
To further validate a predicted genetic code, some approaches are:
- Try predicting the genetic code again with different Codetta parameters or with a different profile HMM database. The predicted genetic code should be robust to such changes.
- Genes for translational components. Does the genome encode tRNAs / release factors consistent with the predicted genetic code?
- High-quality multiple sequence alignments of single-copy conserved genes. Does the codon occur at positions where the predicted amino acid occurs in outgroups?
- Experimental validation. If possible, proteomic mass spectrometry or other methods that show that the predicted amino acid is inserted at codon positions would provide definitive validation.
Y. Shulgina, S. R. Eddy. (2023). Codetta: predicting the genetic code from nucleotide sequence. Bioinformatics, 39:btac802. doi: 10.1093/bioinformatics/btac802
Y. Shulgina, S. R. Eddy. (2021). A computational screen for alternative genetic codes in over 250,000 genomes. eLife, 10:e71402. doi: 10.7554/eLife.71402
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent