Thus the codon tag is clearly very different from aaa, it is as if tag roughly marks the boundaries of regions, and aaa is just scattered throughout. Are any other codons similar to tag? The frequency analysis blurs counts so more subtle differences not as obvious have to run gap counter for each to directly see, and how to easily “see”? Notice how the tag distribution has a long tail. The gap bins only go to 9 in the figures, but for the full dataset the last nonzero gap bin for “tag” is at a remarkable bin 70. For “aaa” the last nonzero bin is much earlier, at bin number 23 (even though there are 10 times as many (aaa) codons as (tag) codons). For “taa” the last nonzero bin is at 60, while for tga it is at 53. The codons taa, tga, and tag are known as the stop codons and the gaps between them are known as ORFs, or ORFs. A subtlety in the statistical analysis is that the stop codons do not have to match to define such anomalously large regions (according to observation). Thus, a biochemical encoding scheme must exist that works with any of the three stop codons seen as equivalent, thus the naming for this group as “stop” codons (and their grouping as such in Figure 3.2). For more nuances of the naming convention “stop” codon when delving into the encoding biochemistry see [1, 3].
Figure 3.2 ORF encoding structure is revealed in the V. cholera genome by gaps between stop codons in the genomic sequence. X‐axis shows the size of the gap in codon count between reference codons (stops for conventional ORFs, or 3com set for comparisons in table), Y‐axis shows the counts.
3.3 ORF Discovery from Long‐Tail Distribution Anomaly
Once codon grouping is revealed, where a frequency analysis on codons on the stop codons (TAA, TAG, TGA) shows they are rare. Focusing on the stop codons it is easily found that the gaps between stop codons can be quite anomalous compared to the gaps between other codons (see prog2.py addendum 6):
------------------- prog2.py addendum 6 --------------------- # need gap stats between codons, the stop codon group (orf_finder), and # the 'common' reference group (corf_finder): def orf_finder ( seq, frame ): gapcounts = {} edgecounts = {} pattern = '[acgtACGT]' result = re.findall(pattern, seq) seqlen = len(seq) output_fh = open("orf_output", 'w') output_fh.close() oldindex=0 oldcodon="" for index in range(frame,seqlen-2): rem = (index+3-frame)%3 if rem!=0: continue codon = result[index]+result[index+1]+result[index+2] if (codon!="TAA" and codon!="TAG" and codon!="TGA"): continue else: gap = index - oldindex if gap%3!=0: print "gap=", gap, "index=", index break quant = 100 bin = gap/quant if oldindex!=0: if bin in gapcounts: gapcounts[bin]+=1 else: gapcounts[bin]=1 if oldcodon!="": edge=oldcodon + codon if edge in edgecounts: edgecounts[edge]+=1 else: edgecounts[edge]=1 slice = result[oldindex: index+2+1] output_fh = open("orf_output", 'a') slicejoin = "" slicejoin = slicejoin.join(slice) orfline = slicejoin + '\n' output_fh.write(orfline) oldindex=index oldcodon=codon npcounts = np.empty((0)) for i in sorted(gapcounts): npcounts = np.append(npcounts,gapcounts[i]+0.0) print "gapbin", i, "count =", gapcounts[i] ecounts = np.empty((0)) for i in sorted(edgecounts): ecounts = np.append(ecounts,edgecounts[i]+0.0) print "edgecodon", i, "count =", edgecounts[i] probs = count_to_freq(npcounts) #usage: orf_finder(EC_sequence,0) # def corf_finder ( seq, frame ): # same except not 'stop' boundariy condition but 'common': # if (codon!="AAA" and codon!="GAA" and codon!="GAT") #usage: #corf_finder(EC_sequence,0) --------------- prog2.py addendum 6 end ---------------------
ORFs are “open reading frames,” where the reference to what is open is lack of encounter with a stop codon when traversing the genome with a particular codon framing, e.g. ORFs are regions devoid of stop codons when traversed with the codon framing choice of the ORF. When referring to ORFs in most of the analysis we refer to ORFs of length 300 bases or greater. The restriction to larger ORFs is due to their highly anomalous occurrences and likely biological encoding origin (see Figure 3.2), e.g. the long ORFs give a strong indication of containing the coding region of a gene. By restricting to transcripts with ORFs >= 300 in length we have a resulting pool of transcripts that are mostly true coding transcripts.
The above example shows a bootstrap finite state automaton (FSA) process on genomic data: first scan through the genomic data base‐by‐base and obtain counts on nucleotide pairs with different gap sizes between the nucleotides observed [1, 3]. This then allows a mutual information analysis on the nucleotide pairs taken at the different gap sizes. What is found for prokaryotic genomes (with their highly dense gene placement), is a clear signal indicating anomalous statistical linkages on bases three apart [1, 3]. What is discovered thereby is codon structure, where the coding information comes in groups of three bases. Knowing this, a bootstrap analysis of the 64 possible 3‐base groupings can then be done, at which point the anomalously low counts on “stop” codons is then observed. Upon identification of the stop codons their placement (topology) in the genome can then be examined and it is found that their counts are anomalously low because there are large stretches of regions with no stop codon (e.g. there are stop codon “voids,” known as “ORFs”). The codon void topologies are examined in a comparative genomic analysis in [1, 3]. As noted previously, the stop codons, which should occur every 21 codons on average if DNA sequence data was random, are sometimes not seen for stretches of several hundred codons (see Figure 3.2).
Not surprisingly, longer genes stand out clearly in this process, since their anomalous, clearly nonrandom DNA sequence, is being maintained as such, and not randomized by mutation, (as this would be selected against in the survival of the organism that is dependent on the gene revealed).
The preceding basic analysis can provide a gene‐finder on prokaryotic genomes that comprises a one‐page Python script that can perform with 90–99% accuracy depending on the prokaryotic genome. A second page of Python coding to introduce a “filter,” along the lines of the bootstrap learning process mentioned above, leads to an ab initio prokaryotic gene‐predictor with 98.0–99.9% accuracy. Python code to accomplish this is shown in what follows (Chapter 4). In this process, all that is used is the raw genomic data (with its highly structured intrinsic statistics) and methods for identifying statistical anomalies and informatics structural anomalies: (i) anomalously high mutual information is identified (revealing codon structure); (ii) anomalously high (or low) statistics on an attribute or event is then identified (low stop codon counts, lengthy stop codon voids); then anomalously high sub‐sequences (binding site motifs) are found in the neighborhood of the identified ORFs (used in the filter).
3.3.1 Ab initio Learning with smORF’s, Holistic Modeling, and Bootstrap Learning
In work on prokaryotic gene prediction (V. cholera in what follows), a program (smORF) was developed for an extended ORF characterization (to characterize “some more ORFs” with different trinucleotide delimiters than stops). Using that software with a simple start‐of‐coding heuristic it was possible to establish good gene prediction for ORFs of length greater than 500 nucleotides. The smORF gene identification was used in a bootstrap gene‐annotation process (where no initial training data was provided). Part of the functionality for smORF is encompassed in prog2.py program described thus far. The strength of the gene identification was then improved