Day 3: Open Reading Frame (ORF) Identification & In Silico Gene Prediction

1. Anatomical Principles of Open Reading Frames (ORFs)

An Open Reading Frame (ORF) is a contiguous stretch of codons in DNA or RNA that begins with a recognized translation initiation (start) codon, continues in-frame with triplet codons encoding amino acids, and terminates at a translation termination (stop) codon, containing zero intermediate stop codons.

2. Distinguishing Genuine Genes from Spurious Random ORFs

In any random sequence of non-coding DNA, a stop codon occurs on average once every $64 / 3 \approx 21.3$ codons ($~64\text{ bp}$). Consequently, random non-coding DNA produces thousands of short, meaningless ORFs. Computational tools apply rigorous statistical filters to separate authentic protein-coding genes from background noise:

• Length Cutoff Threshold: Spurious random ORFs rarely exceed 100 codons (300 bp). Applying a minimum threshold of $\ge 75-100\text{ amino acids}$ eliminates over 95% of random background ORFs.
• Ribosome Binding Site (Shine-Dalgarno in Prokaryotes): In bacteria, genuine start codons are preceded 6–10 nucleotides upstream by the polypurine consensus sequence $5'\text{-AGGAGG-}3'$ (the Shine-Dalgarno sequence), complementary to the 3' end of 16S rRNA.
• Kozak Consensus Context (Eukaryotes): In eukaryotes, efficient translation initiation requires the Kozak sequence context: $$\text{5'-(gcc)gccRccAUGG-3'}$$ where $R$ is a purine (A or G, most crucially an A at position $-3$ and a G at position $+4$ relative to the A of the ATG).
• Codon Usage Bias & Hexamer Frequencies: Authentic coding sequences exhibit distinct synonymous codon preferences matching the organism's tRNA pool, quantified by the Codon Adaptation Index (CAI) and Hexamer frequency scoring (Markov models).

3. Production Python Script: Six-Frame ORF Extractor & Peptide Classifier

#!/usr/bin/env python3
"""
Day 3 Practical: In Silico Open Reading Frame (ORF) Scanner
Detects all forward and reverse-strand ORFs exceeding length thresholds.
"""

from Bio.Seq import Seq
import re

def find_orfs_in_strand(dna_seq: Seq, strand="+", min_protein_len=30, table=1):
    """Identify all ORFs starting with ATG and ending at standard stop codons."""
    seq_str = str(dna_seq).upper()
    seq_len = len(seq_str)
    orfs = []
    
    # Check all three frames on this strand
    for frame in range(3):
        # Scan for ATG start codons
        for start_match in re.finditer(r"ATG", seq_str[frame:]):
            start_pos = frame + start_match.start()
            
            # Translate from this ATG in steps of 3 until stop codon
            subseq = dna_seq[start_pos:]
            # Trim to full codons
            trimmed = subseq[:len(subseq) - (len(subseq) % 3)]
            protein = str(trimmed.translate(table=table))
            
            # Find first in-frame stop codon (*)
            stop_idx = protein.find("*")
            if stop_idx != -1:
                orf_protein = protein[:stop_idx]
                orf_dna_len = (stop_idx + 1) * 3
                end_pos = start_pos + orf_dna_len
                
                if len(orf_protein) >= min_protein_len:
                    orfs.append({
                        "strand": strand,
                        "frame": frame + 1 if strand == "+" else -(frame + 1),
                        "start_nt": start_pos + 1,
                        "end_nt": end_pos,
                        "length_nt": orf_dna_len,
                        "length_aa": len(orf_protein),
                        "protein_seq": orf_protein
                    })
    return orfs

def scan_all_orfs(dna_string: str, min_aa=40):
    """Scan both forward and reverse-complement strands for valid ORFs."""
    dna = Seq(dna_string.upper())
    forward_orfs = find_orfs_in_strand(dna, strand="+", min_protein_len=min_aa)
    rev_dna = dna.reverse_complement()
    reverse_orfs = find_orfs_in_strand(rev_dna, strand="-", min_protein_len=min_aa)
    
    all_orfs = forward_orfs + reverse_orfs
    # Sort descending by amino acid length
    all_orfs.sort(key=lambda x: x["length_aa"], reverse=True)
    return all_orfs

if __name__ == "__main__":
    # Test on full-length synthetic construct
    test_sequence = (
        "GCGCGCGCAAAATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCAGCCGCAGCCTTT"
        "GTGAACCAACACCTGTGCGGCTCACACCTGGTGGAAGCTCTCTACCTAGTGTGCGGGGAACGAGGCTTCTTCTACACACCCAAGACC"
        "CGCCGGGAGGCAGAGGACCTGCAGGTGGGGCAGGTGGAGCTGGGCGGCGGCCCTGGTGCAGGCAGCCTGCAGCCCTTGGCCCTGGAG"
        "GGGTCCCTGCAGAAGCGTGGCATTGTGGAACAATGCTGTACCAGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAACTAGCCC"
        "CCCGGGATCGATCGATCGATC"
    )

    detected = scan_all_orfs(test_sequence, min_aa=30)
    print(f"Total ORFs Detected (>=30 aa): {len(detected)}\n")
    for i, orf in enumerate(detected, 1):
        print(f"ORF #{i} | Strand: {orf['strand']} | Frame: {orf['frame']} | Length: {orf['length_aa']} AA ({orf['length_nt']} nt)")
        print(f"Coordinates: {orf['start_nt']}..{orf['end_nt']}")
        print(f"Protein: {orf['protein_seq'][:40]}...\n")