omics.co.in

6: Object-Oriented Programming

🐍 Downloadable Lesson Python Script 6.4 KB

ch06_dna_oop.py

DnaSequence, ProteinSequence, and GenomicVariant classes implementing object inheritance, molecular weight calculations, and variant mutagenesis.

$ python ch06_dna_oop.py

Chapter 6: Object-Oriented Programming (OOP) in Python for Biologists

In computational biology, biological entities possess both state (e.g. sequence string, accession ID, chromosome coordinates, organism) and behaviors (e.g. transcribing, reverse complementing, calculating GC content, translating). Object-Oriented Programming (OOP) allows us to encapsulate biological data with biological logic into reusable, modular Python classes.

Why OOP for Computational Biology?

Procedural scripts tracking loose string variables quickly become error-prone when processing thousands of genomic records. Encapsulating sequence state inside classes guarantees integrity, coordinates validation, and enables polymorphic pipelines across DNA, RNA, and protein models.

1. Classes (Defining Biological Blueprints)

A class defines the blueprint for a biological entity. The constructor method __init__ initializes instance variables and validates input sequence characters against standard IUPAC alphabets.

class DnaSequence:
    """Represents a double-stranded genomic DNA sequence."""
    
    def __init__(self, seq_id, sequence, organism="Homo sapiens"):
        self.seq_id = str(seq_id)
        self.sequence = sequence.upper().strip()
        self.organism = organism
        self.length = len(self.sequence)
        
        # Validate IUPAC DNA alphabet
        valid_bases = set("ATCGN")
        if not set(self.sequence).issubset(valid_bases):
            invalid = set(self.sequence) - valid_bases
            raise ValueError(f"Invalid DNA base(s) {invalid} in {self.seq_id}")

    def __repr__(self):
        return f"<DnaSequence {self.seq_id} | {self.length} bp | {self.organism}>"

2. Objects (Instantiating Real Genomic Records)

An object is an active instance of a class loaded into memory with its unique biological data.

# Instantiating gene objects
tp53 = DnaSequence("TP53_exon1", "ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGA", "Homo sapiens")
egfr = DnaSequence("EGFR_promoter", "CGCGGGAACAGCGTGCCCGGAGCCCG", "Homo sapiens")

print(tp53)
# Output: <DnaSequence TP53_exon1 | 33 bp | Homo sapiens>
print(f"Gene: {tp53.seq_id} | Length: {tp53.length} bp | Host: {tp53.organism}")

3. Methods (Encapsulating Biological Operations)

Methods are functions defined inside a class that operate directly on the object’s internal sequence data.

class BiologicalSequence(DnaSequence):
    
    def gc_content(self):
        """Calculates GC percentage of the sequence."""
        g = self.sequence.count("G")
        c = self.sequence.count("C")
        return round((g + c) / self.length * 100, 2) if self.length > 0 else 0.0

    def reverse_complement(self):
        """Generates the antiparallel reverse complement strand (5′ → 3′)."""
        trans = str.maketrans("ATCGNatcgn", "TAGCNtagcn")
        rc_seq = self.sequence.translate(trans)[::-1]
        return BiologicalSequence(f"{self.seq_id}_RC", rc_seq, self.organism)

    def transcribe(self):
        """Transcribes genomic DNA into messenger RNA (T -> U)."""
        return self.sequence.replace("T", "U")

    def translate(self):
        """Translates codons into single-letter amino acid peptide sequence."""
        codon_table = {
            "ATA":"I", "ATC":"I", "ATT":"I", "ATG":"M",
            "ACA":"T", "ACC":"T", "ACG":"T", "ACT":"T",
            "AAC":"N", "AAT":"N", "AAA":"K", "AAG":"K",
            "AGC":"S", "AGT":"S", "AGA":"R", "AGG":"R",
            "CTA":"L", "CTC":"L", "CTG":"L", "CTT":"L",
            "CCA":"P", "CCC":"P", "CCG":"P", "CCT":"P",
            "CAC":"H", "CAT":"H", "CAA":"Q", "CAG":"Q",
            "CGA":"R", "CGC":"R", "CGG":"R", "CGT":"R",
            "GTA":"V", "GTC":"V", "GTG":"V", "GTT":"V",
            "GCA":"A", "GCC":"A", "GCG":"A", "GCT":"A",
            "GAC":"D", "GAT":"D", "GAA":"E", "GAG":"E",
            "GGA":"G", "GGC":"G", "GGG":"G", "GGT":"G",
            "TCA":"S", "TCC":"S", "TCG":"S", "TCT":"S",
            "TTC":"F", "TTT":"F", "TTA":"L", "TTG":"L",
            "TAC":"Y", "TAT":"Y", "TAA":"*", "TAG":"*",
            "TGC":"C", "TGT":"C", "TGA":"*", "TGG":"W",
        }
        protein = []
        for i in range(0, self.length - 2, 3):
            codon = self.sequence[i:i+3]
            protein.append(codon_table.get(codon, "X"))
        return "".join(protein)

gene = BiologicalSequence("BRCA1_partial", "ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTCATTAATGCTATGCAGAAA")
print("GC Content:", gene.gc_content(), "%")
print("mRNA Transcript:", gene.transcribe())
print("Translated Peptide:", gene.translate())
print("Reverse Complement:", gene.reverse_complement())

4. Attributes & State Management

Attributes store essential metadata such as genomic coordinates, strand orientation, and Phred quality scores.

class FastqRead:
    """Encapsulates an Illumina or Nanopore sequencing read with Phred scores."""
    
    def __init__(self, read_id, sequence, quality_str):
        self.read_id = read_id
        self.sequence = sequence.strip()
        self.quality_str = quality_str.strip()
        # Decode Phred+33 ASCII scores to numeric integers
        self.phred_scores = [ord(char) - 33 for char in self.quality_str]

    @property
    def mean_quality(self):
        return sum(self.phred_scores) / len(self.phred_scores) if self.phred_scores else 0

    def is_high_quality(self, min_q=30):
        """Returns True if the read meets clinical Q30 benchmark."""
        return self.mean_quality >= min_q

read = FastqRead("@SRR123456.1", "GATCGATCGATCGATC", "IIIIIIIIIIIIIIII")
print(f"Read Mean Phred: Q{read.mean_quality:.1f} | High Quality (Q30): {read.is_high_quality(30)}")
← Return to Home Hub Scroll to Top ↑