Codon Bias

The GeneMarkS-T and GeneMarkS (7 (link)) algorithms share the following: (i) the heuristic method of initialization of the hidden semi-Markov model (HSMM) parameters (13 (link)), (ii) the Viterbi algorithm that finds maximum likelihood parse of transcript sequence into coding and non-coding regions and (iii) the concept of iterative self-training (7 (link)).
Important differences are as follows. Unlike rather homogeneous G + C content of prokaryotic genomes, variation in local G + C content across much longer eukaryotic genomes may reach 30–40%. It was shown that genomic sequence G+C content is one of the major factors driving the genome-wide pattern of codon usage (13 (link),14 (link)). Therefore, GeneMarkS-T attempts to group transcripts by G + C content (Figure 1). The number of groups (clusters) depends on how wide is the distribution of G + C composition of the whole set of transcripts. By adjusting cluster borders we place the same volume of transcript sequence into each cluster. The iterative self-training on sequences of each cluster runs similarly to that described for GeneMarkS (7 (link)). The procedure starts with initialization of the cluster-specific ‘heuristic’ model (13 (link)). Then rounds of (i) predictions of protein-coding regions, (ii) selecting a new set of sequences of predicted genes for training and iii/ re-estimation of parameters, follow until convergence, i.e. the set of predicted genes in the last iteration should be the same as in the previous iteration (Figure 1).
The total volume of transcript data may vary. If the input data is not large enough for self-training, the ‘heuristic’ parameters used in initialization (13 (link)) are accepted as the final set of parameters and predictions made with this parameter set are considered final. The rationale for this approach is the earlier demonstration that the ‘heuristic’ parameters give sufficiently accurate predictions of continuous protein-coding regions in short prokaryotic sequences, e.g. in metagenomic sequences (13 (link),15 (link)).
GeneMarkS-T derives in iterations the species-specific Kozak pattern, positional frequency model of the sequence near TIS (10 (link)). The frequencies are determined from the multiple alignment of 12 bp-long fragments surrounding predicted TISs with nucleotides A of start codons situated in position 7.
Recently introduced strand-specific RNA-Seq technology (16 (link)) determines which DNA strand served as a template for transcription. If this information is available GeneMarkS-T changes the HSMM architecture and eliminates states related to the non-transcribed DNA strand; this change reduces the rate of false positive predictions. In what follows GeneMarkS-T version with the strand specific HSMM is designated as GeneMarkS-T(S).
For genes predicted in each transcript GeneMarkS-T assigns log-odds scores computed as log of the ratio of probability of a sequence given the coding model to probability of the same sequence given the non-coding model. The distribution of lengths of protein coding region and non-coding sequences is taken into account; these distributions are modelled as the gamma distribution the exponential distribution, respectively (17 (link)).

To distinguish protein-coding sequences from the non-coding sequences, we extracted five features, i.e. the length and S-score of MLCDS, length-percentage, score-distance and codon-bias. The length and S-score of MLCDS were used as the first two features, which assess the extent and quality of the MLCDS, respectively (Supplementary Table S3). Moreover, as demonstrated earlier in the text, protein-coding transcripts possess a special reading frame obviously distinct from the other five in the distribution of ANT. We analyzed six MLCDS candidates outputted by dynamic programming of the six reading frames for each transcript, with the assumption that there must exist one best MLCDS (as described earlier in the text); however, this phenomenon does not generally exist for non-coding transcripts. Thus, we defined other two features, length-percentage and score-distance, as follows:

Where Ml is the length of the best MLCDS (according to S-score value) among that of six reading frames, and Y_i represents the length of each six of the MLCDS.

Where S is the S-score of the best MLCDS, and Ej represents the S-score of the other five MLCDS (Supplementary Table S3).
All aforementioned four selected features could, to some extent, distinguish the protein-coding and non-coding sequences and were concordantly higher in protein-coding transcripts and lower in non-coding transcripts (Supplementary Figure S4). Finally, we included the fifth feature, the frequency of single nucleotide triplets, in the MLCDS as the last feature to complement the construction of a classification model. This feature was defined as codon-bias, which evaluated the coding-non-coding bias for each of the 61 kinds of codons (the three stop codons were ruled out) (Supplementary Figure S5).
To get the positive and negative training sets, we extracted the five features for each best MLCDS from the known protein-coding and non-coding transcript data sets, respectively. We then incorporated these two training sets into a support vector machine (SVM) as a model construction (Figure 1c). We used the A Library for Support Vector Machines (LIBSVM) (13 ) to train an SVM model using the standard radial basis function kernel, where the C and gamma parameters were set by default.

Each of 20 amino acids in a protein can be encoded by 1 to 6 codons among a total of 61 sense codons. Each codon is made of 3 nucleotides. Therefore, the abundance of each amino acid in a protein is overall pertaining to the abundance of its codon(s) as a result of evolution and variation of the nucleotide sequence making up all these 61 codons. However, codon bias is widely present across diverse species [13 (link)], which would contribute to considerable variations from codon-based predictions. On the other hand, amino acid relative abundance is also affected by both metabolic (production) cost and amino acid decay rates [14 (link)]. For this study, the expected percentage of an amino acid in a protein was assessed by both the genetic code model and the proteome analysis [14 (link)]. The actual percentage of a residue in the protein is the number of this residue divided by total residue (amino acid) number of this protein. The degree of discrepancy between the expected abundance and the actual abundance indicates the bias or enrichment of a particular residue in a protein by evolution.

The aspartate protease gene was obtained from NCBI (National Center for Biotechnology Information) with the sequence number XP_001401093.1 and base-substituted using SnapGene 3.2.1 according to P. pastoris codon preference. The optimized gene (apa1) was synthesized by Sangon (Shanghai, China). The synthesized apa1 gene was amplified using PCR with forward primer (GCTCCAGCTCCAACTAGAAAG) and reverse primer (AGCTTGAGCAGCAAAACCC) specific to its sequence. The truncated apa1 gene without the signal peptide coding sequence was cloned into the pPICZαA vector using a one-step cloning ligation. The ligated product was identified by agarose gel electrophoresis and purified by a gel recovery kit. The pPICZαA/apa1 plasmid was validated by sequencing.

Evidence for evolutionary trends and conservation in specific organisms was derived from predicted protein sequences from 102 genomes, with broad coverage of eukaryotes and comprehensive coverage of Discoba lineages. This was supplemented with ten transcriptomes from Excavata lineages with poor genomic coverage, from which protein sequences were predicted using TransDecoder v5.5.0 (LongOrfs)⁸⁹ (Supplementary Table 4).
Orthologous groups were determined Orthofinder v2.3.12 (refs. ^{90 (link),91 (link)}) using default settings using diamond v2.0.5 and FastME 2.1.4. Reciprocal best BLAST (RBB) hits to T. brucei TREU927 proteins were identified using National Center for Biotechnology Information BLAST 2.9.0+ (ref. ⁹² ), reciprocal hits irrespective of forward and reverse search e-value and accepting reciprocal hits that identified a T. brucei gene with an identical sequence to the starting gene. For our analyses, a protein in a different species was defined as ‘the’ orthologue of a T. brucei gene if it was either the RBB or was the only orthogroup member in that species. We have used current bioinformatic approaches for all protein–protein orthologue analyses, but these are limited by the power of such computational comparative approaches.
Gain in complexity at different evolutionary distances was carried out using National Center for Biotechnology Information Taxonomy species classifications. We scored the proportion of proteins localizing to a particular organelle that had an orthologue in at least one species at that evolutionary distance (Supplementary Table 4) using a hypergeometric test to detect over-enrichment.
To determine the ratio of number of non-synonymous mutations (K_A) to synonymous mutations (K_S) RBB protein sequences were aligned using Clustal Omega^{93 (link)}. The corresponding coding sequences were mapped to the protein sequence alignment and scored for identical codons (no mutation), synonymous mutation, non-symonymous mutation or indel mutation (alignment gap). K_A/K_S was calculated per codon treating gaps as non-synonymous mutations. K_A/K_S was calculated without any codon bias correction for T. brucei brucei TREU927 against each Trypanozoon (African trypanosome) species for each reciprocal best BLAST orthologue, and averaged.