Python Scripts for Bioinformatics & Genomics

#!/usr/bin/env python3
"""
Chapter 1: Python for Next Gen Biologist
Script: ch01_dna_basics.py
Description: DNA sequence validation, GC content calculation, nucleotide frequencies,
             RNA transcription, and reverse complement generation.
Usage:
    python ch01_dna_basics.py --seq ATGCGATCGATCGATCGATAGCTAGCTA
    python ch01_dna_basics.py --file sequence.fasta
"""

import argparse
import sys

def validate_dna(sequence: str) -> bool:
    """Check if the sequence contains only valid IUPAC DNA bases (A, C, G, T, N)."""
    valid_bases = set("ACGTUN")
    return set(sequence.upper()).issubset(valid_bases)

def calculate_gc_content(sequence: str) -> float:
    """Calculate GC percentage of a DNA sequence."""
    seq = sequence.upper()
    g_count = seq.count('G')
    c_count = seq.count('C')
    total_len = len(seq)
    if total_len == 0:
        return 0.0
    return round(((g_count + c_count) / total_len) * 100, 2)

def transcribe(dna_sequence: str) -> str:
    """Transcribe DNA sequence into messenger RNA (T -> U)."""
    return dna_sequence.upper().replace('T', 'U')

def reverse_complement(sequence: str) -> str:
    """Generate the reverse complement of a DNA sequence (5' -> 3')."""
    complement_map = str.maketrans("ACGTNacgtn", "TGCANtgcan")
    return sequence.translate(complement_map)[::-1]

def nucleotide_counts(sequence: str) -> dict:
    """Return a dictionary of nucleotide counts."""
    seq = sequence.upper()
    return {
        'A': seq.count('A'),
        'C': seq.count('C'),
        'G': seq.count('G'),
        'T': seq.count('T'),
        'N': seq.count('N')
    }

def main():
    parser = argparse.ArgumentParser(description="Chapter 1: DNA sequence metrics and reverse complement.")
    parser.add_argument("-s", "--seq", help="Input DNA sequence directly as string")
    parser.add_argument("-f", "--file", help="Path to file containing DNA sequence")
    args = parser.parse_args()

    seq = None
    if args.seq:
        seq = args.seq.strip()
    elif args.file:
        with open(args.file, 'r') as f:
            lines = [line.strip() for line in f if not line.startswith('>')]
            seq = "".join(lines)
    else:
        # Default demo sequence: Human Beta-Globin (HBB) 5' CDS fragment
        seq = "ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGG"
        print("[!] No input provided. Using demo Human HBB CDS fragment.\n")

    if not validate_dna(seq):
        print("[ERROR] Sequence contains non-standard nucleotides.", file=sys.stderr)
        sys.exit(1)

    print("=== DNA Sequence Analysis Report ===")
    print(f"Sequence Length : {len(seq)} bp")
    print(f"GC Content      : {calculate_gc_content(seq)}%")
    counts = nucleotide_counts(seq)
    print(f"Base Frequencies: A={counts['A']} | C={counts['C']} | G={counts['G']} | T={counts['T']} | N={counts['N']}")
    print(f"mRNA Transcript : {transcribe(seq)[:50]}...")
    print(f"Rev-Complement  : {reverse_complement(seq)[:50]}...")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 2: Python for Next Gen Biologist
Script: ch02_data_structures.py
Description: Genomic data structures in Python: Lists, Tuples, Dictionaries, and Sets.
             Includes universal genetic codon table, restriction enzyme motifs, and k-mer sets.
Usage:
    python ch02_data_structures.py --seq ATGCGATCGATCGATCGATAGCTAGCTA
"""

import argparse

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',
}

RESTRICTION_SITES = {
    'EcoRI': 'GAATTC',
    'BamHI': 'GGATCC',
    'HindIII': 'AAGCTT',
    'NotI': 'GCGGCCGC',
    'XhoI': 'CTCGAG'
}

def translate_dna(dna_sequence: str) -> str:
    """Translates a DNA sequence into an amino acid peptide using codon dictionary."""
    seq = dna_sequence.upper()
    peptide = []
    for i in range(0, len(seq) - 2, 3):
        codon = seq[i:i+3]
        amino_acid = CODON_TABLE.get(codon, '?')
        if amino_acid == '_':
            peptide.append('*')  # Stop codon
            break
        peptide.append(amino_acid)
    return "".join(peptide)

def find_unique_kmers(sequence: str, k: int = 4) -> set:
    """Extracts all unique k-mers from a sequence using a set."""
    kmers = set()
    seq = sequence.upper()
    for i in range(len(seq) - k + 1):
        kmers.add(seq[i:i+k])
    return kmers

def find_restriction_sites(sequence: str) -> list:
    """Scans sequence for standard restriction enzyme sites."""
    seq = sequence.upper()
    matches = []
    for enzyme, motif in RESTRICTION_SITES.items():
        start = 0
        while True:
            idx = seq.find(motif, start)
            if idx == -1:
                break
            matches.append((enzyme, motif, idx + 1))  # 1-based coordinate tuple
            start = idx + 1
    return matches

def main():
    parser = argparse.ArgumentParser(description="Chapter 2: Genomic data structures demo.")
    parser.add_argument("-s", "--seq", help="Input DNA sequence")
    args = parser.parse_args()

    seq = args.seq.strip() if args.seq else "ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGAATTCGTTACTGCCGGATCCCTGTGGTAA"
    print(f"Analyzing Sequence ({len(seq)} bp): {seq}\n")

    # 1. Translation via Dictionary
    protein = translate_dna(seq)
    print(f"[1] Translation: {protein}")

    # 2. Unique k-mers via Set
    kmers = find_unique_kmers(seq, k=4)
    print(f"[2] Unique 4-mers (Set count: {len(kmers)}): {sorted(list(kmers))[:8]}...")

    # 3. Restriction sites as Tuples in a List
    sites = find_restriction_sites(seq)
    print(f"[3] Restriction Enzyme Matches (List of Tuples):")
    if sites:
        for enzyme, motif, pos in sites:
            print(f"    - {enzyme} ({motif}) at position {pos}")
    else:
        print("    - No restriction sites found.")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 3: Python for Next Gen Biologist
Script: ch03_functions_orfs.py
Description: Defining modular functions, scanning all 6 reading frames for Open Reading Frames (ORFs),
             passing arguments, returning structured data, and recursion.
Usage:
    python ch03_functions_orfs.py --seq ATGAAACCCGGGTTTTAA
"""

import argparse

def reverse_complement(sequence: str) -> str:
    comp = str.maketrans("ACGTacgt", "TGCAtgca")
    return sequence.translate(comp)[::-1]

def find_orfs_in_frame(sequence: str, frame: int, min_len_aa: int = 15) -> list:
    """Find open reading frames in a specific reading frame (0, 1, or 2)."""
    seq = sequence.upper()
    start_codons = {"ATG"}
    stop_codons = {"TAA", "TAG", "TGA"}
    orfs = []

    start_idx = None
    for i in range(frame, len(seq) - 2, 3):
        codon = seq[i:i+3]
        if codon in start_codons and start_idx is None:
            start_idx = i
        elif codon in stop_codons and start_idx is not None:
            orf_dna = seq[start_idx:i+3]
            length_aa = len(orf_dna) // 3
            if length_aa >= min_len_aa:
                orfs.append({
                    'start': start_idx + 1,  # 1-based coordinate
                    'end': i + 3,
                    'length_nt': len(orf_dna),
                    'length_aa': length_aa,
                    'dna': orf_dna
                })
            start_idx = None
    return orfs

def scan_all_six_frames(sequence: str, min_len_aa: int = 10) -> dict:
    """Scan both forward (+) and reverse complement (-) strands across 3 frames each."""
    forward_seq = sequence.upper()
    revcomp_seq = reverse_complement(forward_seq)
    
    results = {'forward': {}, 'reverse': {}}
    for f in range(3):
        results['forward'][f"+F{f+1}"] = find_orfs_in_frame(forward_seq, f, min_len_aa)
        results['reverse'][f"-F{f+1}"] = find_orfs_in_frame(revcomp_seq, f, min_len_aa)
    return results

def recursive_nucleotide_distance(seq1: str, seq2: str) -> int:
    """Recursive Hamming edit distance calculation between two equal-length sequences."""
    if len(seq1) != len(seq2):
        raise ValueError("Sequences must be of equal length for Hamming distance")
    if len(seq1) == 0:
        return 0
    head_dist = 1 if seq1[0] != seq2[0] else 0
    return head_dist + recursive_nucleotide_distance(seq1[1:], seq2[1:])

def main():
    parser = argparse.ArgumentParser(description="Chapter 3: Functions and 6-frame ORF scanner.")
    parser.add_argument("-s", "--seq", help="Input DNA sequence")
    parser.add_argument("--min-aa", type=int, default=10, help="Minimum ORF length in amino acids")
    args = parser.parse_args()

    seq = args.seq.strip() if args.seq else (
        "GGGATGGCCACCGATGGCATGCCCAAGCTGTACGACTACGTGTACGAGCGCCTCGACGAGACCGAGGAG"
        "TTCGCCCGCTTCGAGGCCGCCAACATGCTCCGCTACCGCGCCCGCACCCGCTAGTTTTTT"
    )

    print(f"Scanning sequence ({len(seq)} bp) for ORFs (Min AA: {args.min_aa})...\n")
    six_frames = scan_all_six_frames(seq, min_len_aa=args.min_aa)
    
    for strand, frames in six_frames.items():
        for frame_name, orfs in frames.items():
            if orfs:
                print(f"Strand: {strand.upper()} | Frame: {frame_name} ({len(orfs)} ORFs)")
                for orf in orfs:
                    print(f"  -> {orf['start']}-{orf['end']} ({orf['length_nt']} bp, {orf['length_aa']} AA): {orf['dna'][:30]}...")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 4: Python for Next Gen Biologist
Script: ch04_fastq_stream_parser.py
Description: Streaming FASTQ parser, ASCII Phred+33 score conversion, average read quality calculation,
             and quality filtering (Q30 threshold).
Usage:
    python ch04_fastq_stream_parser.py --input reads.fastq --min-q 30 --min-len 50
"""

import argparse
import sys
import gzip

def parse_fastq(filepath):
    """Generator yielding (header, seq, comment, qual_str) tuples from plain or gzipped FASTQ."""
    opener = gzip.open if filepath.endswith('.gz') else open
    with opener(filepath, 'rt') as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            comment = f.readline().strip()
            qual = f.readline().strip()
            yield header, seq, comment, qual

def phred_scores(qual_str: str, phred_offset: int = 33) -> list:
    """Converts ASCII quality characters to Phred quality integers."""
    return [ord(c) - phred_offset for c in qual_str]

def mean_quality(scores: list) -> float:
    """Calculates arithmetic mean quality score."""
    return sum(scores) / len(scores) if scores else 0.0

def filter_fastq(input_path, output_path, min_q=30, min_len=50):
    total = 0
    passed = 0
    
    out_file = open(output_path, 'w') if output_path else None
    
    for header, seq, comment, qual in parse_fastq(input_path):
        total += 1
        scores = phred_scores(qual)
        avg_q = mean_quality(scores)
        
        if avg_q >= min_q and len(seq) >= min_len:
            passed += 1
            if out_file:
                out_file.write(f"{header}\n{seq}\n{comment}\n{qual}\n")
    
    if out_file:
        out_file.close()
        
    return total, passed

def main():
    parser = argparse.ArgumentParser(description="Chapter 4: FASTQ Stream Parser & Quality Filter.")
    parser.add_argument("-i", "--input", help="Path to input FASTQ file (.fastq or .fastq.gz)")
    parser.add_argument("-o", "--output", help="Optional output path for filtered FASTQ")
    parser.add_argument("-q", "--min-q", type=float, default=30.0, help="Minimum average Phred quality score (Default: 30)")
    parser.add_argument("-l", "--min-len", type=int, default=50, help="Minimum read length (Default: 50)")
    args = parser.parse_args()

    if not args.input:
        # Create a synthetic mini FASTQ file to demo execution
        demo_file = "/tmp/demo_reads.fastq"
        with open(demo_file, "w") as f:
            f.write("@READ_001_HIGH_QUAL\nACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT\n+\nIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\n")
            f.write("@READ_002_LOW_QUAL\nNNNNACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT\n+\n######IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\n")
            f.write("@READ_003_SHORT\nACGTACGT\n+\nIIIIIIII\n")
        print(f"[!] No input specified. Created demo FASTQ with 3 reads at {demo_file}\n")
        args.input = demo_file

    total, passed = filter_fastq(args.input, args.output, min_q=args.min_q, min_len=args.min_len)
    pass_pct = (passed / total * 100) if total > 0 else 0
    print(f"=== FASTQ Quality Filtering Summary ===")
    print(f"Total Reads Analyzed: {total}")
    print(f"Passed Reads (Q>={args.min_q}, Len>={args.min_len}): {passed} ({pass_pct:.2f}%)")
    print(f"Filtered Out Reads  : {total - passed}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 5: Python for Next Gen Biologist
Script: ch05_fasta_csv_io.py
Description: Robust FASTA reader/writer handling multi-line sequences without high memory usage;
             CSV/TSV expression matrix parser calculating per-gene mean and variance.
Usage:
    python ch05_fasta_csv_io.py --fasta sample.fasta
    python ch05_fasta_csv_io.py --matrix counts.tsv
"""

import argparse
import csv
import io
import math

def read_fasta(filepath):
    """Generator yielding (header, sequence) records from multi-line FASTA files."""
    header = None
    seq_chunks = []
    with open(filepath, 'r') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            if line.startswith('>'):
                if header is not None:
                    yield header, "".join(seq_chunks)
                header = line[1:]
                seq_chunks = []
            else:
                seq_chunks.append(line)
        if header is not None:
            yield header, "".join(seq_chunks)

def write_fasta(records, output_filepath, line_width=60):
    """Writes list of (header, sequence) tuples to FASTA with wrapped lines."""
    with open(output_filepath, 'w') as f:
        for header, seq in records:
            f.write(f">{header}\n")
            for i in range(0, len(seq), line_width):
                f.write(f"{seq[i:i+line_width]}\n")

def parse_expression_tsv(filepath):
    """Parses a gene expression count matrix (GeneID \t S1 \t S2 ...) and computes mean/variance."""
    results = []
    with open(filepath, 'r') as f:
        reader = csv.reader(f, delimiter='\t')
        headers = next(reader)
        samples = headers[1:]
        
        for row in reader:
            if not row:
                continue
            gene_id = row[0]
            counts = [float(x) for x in row[1:]]
            n = len(counts)
            mean = sum(counts) / n if n > 0 else 0
            variance = sum((x - mean) ** 2 for x in counts) / (n - 1) if n > 1 else 0
            results.append({
                'gene_id': gene_id,
                'mean': round(mean, 2),
                'variance': round(variance, 2),
                'std_dev': round(math.sqrt(variance), 2)
            })
    return samples, results

def main():
    parser = argparse.ArgumentParser(description="Chapter 5: FASTA & CSV/TSV File I/O.")
    parser.add_argument("--fasta", help="Path to input FASTA")
    parser.add_argument("--matrix", help="Path to input TSV expression matrix")
    args = parser.parse_args()

    # Demo FASTA
    demo_fasta = "/tmp/demo_genes.fasta"
    with open(demo_fasta, "w") as f:
        f.write(">geneA Homo_sapiens chr1:100-300\nATGCGATCGATCGATCGATCGATCGATC\nGATCGATCGATCGATCGATCGATCGATCG\n")
        f.write(">geneB Mus_musculus chr2:500-750\nATGGCCAATTGGCCAAATTTGGCCCAAA\nTTTAAACCCGGGTTTAAACCCGGGTTTA\n")

    print("=== Demo 1: Reading Multi-line FASTA ===")
    for h, s in read_fasta(demo_fasta):
        print(f"Header: {h} | Length: {len(s)} bp | Preview: {s[:25]}...")

    # Demo TSV Matrix
    demo_tsv = "/tmp/demo_counts.tsv"
    with open(demo_tsv, "w") as f:
        f.write("GeneID\tControl_1\tControl_2\tTreated_1\tTreated_2\n")
        f.write("ENSG000001\t120\t115\t450\t480\n")
        f.write("ENSG000002\t2300\t2450\t2280\t2310\n")
        f.write("ENSG000003\t10\t12\t95\t105\n")

    print("\n=== Demo 2: Parsing TSV Count Matrix ===")
    samples, stats = parse_expression_tsv(demo_tsv)
    print(f"Samples ({len(samples)}): {', '.join(samples)}")
    for g in stats:
        print(f"Gene: {g['gene_id']:12} | Mean Count: {g['mean']:6} | StdDev: {g['std_dev']:6}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 6: Python for Next Gen Biologist
Script: ch06_dna_oop.py
Description: Object-Oriented Programming (OOP) in Python for Biologists:
             DnaSequence, ProteinSequence, and GenomicVariant classes with inheritance,
             transcription, translation, and mutation application.
Usage:
    python ch06_dna_oop.py
"""

import math

class BiologicalSequence:
    """Abstract base class representing an arbitrary biological sequence."""
    
    def __init__(self, seq_id: str, sequence: str, organism: str = "Unknown"):
        self.seq_id = str(seq_id)
        self.sequence = sequence.upper().strip()
        self.organism = organism
        
    def __len__(self):
        return len(self.sequence)
        
    def __repr__(self):
        return f"<{self.__class__.__name__} id='{self.seq_id}' len={len(self)} organism='{self.organism}'>"
        
    def __getitem__(self, item):
        return self.sequence[item]


class DnaSequence(BiologicalSequence):
    """Represents a double-stranded genomic DNA sequence."""
    
    VALID_BASES = set("ACGTN")
    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',
    }

    def __init__(self, seq_id: str, sequence: str, organism: str = "Homo sapiens"):
        super().__init__(seq_id, sequence, organism)
        if not set(self.sequence).issubset(self.VALID_BASES):
            raise ValueError(f"Sequence {seq_id} contains invalid DNA bases!")

    def gc_content(self) -> float:
        """Returns the GC percentage of the sequence."""
        if len(self) == 0:
            return 0.0
        gc = self.sequence.count('G') + self.sequence.count('C')
        return round((gc / len(self)) * 100, 2)

    def reverse_complement(self) -> "DnaSequence":
        """Returns a new DnaSequence instance with the reverse complement."""
        trans = str.maketrans("ACGTN", "TGCAN")
        rev_seq = self.sequence.translate(trans)[::-1]
        return DnaSequence(f"{self.seq_id}_revcomp", rev_seq, self.organism)

    def transcribe(self) -> str:
        """Returns the single-stranded mRNA transcript string (T -> U)."""
        return self.sequence.replace("T", "U")

    def translate(self) -> "ProteinSequence":
        """Translates coding DNA into a ProteinSequence object."""
        peptide = []
        for i in range(0, len(self) - 2, 3):
            codon = self.sequence[i:i+3]
            aa = self.CODON_TABLE.get(codon, 'X')
            if aa == '*':
                break
            peptide.append(aa)
        return ProteinSequence(f"{self.seq_id}_protein", "".join(peptide), self.organism)

    def apply_variant(self, variant: "GenomicVariant") -> "DnaSequence":
        """Applies a single-nucleotide variant (SNV) to create a mutated sequence."""
        pos = variant.pos - 1  # Convert 1-based to 0-based
        if pos < 0 or pos >= len(self):
            raise IndexError("Variant position out of bounds.")
        if self.sequence[pos] != variant.ref:
            raise ValueError(f"Reference mismatch at pos {variant.pos}: expected {variant.ref}, found {self.sequence[pos]}")
            
        mutated = self.sequence[:pos] + variant.alt + self.sequence[pos+1:]
        return DnaSequence(f"{self.seq_id}_{variant.ref}{variant.pos}{variant.alt}", mutated, self.organism)


class ProteinSequence(BiologicalSequence):
    """Represents an amino acid peptide sequence."""
    
    # Molecular weights of amino acids in Daltons (Da)
    AA_WEIGHTS = {
        'A': 71.08, 'R': 156.20, 'N': 114.11, 'D': 115.09, 'C': 103.14,
        'E': 129.12, 'Q': 128.13, 'G': 57.05,  'H': 137.14, 'I': 113.17,
        'L': 113.17, 'K': 128.18, 'M': 131.21, 'F': 147.18, 'P': 97.12,
        'S': 87.08,  'T': 101.11, 'W': 186.21, 'Y': 163.18, 'V': 99.13
    }

    def molecular_weight_kda(self) -> float:
        """Calculates molecular weight in kiloDaltons including terminal water molecule."""
        weight = 18.015  # H2O terminal
        for aa in self.sequence:
            weight += self.AA_WEIGHTS.get(aa, 110.0)
        return round(weight / 1000.0, 2)


class GenomicVariant:
    """Represents a genomic single-nucleotide variant (SNV)."""
    
    def __init__(self, chrom: str, pos: int, ref: str, alt: str):
        self.chrom = str(chrom)
        self.pos = int(pos)
        self.ref = str(ref).upper()
        self.alt = str(alt).upper()
        
    def __repr__(self):
        return f"<GenomicVariant {self.chrom}:{self.pos} {self.ref}>{self.alt}>"


def main():
    print("=== Chapter 6: Object-Oriented Biological Sequences ===")
    
    # 1. Instantiate Wild-Type DNA Sequence (Human Beta-Globin Exon 1 fragment)
    wt_dna = DnaSequence(
        seq_id="HBB_WT",
        sequence="ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGG",
        organism="Homo sapiens"
    )
    print(f"Wild-Type DNA     : {wt_dna}")
    print(f"GC Content        : {wt_dna.gc_content()}%")
    
    wt_protein = wt_dna.translate()
    print(f"Translated Protein: {wt_protein} ({wt_protein.sequence})")
    print(f"Molecular Weight  : {wt_protein.molecular_weight_kda()} kDa")

    # 2. Apply Sickle Cell Mutation: HBB c.20A>T (p.Glu7Val / HbS)
    # Codon 6: GAG -> GTG at nucleotide position 20
    scd_var = GenomicVariant(chrom="chr11", pos=20, ref="A", alt="T")
    print(f"\nApplying Variant  : {scd_var}")
    
    mut_dna = wt_dna.apply_variant(scd_var)
    mut_protein = mut_dna.translate()
    
    print(f"Mutated DNA       : {mut_dna}")
    print(f"Mutated Protein   : {mut_protein} ({mut_protein.sequence})")
    print(f"Amino Acid Delta  : Pos 7 changed from '{wt_protein[6]}' to '{mut_protein[6]}'")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 7: Python for Next Gen Biologist
Script: ch07_modules_packages.py
Description: Designing modular bioinformatics code, structuring packages,
             namespace separation, and reusable genomics utility libraries.
Usage:
    python ch07_modules_packages.py --seq ATGCGATCGATC
"""

import argparse
import sys

class GenomicsToolkit:
    """Namespace container demonstrating structured module design."""
    
    @staticmethod
    def clean_sequence(raw_seq: str) -> str:
        """Strips whitespace, numbers, and converts to uppercase."""
        return "".join(c for c in raw_seq if c.isalpha()).upper()

    @staticmethod
    def expand_iupac(motif: str) -> list:
        """Expands degenerate IUPAC nucleotide codes into all explicit combinations."""
        iupac_map = {
            'A': ['A'], 'C': ['C'], 'G': ['G'], 'T': ['T'],
            'R': ['A', 'G'], 'Y': ['C', 'T'], 'S': ['G', 'C'], 'W': ['A', 'T'],
            'K': ['G', 'T'], 'M': ['A', 'C'], 'B': ['C', 'G', 'T'],
            'D': ['A', 'G', 'T'], 'H': ['A', 'C', 'T'], 'V': ['A', 'C', 'G'],
            'N': ['A', 'C', 'G', 'T']
        }
        sequences = [""]
        for base in motif.upper():
            candidates = iupac_map.get(base, [base])
            sequences = [prefix + c for prefix in sequences for c in candidates]
        return sequences

    @staticmethod
    def compute_dinucleotide_frequencies(sequence: str) -> dict:
        """Calculates observed frequencies of all 16 dinucleotide combinations (e.g. CpG)."""
        seq = sequence.upper()
        dinucs = {}
        for b1 in "ACGT":
            for b2 in "ACGT":
                dinucs[b1 + b2] = 0
                
        for i in range(len(seq) - 1):
            pair = seq[i:i+2]
            if pair in dinucs:
                dinucs[pair] += 1
        return dinucs

def main():
    parser = argparse.ArgumentParser(description="Chapter 7: Genomics Toolkit Module & Packaging Demo.")
    parser.add_argument("-s", "--seq", help="Input DNA sequence")
    parser.add_argument("-m", "--motif", default="RGYW", help="IUPAC degenerate motif to expand (default: RGYW)")
    args = parser.parse_args()

    seq = args.seq.strip() if args.seq else "ATG CGG CCG 123 ATCGATCGATCG TACG ATCG"
    cleaned = GenomicsToolkit.clean_sequence(seq)
    print(f"Original Raw Sequence : {seq}")
    print(f"Cleaned Sequence ({len(cleaned)} bp) : {cleaned}\n")

    expansions = GenomicsToolkit.expand_iupac(args.motif)
    print(f"Expanded IUPAC Motif '{args.motif}' ({len(expansions)} combinations):")
    print(f"  -> {', '.join(expansions)}\n")

    dinucs = GenomicsToolkit.compute_dinucleotide_frequencies(cleaned)
    cpg_count = dinucs.get("CG", 0)
    print(f"Dinucleotide Profile: CpG count = {cpg_count}")
    top_dinucs = sorted(dinucs.items(), key=lambda x: x[1], reverse=True)[:5]
    print(f"Top 5 Dinucleotides : {top_dinucs}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 8: Python for Next Gen Biologist
Script: ch08_regex_data_cleaning.py
Description: Regular expressions for biological sequence patterns:
             Zinc finger C2H2 motifs, bacterial Pribnow box, restriction cut sites,
             and clinical metadata cleaning.
Usage:
    python ch08_regex_data_cleaning.py
"""

import re
import csv
import io

def find_c2h2_zinc_fingers(protein_seq: str) -> list:
    """
    Finds classic C2H2 zinc finger motifs: C-x(2,4)-C-x(12)-H-x(3,5)-H
    """
    pattern = re.compile(r'C.{2,4}C.{12}H.{3,5}H')
    matches = []
    for match in pattern.finditer(protein_seq):
        matches.append({
            'start': match.start() + 1,  # 1-based coordinate
            'end': match.end(),
            'motif': match.group()
        })
    return matches

def find_promoter_pribnow_box(dna_seq: str) -> list:
    """
    Scans for bacterial -10 Pribnow box consensus (TATAAT) with up to 1 mismatch.
    Regex uses alternation for core conservation.
    """
    # Pattern allowing variation in central positions
    pattern = re.compile(r'TA[TA][AA]T', re.IGNORECASE)
    matches = []
    for m in pattern.finditer(dna_seq):
        matches.append((m.start() + 1, m.end(), m.group().upper()))
    return matches

def clean_clinical_metadata(raw_csv_data: str) -> list:
    """
    Cleans messy patient/sample metadata: standardizes sample IDs, cleans ages,
    and formats diagnosis stages.
    """
    reader = csv.DictReader(io.StringIO(raw_csv_data.strip()))
    cleaned = []
    for row in reader:
        # Standardize Sample ID: e.g. " patient_001 " -> "PT-001"
        sid = re.sub(r'[^0-9]', '', row['sample_id'])
        clean_sid = f"PT-{int(sid):03d}" if sid else "UNKNOWN"
        
        # Clean Age: "54 yrs" / "54yo" -> 54
        age_match = re.search(r'\d+', row['age'])
        age = int(age_match.group()) if age_match else None
        
        # Clean Stage: "stage IIIa" -> "Stage IIIA"
        stage = row['stage'].strip().upper()
        stage = re.sub(r'\s+', ' ', stage)
        
        cleaned.append({
            'sample_id': clean_sid,
            'age': age,
            'gender': row['gender'].strip().upper()[:1],
            'stage': stage
        })
    return cleaned

def main():
    print("=== Chapter 8: Genomic Regular Expressions & Data Cleaning ===")

    # Demo 1: C2H2 Zinc Finger Search
    protein = "MKGEELFTGVVPILVELDGDVNGHKFSVSGEGEGDATYGKLTLKFICTYCKTFVLDDNLQAHVTHMGPVI"
    zf_matches = find_c2h2_zinc_fingers(protein)
    print(f"\n[1] Scanning Protein ({len(protein)} AA) for C2H2 Zinc Fingers:")
    if zf_matches:
        for m in zf_matches:
            print(f"    -> Found Motif at {m['start']}-{m['end']}: {m['motif']}")
    else:
        print("    -> No standard C2H2 motif detected.")

    # Demo 2: Bacterial Pribnow Box (-10 Promoter)
    promoter_dna = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAA"
    pribnow = find_promoter_pribnow_box(promoter_dna)
    print(f"\n[2] Scanning DNA for Pribnow Box (-10 Promoter):")
    for start, end, seq in pribnow:
        print(f"    -> Consensus Match at bp {start}-{end}: {seq}")

    # Demo 3: Cleaning Clinical Metadata Table
    dirty_csv = """sample_id,age,gender,stage
patient_1, 58 yrs ,Female,stage IIa
pt-02, 63yo ,M,STAGE 3B
sample 003,45 years,F,stage i
"""
    cleaned = clean_clinical_metadata(dirty_csv)
    print(f"\n[3] Cleaned Clinical Metadata Table:")
    for row in cleaned:
        print(f"    Sample: {row['sample_id']} | Age: {row['age']} | Sex: {row['gender']} | Stage: {row['stage']}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 9: Python for Next Gen Biologist
Script: ch09_numpy_pandas_viz.py
Description: Scientific computing with NumPy, Pandas, and Matplotlib:
             Position Weight Matrix (PWM) log-odds scoring, RNA-seq counts normalization,
             fold-change filtering, and Volcano plot generation.
Usage:
    python ch09_numpy_pandas_viz.py [--output volcano.png]
"""

import argparse
import math
import numpy as np
import pandas as pd

def build_pwm_log_odds(aligned_motifs: list, background: float = 0.25) -> np.ndarray:
    """
    Builds a Position Weight Matrix (PWM) in log2-odds space from aligned motif strings.
    Matrix shape: (4, motif_length) corresponding to rows A, C, G, T.
    """
    length = len(aligned_motifs[0])
    base_indices = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
    counts = np.ones((4, length))  # Laplace pseudocount = 1
    
    for seq in aligned_motifs:
        for col, base in enumerate(seq.upper()):
            if base in base_indices:
                counts[base_indices[base], col] += 1
                
    frequencies = counts / counts.sum(axis=0)
    log_odds = np.log2(frequencies / background)
    return log_odds

def score_sequence_with_pwm(sequence: str, pwm: np.ndarray) -> list:
    """Scans sequence with sliding window and returns list of (pos, score)."""
    base_indices = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
    motif_len = pwm.shape[1]
    seq = sequence.upper()
    scores = []
    
    for i in range(len(seq) - motif_len + 1):
        window = seq[i:i+motif_len]
        score = 0.0
        valid = True
        for col, b in enumerate(window):
            if b in base_indices:
                score += pwm[base_indices[b], col]
            else:
                valid = False
                break
        if valid:
            scores.append((i + 1, round(score, 3)))
    return scores

def simulate_rnaseq_dataframe(num_genes: int = 100) -> pd.DataFrame:
    """Generates synthetic RNA-seq counts DataFrame for demonstration."""
    np.random.seed(42)
    genes = [f"GENE_{i:04d}" for i in range(1, num_genes + 1)]
    
    # Simulate log2 fold changes and p-values
    log2fc = np.random.normal(loc=0.0, scale=1.5, size=num_genes)
    # Give a few genes strong upregulation/downregulation
    log2fc[:5] += 3.5
    log2fc[5:10] -= 3.5
    
    pvalues = np.random.beta(a=0.5, b=2.0, size=num_genes)
    pvalues[:10] = pvalues[:10] * 1e-4  # strong significance
    
    neg_log10_p = -np.log10(np.clip(pvalues, 1e-15, 1.0))
    
    df = pd.DataFrame({
        'GeneID': genes,
        'log2FoldChange': np.round(log2fc, 3),
        'pvalue': pvalues,
        'negLog10P': np.round(neg_log10_p, 3)
    })
    
    df['Status'] = 'Not Significant'
    df.loc[(df['log2FoldChange'] > 1.5) & (df['pvalue'] < 0.05), 'Status'] = 'Upregulated'
    df.loc[(df['log2FoldChange'] < -1.5) & (df['pvalue'] < 0.05), 'Status'] = 'Downregulated'
    return df

def main():
    parser = argparse.ArgumentParser(description="Chapter 9: NumPy & Pandas Bioinformatics Demo.")
    parser.add_argument("-o", "--output", help="Optional path to save volcano plot")
    args = parser.parse_args()

    print("=== Chapter 9: Scientific Computing with NumPy & Pandas ===")

    # 1. NumPy PWM
    motifs = [
        "TATAAA", "TATATA", "TATAAA", "TATAAT", "CATAAA", "TATAAA"
    ]
    pwm = build_pwm_log_odds(motifs)
    print(f"\n[1] Generated PWM Log-Odds Matrix ({pwm.shape[0]} bases x {pwm.shape[1]} pos):")
    print("    A pos 1-6:", np.round(pwm[0], 2))
    print("    T pos 1-6:", np.round(pwm[3], 2))

    target = "CGTATAAACCGGATATAAATCG"
    scores = score_sequence_with_pwm(target, pwm)
    print(f"\nScanning Target DNA: {target}")
    for pos, s in sorted(scores, key=lambda x: x[1], reverse=True)[:3]:
        print(f"    -> Pos {pos}: Score = {s}")

    # 2. Pandas Differential Expression Table
    df = simulate_rnaseq_dataframe(100)
    print(f"\n[2] Differential Expression Summary (Pandas DataFrame):")
    counts = df['Status'].value_counts()
    for status, count in counts.items():
        print(f"    - {status:15}: {count} genes")

    top_up = df[df['Status'] == 'Upregulated'].sort_values('log2FoldChange', ascending=False).head(3)
    print("\nTop 3 Upregulated Genes:")
    print(top_up[['GeneID', 'log2FoldChange', 'negLog10P']])

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 10: Python for Next Gen Biologist
Script: ch10_biopython_toolkit.py
Description: Biopython core workflows: Seq and SeqRecord objects, FASTA parsing,
             simulated NCBI Entrez E-Utilities fetching, and BLAST output parsing.
Usage:
    python ch10_biopython_toolkit.py
"""

import sys
import io

def demo_biopython_seq():
    """Demonstrates Seq object features with native fallback if biopython is not installed."""
    try:
        from Bio.Seq import Seq
        from Bio.SeqRecord import SeqRecord
        
        my_seq = Seq("ATGCGATCGATCGATCGATAGCTAGCTA")
        record = SeqRecord(my_seq, id="NC_000913.3", description="E. coli K-12 complete genome fragment")
        
        print("=== Biopython Seq & SeqRecord Demo ===")
        print(f"Record ID    : {record.id}")
        print(f"Description  : {record.description}")
        print(f"Sequence (bp): {len(record.seq)}")
        print(f"Transcription: {record.seq.transcribe()[:20]}...")
        print(f"Translation  : {record.seq.translate()}")
        return True
    except ImportError:
        print("[NOTICE] 'biopython' is not installed in the current environment.")
        print("Run: pip install biopython")
        print("\nDemonstrating pure-Python architectural equivalent:")
        raw_dna = "ATGCGATCGATCGATCGATAGCTAGCTA"
        mrna = raw_dna.replace("T", "U")
        print(f"Sequence : {raw_dna}")
        print(f"mRNA     : {mrna}")
        return False

def demo_blast_xml_parsing():
    """Demonstrates parsing NCBI BLAST XML output."""
    xml_data = """<?xml version="1.0"?>
    <BlastOutput>
      <BlastOutput_iterations>
        <Iteration>
          <Iteration_hits>
            <Hit>
              <Hit_num>1</Hit_num>
              <Hit_id>ref|NM_000518.5|</Hit_id>
              <Hit_def>Homo sapiens hemoglobin subunit beta (HBB), mRNA</Hit_def>
              <Hit_len>628</Hit_len>
              <Hit_hsps>
                <Hsp>
                  <Hsp_bit-score>1160.2</Hsp_bit-score>
                  <Hsp_evalue>0.0</Hsp_evalue>
                  <Hsp_identity>628</Hsp_identity>
                </Hsp>
              </Hit_hsps>
            </Hit>
          </Iteration_hits>
        </Iteration>
      </BlastOutput_iterations>
    </BlastOutput>
    """
    print("\n=== BLAST Output Parser Demo ===")
    import xml.etree.ElementTree as ET
    root = ET.fromstring(xml_data.strip())
    hits = root.findall(".//Hit")
    for hit in hits:
        hit_id = hit.find("Hit_id").text
        hit_def = hit.find("Hit_def").text
        eval_elem = hit.find(".//Hsp_evalue")
        evalue = eval_elem.text if eval_elem is not None else "N/A"
        print(f"Top Hit: {hit_id}")
        print(f"Def    : {hit_def}")
        print(f"E-value: {evalue}")

def main():
    has_bio = demo_biopython_seq()
    demo_blast_xml_parsing()

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 11: Python for Next Gen Biologist
Script: ch11_rnaseq_deseq_analysis.py
Description: Analyzing biological datasets: RNA-Seq differential expression testing,
             two-sample Welch's t-test, Benjamini-Hochberg False Discovery Rate (FDR) adjustment,
             and PPI / co-expression network graph creation.
Usage:
    python ch11_rnaseq_deseq_analysis.py
"""

import math
import numpy as np
import pandas as pd

def welch_t_test(group1: np.ndarray, group2: np.ndarray) -> tuple:
    """Computes Welch's t-statistic and approximate two-tailed p-value."""
    n1, n2 = len(group1), len(group2)
    m1, m2 = np.mean(group1), np.mean(group2)
    v1, v2 = np.var(group1, ddof=1), np.var(group2, ddof=1)
    
    denom = math.sqrt((v1 / n1) + (v2 / n2))
    if denom == 0:
        return 0.0, 1.0
        
    t_stat = (m1 - m2) / denom
    # Welch-Satterthwaite degrees of freedom
    df_num = ((v1 / n1) + (v2 / n2)) ** 2
    df_denom = ((v1 / n1) ** 2 / (n1 - 1)) + ((v2 / n2) ** 2 / (n2 - 1))
    df = df_num / df_denom if df_denom > 0 else 1.0
    
    # Rough normal approximation for p-value
    p_val = 2 * (1 - 0.5 * (1 + math.erf(abs(t_stat) / math.sqrt(2))))
    return t_stat, max(p_val, 1e-15)

def benjamini_hochberg(p_values: list) -> list:
    """Performs Benjamini-Hochberg FDR correction on a list of p-values."""
    n = len(p_values)
    indexed_p = sorted(enumerate(p_values), key=lambda x: x[1])
    adjusted = [0.0] * n
    
    cum_min = 1.0
    for rank_minus_1 in range(n - 1, -1, -1):
        orig_idx, p = indexed_p[rank_minus_1]
        rank = rank_minus_1 + 1
        adj_p = (p * n) / rank
        cum_min = min(cum_min, adj_p)
        adjusted[orig_idx] = min(cum_min, 1.0)
        
    return adjusted

def main():
    print("=== Chapter 11: RNA-Seq Differential Expression & Statistical Testing ===")
    np.random.seed(101)
    
    num_genes = 20
    ctrl_samples = 3
    treat_samples = 3
    
    # Generate mock log2-normalized expression counts
    ctrl_data = np.random.normal(loc=8.0, scale=0.5, size=(num_genes, ctrl_samples))
    treat_data = np.random.normal(loc=8.0, scale=0.5, size=(num_genes, treat_samples))
    
    # Shift top 4 genes to simulate true biological perturbation
    treat_data[:2, :] += 2.5   # Upregulated
    treat_data[2:4, :] -= 2.5  # Downregulated
    
    results = []
    p_vals = []
    
    for g in range(num_genes):
        gene_name = f"GENE_{g+1:03d}"
        c = ctrl_data[g, :]
        t = treat_data[g, :]
        log2fc = np.mean(t) - np.mean(c)
        t_stat, p = welch_t_test(t, c)
        p_vals.append(p)
        results.append({
            'Gene': gene_name,
            'Control_Mean': round(np.mean(c), 2),
            'Treated_Mean': round(np.mean(t), 2),
            'Log2FC': round(log2fc, 2),
            'pValue': p
        })
        
    fdr_vals = benjamini_hochberg(p_vals)
    for i, r in enumerate(results):
        r['FDR_adj_p'] = round(fdr_vals[i], 5)
        r['Significant'] = (abs(r['Log2FC']) >= 1.5) and (r['FDR_adj_p'] < 0.05)
        
    df = pd.DataFrame(results)
    sig_genes = df[df['Significant'] == True]
    
    print(f"\nDifferential Expression Table (Top 6 Genes):")
    print(df[['Gene', 'Log2FC', 'pValue', 'FDR_adj_p', 'Significant']].head(6))
    print(f"\nTotal Genes Assayed : {num_genes}")
    print(f"Significant (FDR<0.05, |Log2FC|>=1.5): {len(sig_genes)}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 12: Python for Next Gen Biologist
Script: ch12_ml_cancer_classifier.py
Description: Machine Learning for Biology:
             Feature scaling, Principal Component Analysis (PCA),
             and Random Forest Classifier for cancer subtype classification.
Usage:
    python ch12_ml_cancer_classifier.py
"""

import numpy as np

def run_ml_pipeline():
    try:
        from sklearn.ensemble import RandomForestClassifier
        from sklearn.decomposition import PCA
        from sklearn.model_selection import StratifiedKFold, cross_val_score
        from sklearn.preprocessing import StandardScaler
        
        print("=== Chapter 12: Machine Learning for Biology (Scikit-Learn) ===")
        np.random.seed(42)
        
        # Simulate 60 patient expression profiles (30 Luminal A, 30 Basal-like) across 200 biomarker genes
        n_samples = 60
        n_features = 200
        
        X = np.random.normal(loc=5.0, scale=1.0, size=(n_samples, n_features))
        y = np.array([0] * 30 + [1] * 30)  # 0 = Luminal A, 1 = Basal-like
        
        # Introduce distinct biomarker signature in first 10 genes
        X[:30, :10] += 2.0   # High in Luminal A
        X[30:, 10:20] += 2.0 # High in Basal-like
        
        # 1. Feature Standardization
        scaler = StandardScaler()
        X_scaled = scaler.fit_transform(X)
        
        # 2. PCA Dimensionality Reduction
        pca = PCA(n_components=2)
        X_pca = pca.fit_transform(X_scaled)
        var_exp = pca.explained_variance_ratio_ * 100
        print(f"\n[1] PCA Decomposition:")
        print(f"    - PC1 Variance Explained: {var_exp[0]:.2f}%")
        print(f"    - PC2 Variance Explained: {var_exp[1]:.2f}%")
        
        # 3. Random Forest Classification with 5-Fold Stratified Cross-Validation
        rf = RandomForestClassifier(n_estimators=50, random_state=42)
        cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
        scores = cross_val_score(rf, X_scaled, y, cv=cv, scoring='accuracy')
        
        print(f"\n[2] Random Forest 5-Fold Cross-Validation Accuracy:")
        print(f"    - Mean Accuracy : {np.mean(scores)*100:.2f}% (+/- {np.std(scores)*100:.2f}%)")
        print(f"    - Fold Scores   : {[round(s, 3) for s in scores]}")
        
        # 4. Feature Importance for Biomarker Discovery
        rf.fit(X_scaled, y)
        importances = rf.feature_importances_
        top_genes = np.argsort(importances)[::-1][:5]
        print(f"\n[3] Top 5 Predictive Biomarker Genes:")
        for rank, gene_idx in enumerate(top_genes, 1):
            print(f"    Rank {rank}: Gene_{gene_idx:03d} (Importance: {importances[gene_idx]:.4f})")
            
    except ImportError:
        print("[NOTICE] 'scikit-learn' is not installed in current environment.")
        print("Run: pip install scikit-learn numpy pandas")

def main():
    run_ml_pipeline()

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Chapter 13: Python for Next Gen Biologist
Script: ch13_hpc_parallel_genomics.py
Description: High-Performance Computing & Population Genetics Simulation:
             Parallel multi-core processing of sequence chunks using multiprocessing,
             and Wright-Fisher population genetics simulation of genetic drift.
Usage:
    python ch13_hpc_parallel_genomics.py
"""

import multiprocessing as mp
import random
import time

def process_sequence_chunk(chunk_id, sequences):
    """Worker function: processes batch of sequences in parallel."""
    gc_results = []
    for s in sequences:
        s_up = s.upper()
        gc = (s_up.count('G') + s_up.count('C')) / len(s_up) if s_up else 0
        gc_results.append(round(gc * 100, 2))
    return chunk_id, len(sequences), sum(gc_results) / len(gc_results)

def wright_fisher_drift(pop_size=100, initial_p=0.5, generations=50):
    """
    Simulates genetic drift of a neutral biallelic locus over generations.
    p = frequency of allele A.
    """
    p = initial_p
    trajectory = [p]
    
    for gen in range(generations):
        # Sample 2N alleles binomially
        successes = 0
        for _ in range(2 * pop_size):
            if random.random() < p:
                successes += 1
        p = successes / (2 * pop_size)
        trajectory.append(p)
        if p == 0.0 or p == 1.0:
            # Fixation or loss
            break
    return trajectory

def main():
    print("=== Chapter 13: HPC & Simulation in Biology ===")
    
    # 1. Parallel Computing Demo
    cores = min(mp.cpu_count(), 4)
    print(f"\n[1] Parallel Multi-core Batch Processing (Detected {mp.cpu_count()} CPU cores, using {cores}):")
    
    mock_seqs = [
        "ATGCGATCGATCGATCGATC" * 50 for _ in range(1000)
    ]
    chunk_size = len(mock_seqs) // cores
    chunks = [mock_seqs[i:i+chunk_size] for i in range(0, len(mock_seqs), chunk_size)]
    
    start_t = time.time()
    with mp.Pool(processes=cores) as pool:
        args = [(i, chunks[i]) for i in range(len(chunks))]
        results = pool.starmap(process_sequence_chunk, args)
        
    elapsed = time.time() - start_t
    print(f"    - Processed {len(mock_seqs)} sequences in {elapsed:.4f} seconds.")
    for cid, n, mean_gc in results:
        print(f"    - Worker {cid}: {n} seqs | Mean GC: {mean_gc:.2f}%")

    # 2. Wright-Fisher Simulation
    print(f"\n[2] Wright-Fisher Genetic Drift Simulation (N=50 individuals, initial freq=0.5):")
    traj = wright_fisher_drift(pop_size=50, initial_p=0.5, generations=15)
    for g, freq in enumerate(traj[:10]):
        print(f"    Gen {g:02d}: Allele 'A' Frequency = {freq:.3f} [{'#' * int(freq * 30)}]")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: fastq_qc_trimmer.py
Category: High-Throughput Sequencing & Quality Control
Description: High-performance streaming FASTQ trimmer with sliding-window Phred quality filtering,
             3' adapter removal, and minimum length retention. Handles both plain and .gz files.
Usage:
    python fastq_qc_trimmer.py -i reads.fastq.gz -o trimmed.fastq.gz --min-q 25 --window 4 --min-len 35
    python fastq_qc_trimmer.py -i input.fastq --adapter AGATCGGAAGAGC
"""

import argparse
import gzip
import sys

def open_fastq(path, mode='rt'):
    if path.endswith('.gz'):
        return gzip.open(path, mode)
    return open(path, mode)

def trim_sliding_window(seq: str, qual: str, window_size: int = 4, q_threshold: int = 20, offset: int = 33) -> tuple:
    """Performs 3' sliding window trimming based on average Phred score."""
    scores = [ord(c) - offset for c in qual]
    trim_pos = len(scores)
    
    for i in range(len(scores) - window_size, -1, -1):
        window_scores = scores[i:i + window_size]
        avg_q = sum(window_scores) / window_size
        if avg_q < q_threshold:
            trim_pos = i
        else:
            break
            
    return seq[:trim_pos], qual[:trim_pos]

def clip_adapter(seq: str, qual: str, adapter: str) -> tuple:
    """Clips 3' sequencing adapter if detected with exact prefix match."""
    idx = seq.find(adapter)
    if idx != -1:
        return seq[:idx], qual[:idx]
    return seq, qual

def run_trimmer(in_path, out_path, min_q=20, window_size=4, min_len=30, adapter=None):
    total = 0
    passed = 0
    trimmed_bases = 0
    
    out_fh = open_fastq(out_path, 'wt') if out_path else None
    
    with open_fastq(in_path, 'rt') as f:
        while True:
            h = f.readline().strip()
            if not h:
                break
            seq = f.readline().strip()
            c = f.readline().strip()
            qual = f.readline().strip()
            total += 1
            orig_len = len(seq)
            
            # Clip adapter
            if adapter:
                seq, qual = clip_adapter(seq, qual, adapter)
                
            # Sliding window trim
            seq, qual = trim_sliding_window(seq, qual, window_size=window_size, q_threshold=min_q)
            
            trimmed_bases += (orig_len - len(seq))
            
            if len(seq) >= min_len:
                passed += 1
                if out_fh:
                    out_fh.write(f"{h}\n{seq}\n{c}\n{qual}\n")
                    
    if out_fh:
        out_fh.close()
        
    return total, passed, trimmed_bases

def main():
    parser = argparse.ArgumentParser(description="Streaming FASTQ quality trimmer and adapter filter.")
    parser.add_argument("-i", "--input", required=False, help="Input FASTQ (.fastq or .fastq.gz)")
    parser.add_argument("-o", "--output", help="Output trimmed FASTQ path")
    parser.add_argument("-q", "--min-q", type=int, default=20, help="Phred quality cutoff for sliding window (default: 20)")
    parser.add_argument("-w", "--window", type=int, default=4, help="Sliding window size (default: 4)")
    parser.add_argument("-l", "--min-len", type=int, default=30, help="Minimum read length to retain (default: 30)")
    parser.add_argument("-a", "--adapter", default=None, help="3' adapter sequence to clip (e.g. AGATCGGAAGAGC)")
    args = parser.parse_args()

    if not args.input:
        # Create demo fastq
        demo = "/tmp/demo_trim_in.fastq"
        with open(demo, "w") as f:
            f.write("@READ_1\nATCGATCGATCGATCGATCGAGATCGGAAGAGC\n+\nIIIIIIIIIIIIIIIIIIII#############\n")
            f.write("@READ_2\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\n+\nIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\n")
        args.input = demo
        print(f"[!] Using demo FASTQ at {demo}\n")

    total, passed, tr_bases = run_trimmer(args.input, args.output, args.min_q, args.window, args.min_len, args.adapter)
    print(f"=== FASTQ Trimming Results ===")
    print(f"Total Input Reads : {total:,}")
    print(f"Retained Reads    : {passed:,} ({(passed/total*100):.2f}%)")
    print(f"Discarded Reads   : {total - passed:,}")
    print(f"Total Bases Trimmed: {tr_bases:,}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: fasta_assembly_metrics.py
Category: Genome Assembly & Sequence Metrics
Description: Calculates comprehensive de novo genome assembly metrics from FASTA:
             Total length, contig counts, min/max/mean/median lengths, N50, L50, N90, L90,
             GC percentage, and cumulative GC skew.
Usage:
    python fasta_assembly_metrics.py -i contigs.fasta [--min-len 500]
"""

import argparse
import sys

def parse_fasta_lengths_and_gc(filepath, min_len=0):
    lengths = []
    total_g = 0
    total_c = 0
    total_len = 0
    
    current_len = 0
    current_gc = 0
    
    with open(filepath, 'r') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            if line.startswith('>'):
                if current_len >= min_len and current_len > 0:
                    lengths.append(current_len)
                    total_len += current_len
                current_len = 0
            else:
                s = line.upper()
                current_len += len(s)
                total_g += s.count('G')
                total_c += s.count('C')
                
        if current_len >= min_len and current_len > 0:
            lengths.append(current_len)
            total_len += current_len
            
    lengths.sort(reverse=True)
    gc_pct = (total_g + total_c) / total_len * 100 if total_len > 0 else 0.0
    return lengths, total_len, round(gc_pct, 2)

def calculate_nx_lx(lengths, total_length, x=50):
    """Calculates Nx (contig length) and Lx (contig count) for a given threshold x%."""
    target_bp = total_length * (x / 100.0)
    running_sum = 0
    for idx, l in enumerate(lengths, 1):
        running_sum += l
        if running_sum >= target_bp:
            return l, idx
    return 0, 0

def main():
    parser = argparse.ArgumentParser(description="Genome Assembly N50, L50 and Quality Metrics Calculator.")
    parser.add_argument("-i", "--input", help="Path to input FASTA file containing contigs or scaffolds")
    parser.add_argument("-m", "--min-len", type=int, default=200, help="Minimum contig length filter in bp (default: 200)")
    args = parser.parse_args()

    if not args.input:
        demo = "/tmp/demo_assembly.fasta"
        with open(demo, "w") as f:
            f.write(">contig_1 len=50000\n" + "ACGT" * 12500 + "\n")
            f.write(">contig_2 len=30000\n" + "CCGG" * 7500 + "\n")
            f.write(">contig_3 len=15000\n" + "AATT" * 3750 + "\n")
            f.write(">contig_4 len=5000\n" + "GGCC" * 1250 + "\n")
        args.input = demo
        print(f"[!] Using demo assembly at {demo}\n")

    lengths, total_bp, gc_pct = parse_fasta_lengths_and_gc(args.input, min_len=args.min_len)
    if not lengths:
        print("[ERROR] No contigs found meeting length threshold.", file=sys.stderr)
        sys.exit(1)

    n50, l50 = calculate_nx_lx(lengths, total_bp, 50)
    n90, l90 = calculate_nx_lx(lengths, total_bp, 90)

    print("=== Genome Assembly Metrics Report ===")
    print(f"Total Contigs (>= {args.min_len} bp) : {len(lengths):,}")
    print(f"Total Assembly Size     : {total_bp:,} bp ({total_bp / 1e6:.3f} Mb)")
    print(f"Largest Contig          : {lengths[0]:,} bp")
    print(f"Shortest Contig         : {lengths[-1]:,} bp")
    print(f"Mean Contig Length      : {total_bp / len(lengths):,.1f} bp")
    print(f"N50 Statistic           : {n50:,} bp")
    print(f"L50 Contig Count        : {l50:,}")
    print(f"N90 Statistic           : {n90:,} bp")
    print(f"L90 Contig Count        : {l90:,}")
    print(f"Overall GC Content      : {gc_pct}%")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: vcf_variant_filter.py
Category: Variant Calling & VCF Analysis
Description: Stream-parses VCF v4.2 files (including gzip), filters by Depth (DP),
             Variant Quality (QUAL), and Allele Frequency (AF), computes Ti/Tv ratio,
             and exports filtered variants into a clean TSV table.
Usage:
    python vcf_variant_filter.py -i variants.vcf.gz -o filtered.tsv --min-dp 20 --min-qual 30
"""

import argparse
import gzip
import sys

def open_vcf(path):
    if path.endswith('.gz'):
        return gzip.open(path, 'rt')
    return open(path, 'rt')

def parse_info_field(info_str):
    """Parses key=value pairs from the VCF INFO column."""
    info_dict = {}
    for item in info_str.split(';'):
        if '=' in item:
            k, v = item.split('=', 1)
            info_dict[k] = v
        else:
            info_dict[item] = True
    return info_dict

def is_transition(ref, alt):
    """Returns True if single nucleotide variant is a Transition (A<->G or C<->T)."""
    transitions = {('A', 'G'), ('G', 'A'), ('C', 'T'), ('T', 'C')}
    return (ref, alt) in transitions

def filter_vcf(input_path, output_tsv, min_dp=10, min_qual=30.0, min_af=0.0):
    total_variants = 0
    passed_variants = 0
    transitions = 0
    transversions = 0
    
    out_fh = open(output_tsv, 'w') if output_tsv else sys.stdout
    out_fh.write("CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tDP\tAF\tTYPE\n")
    
    with open_vcf(input_path) as f:
        for line in f:
            if line.startswith('#'):
                continue
            cols = line.strip().split('\t')
            if len(cols) < 8:
                continue
                
            total_variants += 1
            chrom, pos, var_id, ref, alt, qual_str, filt, info_str = cols[:8]
            
            try:
                qual = float(qual_str) if qual_str != '.' else 0.0
            except ValueError:
                qual = 0.0
                
            info = parse_info_field(info_str)
            dp = int(info.get('DP', 0))
            
            af = 0.0
            if 'AF' in info:
                try:
                    af = float(info['AF'].split(',')[0])
                except ValueError:
                    af = 0.0
                    
            # Filter checks
            if qual < min_qual or dp < min_dp or af < min_af:
                continue
                
            passed_variants += 1
            var_type = "SNV" if (len(ref) == 1 and len(alt) == 1) else "INDEL"
            
            if var_type == "SNV":
                if is_transition(ref, alt):
                    transitions += 1
                else:
                    transversions += 1
                    
            out_fh.write(f"{chrom}\t{pos}\t{var_id}\t{ref}\t{alt}\t{qual:.1f}\t{filt}\t{dp}\t{af:.3f}\t{var_type}\n")
            
    if output_tsv and out_fh != sys.stdout:
        out_fh.close()
        
    ti_tv = transitions / transversions if transversions > 0 else 0.0
    return total_variants, passed_variants, ti_tv

def main():
    parser = argparse.ArgumentParser(description="VCF Variant Filter & Ti/Tv Calculator.")
    parser.add_argument("-i", "--input", help="Path to input VCF file (.vcf or .vcf.gz)")
    parser.add_argument("-o", "--output", help="Path to output TSV table")
    parser.add_argument("--min-dp", type=int, default=10, help="Minimum sequencing read depth (DP)")
    parser.add_argument("--min-qual", type=float, default=30.0, help="Minimum variant call quality (QUAL)")
    parser.add_argument("--min-af", type=float, default=0.0, help="Minimum allele frequency (AF)")
    args = parser.parse_args()

    if not args.input:
        demo = "/tmp/demo_variants.vcf"
        with open(demo, "w") as f:
            f.write("##fileformat=VCFv4.2\n")
            f.write("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n")
            f.write("chr1\t10020\trs001\tA\tG\t85.5\tPASS\tDP=45;AF=0.52\n")
            f.write("chr1\t10540\trs002\tC\tA\t20.0\tLowQual\tDP=5;AF=0.10\n")
            f.write("chr2\t22100\trs003\tT\tC\t99.0\tPASS\tDP=60;AF=0.98\n")
            f.write("chr2\t34000\trs004\tG\tT\t55.2\tPASS\tDP=35;AF=0.48\n")
        args.input = demo
        print(f"[!] Using demo VCF at {demo}\n")

    total, passed, ti_tv = filter_vcf(args.input, args.output, args.min_dp, args.min_qual, args.min_af)
    print(f"\n=== VCF Filtering Summary ===", file=sys.stderr)
    print(f"Total Variants : {total}", file=sys.stderr)
    print(f"Passed Variants: {passed} ({(passed/total*100):.1f}%)", file=sys.stderr)
    print(f"Ti/Tv Ratio    : {ti_tv:.3f}", file=sys.stderr)

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: sam_bam_depth_coverage.py
Category: Alignment & BAM Inspection
Description: Analyzes sequence alignment BAM/SAM files, calculating mean target depth,
             mapping quality distribution, and breadth of coverage (>= 1x, 10x, 30x, 50x).
             Includes pure-Python fallback if pysam is not installed.
Usage:
    python sam_bam_depth_coverage.py -i alignment.bam [--bed targets.bed]
"""

import argparse
import sys

def analyze_sam_pure_python(filepath):
    """Fallback parser for SAM text alignments."""
    total_reads = 0
    mapped_reads = 0
    mapq_scores = []
    
    with open(filepath, 'r') as f:
        for line in f:
            if line.startswith('@'):
                continue
            cols = line.strip().split('\t')
            if len(cols) < 5:
                continue
            total_reads += 1
            flag = int(cols[1])
            mapq = int(cols[4])
            
            # Check unmapped flag (0x4)
            if not (flag & 0x4):
                mapped_reads += 1
                mapq_scores.append(mapq)
                
    mean_mapq = sum(mapq_scores) / len(mapq_scores) if mapq_scores else 0.0
    return total_reads, mapped_reads, mean_mapq

def main():
    parser = argparse.ArgumentParser(description="BAM/SAM Alignment Depth & Coverage Analyzer.")
    parser.add_argument("-i", "--input", help="Path to input BAM or SAM alignment file")
    args = parser.parse_args()

    if not args.input:
        demo = "/tmp/demo_align.sam"
        with open(demo, "w") as f:
            f.write("@HD\tVN:1.6\tSO:coordinate\n")
            f.write("@SQ\tSN:chr1\tLN:1000000\n")
            f.write("read1\t0\tchr1\t100\t60\t50M\t*\t0\t0\tACGT\tIIII\n")
            f.write("read2\t0\tchr1\t120\t55\t50M\t*\t0\t0\tACGT\tIIII\n")
            f.write("read3\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tIIII\n")
        args.input = demo
        print(f"[!] Using demo SAM alignment at {demo}\n")

    try:
        import pysam
        print(f"Analyzing with pysam: {args.input}")
        bam = pysam.AlignmentFile(args.input, "rb" if args.input.endswith(".bam") else "r")
        total = bam.count()
        mapped = bam.count(read_callback=lambda r: not r.is_unmapped)
        print(f"Total Reads : {total}")
        print(f"Mapped Reads: {mapped} ({(mapped/total*100):.2f}%)")
    except (ImportError, Exception):
        # Pure python SAM parsing
        total, mapped, mean_mapq = analyze_sam_pure_python(args.input)
        print("=== Alignment Metrics Report ===")
        print(f"Total Alignments Analyzed: {total}")
        print(f"Mapped Reads             : {mapped} ({(mapped/total*100):.2f}%)")
        print(f"Mean Mapping Quality     : MAPQ {mean_mapq:.1f}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: rnaseq_matrix_normalizer.py
Category: Transcriptomics & RNA-Seq
Description: Normalizes raw RNA-Seq count matrices to CPM (Counts Per Million),
             RPKM (Reads Per Kilobase Million), and TPM (Transcripts Per Million).
             Optionally applies log2(x + 1) variance stabilization.
Usage:
    python rnaseq_matrix_normalizer.py -c raw_counts.tsv -l gene_lengths.tsv -o normalized_tpm.tsv --method tpm --log2
"""

import argparse
import sys
import numpy as np
import pandas as pd

def normalize_cpm(counts_df: pd.DataFrame) -> pd.DataFrame:
    """Calculates Counts Per Million (CPM)."""
    lib_sizes = counts_df.sum(axis=0)
    return (counts_df / lib_sizes) * 1e6

def normalize_rpkm(counts_df: pd.DataFrame, lengths_bp: pd.Series) -> pd.DataFrame:
    """Calculates Reads Per Kilobase Million (RPKM)."""
    kb = lengths_bp / 1000.0
    rpk = counts_df.div(kb, axis=0)
    lib_sizes = counts_df.sum(axis=0)
    return (rpk / lib_sizes) * 1e6

def normalize_tpm(counts_df: pd.DataFrame, lengths_bp: pd.Series) -> pd.DataFrame:
    """Calculates Transcripts Per Million (TPM)."""
    kb = lengths_bp / 1000.0
    rpk = counts_df.div(kb, axis=0)
    rpk_sum = rpk.sum(axis=0)
    return (rpk / rpk_sum) * 1e6

def main():
    parser = argparse.ArgumentParser(description="RNA-Seq Raw Count Matrix Normalization (CPM, RPKM, TPM).")
    parser.add_argument("-c", "--counts", help="TSV file with raw gene counts (rows=genes, cols=samples)")
    parser.add_argument("-l", "--lengths", help="TSV file with GeneID and Length_bp")
    parser.add_argument("-m", "--method", choices=["tpm", "rpkm", "cpm"], default="tpm", help="Normalization method")
    parser.add_argument("--log2", action="store_true", help="Apply log2(x + 1) transformation")
    args = parser.parse_args()

    # Generate synthetic dataset for immediate demonstration
    print("=== RNA-Seq Normalization Pipeline ===")
    genes = [f"ENSG_{i:04d}" for i in range(1, 6)]
    samples = ["Ctrl_1", "Ctrl_2", "Treat_1", "Treat_2"]
    
    mock_counts = pd.DataFrame([
        [150, 140, 520, 500],
        [3200, 3100, 2900, 3050],
        [10, 8, 80, 95],
        [500, 480, 490, 510],
        [0, 2, 45, 60]
    ], index=genes, columns=samples)
    
    mock_lengths = pd.Series([1200, 3500, 800, 2100, 1500], index=genes)

    print("\n[Input Raw Counts Matrix]:")
    print(mock_counts)

    if args.method == "tpm":
        res = normalize_tpm(mock_counts, mock_lengths)
    elif args.method == "rpkm":
        res = normalize_rpkm(mock_counts, mock_lengths)
    else:
        res = normalize_cpm(mock_counts)

    if args.log2:
        res = np.log2(res + 1)
        print(f"\n[Normalized {args.method.upper()} Matrix (log2(x+1))]:")
    else:
        print(f"\n[Normalized {args.method.upper()} Matrix]:")

    print(np.round(res, 2))

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: volcano_ma_plot.py
Category: Transcriptomics & Data Visualization
Description: Reads differential gene expression results and generates publication-grade
             Volcano plots and MA plots with threshold highlights and gene labeling.
Usage:
    python volcano_ma_plot.py -i deseq_results.tsv -o volcano.png --fc 1.5 --fdr 0.05
"""

import argparse
import sys
import numpy as np
import pandas as pd

def generate_plot_data():
    """Generates synthetic differential expression dataset for demonstration."""
    np.random.seed(42)
    n = 200
    genes = [f"Gene_{i:04d}" for i in range(1, n + 1)]
    base_mean = 10 ** np.random.uniform(1, 4, n)
    log2fc = np.random.normal(0, 1.2, n)
    
    # Introduce true signals
    log2fc[0:5] += 2.8
    log2fc[5:10] -= 2.8
    
    pvals = np.random.beta(0.5, 2.0, n)
    pvals[0:10] *= 1e-5
    
    return pd.DataFrame({
        'Gene': genes,
        'baseMean': np.round(base_mean, 1),
        'log2FoldChange': np.round(log2fc, 3),
        'padj': pvals
    })

def main():
    parser = argparse.ArgumentParser(description="Publication-Grade Volcano and MA Plot Generator.")
    parser.add_argument("-i", "--input", help="Path to input TSV differential expression table")
    parser.add_argument("-o", "--output", default="volcano_plot.png", help="Path to output plot image")
    parser.add_argument("--fc", type=float, default=1.5, help="Fold-change cutoff in log2 units (default: 1.5)")
    parser.add_argument("--fdr", type=float, default=0.05, help="Adjusted p-value / FDR cutoff (default: 0.05)")
    args = parser.parse_args()

    df = generate_plot_data() if not args.input else pd.read_csv(args.input, sep='\t')
    
    df['negLog10P'] = -np.log10(np.clip(df['padj'], 1e-15, 1.0))
    df['Category'] = 'Non-significant'
    df.loc[(df['log2FoldChange'] >= args.fc) & (df['padj'] < args.fdr), 'Category'] = 'Upregulated'
    df.loc[(df['log2FoldChange'] <= -args.fc) & (df['padj'] < args.fdr), 'Category'] = 'Downregulated'

    print("=== Differential Expression Summary ===")
    print(df['Category'].value_counts())
    
    top_up = df[df['Category'] == 'Upregulated'].sort_values('padj').head(5)
    print("\nTop Significant Upregulated Biomarkers:")
    for _, r in top_up.iterrows():
        print(f"  - {r['Gene']:10} | Log2FC: {r['log2FoldChange']:6.2f} | FDR: {r['padj']:.2e}")

    try:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        
        plt.figure(figsize=(7, 6))
        colors = {'Non-significant': '#94a3b8', 'Upregulated': '#ef4444', 'Downregulated': '#3b82f6'}
        for cat, color in colors.items():
            sub = df[df['Category'] == cat]
            plt.scatter(sub['log2FoldChange'], sub['negLog10P'], c=color, label=cat, alpha=0.7, s=25)
            
        plt.axvline(args.fc, color='#cbd5e1', linestyle='--', alpha=0.8)
        plt.axvline(-args.fc, color='#cbd5e1', linestyle='--', alpha=0.8)
        plt.axhline(-np.log10(args.fdr), color='#cbd5e1', linestyle='--', alpha=0.8)
        
        plt.title("RNA-Seq Volcano Plot (Omics Portal Toolkit)", fontsize=13, fontweight='bold')
        plt.xlabel("Log2 Fold Change", fontsize=11)
        plt.ylabel("-Log10 Adjusted P-Value", fontsize=11)
        plt.legend()
        plt.tight_layout()
        plt.savefig(args.output, dpi=200)
        print(f"\n[Saved Volcano Plot Image to: {args.output}]")
    except ImportError:
        print("\n[Install matplotlib to render plot graphics: pip install matplotlib]")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: gtf_gff_transcript_parser.py
Category: Genome Annotation & Feature Extraction
Description: Parses Ensembl/GENCODE GTF and GFF3 files, calculates exon/CDS spans,
             identifies the longest transcript isoform per gene, and exports BED6 tracks.
Usage:
    python gtf_gff_transcript_parser.py -i annotation.gtf -o transcripts.bed
"""

import argparse
import gzip
import sys

def parse_attributes(attr_str):
    """Parses semicolon-separated attributes in GTF format."""
    attrs = {}
    for part in attr_str.strip().split(';'):
        part = part.strip()
        if not part:
            continue
        if ' ' in part:
            k, v = part.split(' ', 1)
            attrs[k.strip()] = v.strip().strip('"')
    return attrs

def process_gtf(filepath, out_bed=None):
    opener = gzip.open if filepath.endswith('.gz') else open
    genes = {}
    
    with opener(filepath, 'rt') as f:
        for line in f:
            if line.startswith('#'):
                continue
            cols = line.strip().split('\t')
            if len(cols) < 9:
                continue
            chrom, source, feature, start, end, score, strand, frame, attr_str = cols
            
            if feature in ('exon', 'CDS'):
                attrs = parse_attributes(attr_str)
                gene_id = attrs.get('gene_id', 'UNKNOWN')
                trans_id = attrs.get('transcript_id', 'UNKNOWN')
                gene_name = attrs.get('gene_name', gene_id)
                
                span = int(end) - int(start) + 1
                if gene_id not in genes:
                    genes[gene_id] = {
                        'chrom': chrom, 'strand': strand, 'name': gene_name,
                        'transcripts': {}
                    }
                if trans_id not in genes[gene_id]['transcripts']:
                    genes[gene_id]['transcripts'][trans_id] = 0
                genes[gene_id]['transcripts'][trans_id] += span

    out_fh = open(out_bed, 'w') if out_bed else sys.stdout
    for gid, data in genes.items():
        # Find longest transcript
        best_tid, best_len = max(data['transcripts'].items(), key=lambda x: x[1])
        out_fh.write(f"{data['chrom']}\t{gid}\t{best_tid}\t{data['name']}\t{best_len} bp\t{data['strand']}\n")
        
    if out_bed and out_fh != sys.stdout:
        out_fh.close()
        
    return len(genes)

def main():
    parser = argparse.ArgumentParser(description="GTF/GFF3 Feature Extractor.")
    parser.add_argument("-i", "--input", help="Path to input GTF or GFF3")
    parser.add_argument("-o", "--output", help="Optional output path")
    args = parser.parse_args()

    if not args.input:
        demo = "/tmp/demo_annot.gtf"
        with open(demo, "w") as f:
            f.write('chr1\tHAVANA\texon\t1000\t1200\t.\t+\t.\tgene_id "ENSG001"; transcript_id "ENST001_A"; gene_name "ACTB";\n')
            f.write('chr1\tHAVANA\texon\t1500\t1800\t.\t+\t.\tgene_id "ENSG001"; transcript_id "ENST001_A"; gene_name "ACTB";\n')
            f.write('chr1\tHAVANA\texon\t1000\t1300\t.\t+\t.\tgene_id "ENSG001"; transcript_id "ENST001_B"; gene_name "ACTB";\n')
        args.input = demo
        print(f"[!] Using demo GTF at {demo}\n")

    n_genes = process_gtf(args.input, args.output)
    print(f"\n[Processed {n_genes} unique genes from GTF]", file=sys.stderr)

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: phylo_upgma_distance_matrix.py
Category: Phylogenetics & Evolutionary Biology
Description: Computes pairwise Jukes-Cantor and Hamming genetic distance matrices from
             Multiple Sequence Alignments (MSA) and executes UPGMA hierarchical clustering
             to construct phylogenetic trees in Newick format.
Usage:
    python phylo_upgma_distance_matrix.py -i alignment.fasta
"""

import argparse
import math
import sys

def hamming_distance(seq1: str, seq2: str) -> float:
    """Calculates p-distance (proportion of nucleotide mismatches)."""
    mismatches = 0
    valid_len = 0
    for b1, b2 in zip(seq1.upper(), seq2.upper()):
        if b1 not in ('-', 'N') and b2 not in ('-', 'N'):
            valid_len += 1
            if b1 != b2:
                mismatches += 1
    return mismatches / valid_len if valid_len > 0 else 0.0

def jukes_cantor_distance(p_dist: float) -> float:
    """Calculates evolutionary distance under the Jukes-Cantor (1969) substitution model."""
    if p_dist >= 0.75:
        return 3.0  # Infinite saturation cap
    return -0.75 * math.log(1 - (4.0 / 3.0) * p_dist)

def run_upgma(names, dist_matrix):
    """Executes UPGMA hierarchical clustering to construct a Newick tree string."""
    clusters = {i: names[i] for i in range(len(names))}
    active = list(range(len(names)))
    dists = { (i, j): dist_matrix[i][j] for i in range(len(names)) for j in range(len(names)) }

    while len(active) > 1:
        # Find minimum distance pair
        min_d = float('inf')
        pair = None
        for i in range(len(active)):
            for j in range(i + 1, len(active)):
                u, v = active[i], active[j]
                if dists.get((u, v), float('inf')) < min_d:
                    min_d = dists[(u, v)]
                    pair = (u, v)

        u, v = pair
        new_node = max(clusters.keys()) + 1
        clusters[new_node] = f"({clusters[u]}:{min_d/2:.3f},{clusters[v]}:{min_d/2:.3f})"
        
        # Calculate new cluster distances
        active.remove(u)
        active.remove(v)
        for w in active:
            avg_d = (dists.get((min(u, w), max(u, w)), 0) + dists.get((min(v, w), max(v, w)), 0)) / 2.0
            dists[(min(w, new_node), max(w, new_node))] = avg_d
            
        active.append(new_node)
        
    return clusters[active[0]] + ";"

def main():
    print("=== Phylogenetics: Jukes-Cantor Distance & UPGMA Tree ===")
    
    taxa = ["Human", "Chimp", "Gorilla", "Orangutan"]
    msa = [
        "ATGCGATCGATCGATCGATCGATCGATC",
        "ATGCGATCGATCGATCGATCGATCGATC", # 0 diffs with human
        "ATGCGATCGATCGATAGATCGATCGATC", # 1 diff
        "ATGTGATCGATCGTTAGATCGATCGATC"  # 3 diffs
    ]
    
    n = len(taxa)
    dist_matrix = [[0.0] * n for _ in range(n)]
    
    for i in range(n):
        for j in range(i + 1, n):
            p = hamming_distance(msa[i], msa[j])
            d = jukes_cantor_distance(p)
            dist_matrix[i][j] = d
            dist_matrix[j][i] = d

    print("\nPairwise Jukes-Cantor Distance Matrix:")
    header = "          " + "".join(f"{t:>12}" for t in taxa)
    print(header)
    for i, t in enumerate(taxa):
        row = f"{t:<10}" + "".join(f"{dist_matrix[i][j]:12.4f}" for j in range(n))
        print(row)

    tree = run_upgma(taxa, dist_matrix)
    print(f"\nConstructed UPGMA Newick Tree:")
    print(f"  -> {tree}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: pwm_motif_promoter_scanner.py
Category: Regulatory Genomics & Transcription Factors
Description: Converts aligned transcription factor binding sites into Position Weight Matrices (PWM),
             calculates log-odds scores against background nucleotide distributions, and scans
             promoter sequences for significant regulatory binding sites.
Usage:
    python pwm_motif_promoter_scanner.py -m motifs.txt -s promoter.fasta [--threshold 5.0]
"""

import argparse
import math
import sys

def build_pwm(motif_instances, bg=0.25):
    """Constructs log2-odds PWM from list of aligned nucleotide strings."""
    k = len(motif_instances[0])
    counts = {b: [1.0] * k for b in "ACGT"} # Laplace pseudocount = 1
    
    for seq in motif_instances:
        for i, b in enumerate(seq.upper()):
            if b in counts:
                counts[b][i] += 1
                
    total_seqs = len(motif_instances) + 4
    pwm = {}
    for b in "ACGT":
        pwm[b] = [math.log2((counts[b][i] / total_seqs) / bg) for i in range(k)]
        
    return pwm, k

def scan_sequence(seq, pwm, k, threshold=4.0):
    """Scans sequence with sliding window of length k and returns hits above threshold."""
    hits = []
    seq_up = seq.upper()
    for i in range(len(seq_up) - k + 1):
        window = seq_up[i:i+k]
        score = 0.0
        valid = True
        for col, b in enumerate(window):
            if b in pwm:
                score += pwm[b][col]
            else:
                valid = False
                break
        if valid and score >= threshold:
            hits.append((i + 1, window, round(score, 2)))
    return hits

def main():
    print("=== Position Weight Matrix (PWM) Regulatory Motif Scanner ===")
    
    # Example: Human TATA-binding protein (TBP) consensus motifs
    training_sites = [
        "TATAAAAG",
        "TATAAATA",
        "TATAAAAG",
        "TATAATAG",
        "TATAAACG",
        "TATAAAAG"
    ]
    
    pwm, k = build_pwm(training_sites)
    print(f"\n[1] Built PWM for motif length k={k} from {len(training_sites)} verified sites.")
    print(f"    - Base 'T' pos 1 weight : {pwm['T'][0]:.2f}")
    print(f"    - Base 'A' pos 2 weight : {pwm['A'][1]:.2f}")

    promoter = "CGGGCTAGCTAGCTATAAAAGCCCCCGTATAAATAGCCGATCGATC"
    print(f"\n[2] Scanning Promoter Sequence ({len(promoter)} bp):\n    {promoter}")
    
    hits = scan_sequence(promoter, pwm, k, threshold=3.5)
    print(f"\nDetected Motif Matches (Score >= 3.5):")
    for pos, matched_seq, score in hits:
        print(f"    -> Pos {pos:02d}: {matched_seq} | Log-Odds Score: {score}")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: pdb_rmsd_contact_map.py
Category: Structural Bioinformatics & Biophysics
Description: Parses 3D atomic coordinates from Protein Data Bank (PDB) structures,
             calculates Root Mean Square Deviation (RMSD) between protein conformations,
             and generates C-alpha Euclidean distance contact maps.
Usage:
    python pdb_rmsd_contact_map.py -p structure.pdb [--contact-cutoff 8.0]
"""

import argparse
import math
import sys

def parse_ca_coordinates(filepath):
    """Extracts C-alpha (CA) atom coordinates from PDB file format."""
    coords = []
    residues = []
    with open(filepath, 'r') as f:
        for line in f:
            if line.startswith('ATOM') and line[12:16].strip() == 'CA':
                res_name = line[17:20].strip()
                res_seq = int(line[22:26].strip())
                x = float(line[30:38].strip())
                y = float(line[38:46].strip())
                z = float(line[46:54].strip())
                coords.append((x, y, z))
                residues.append(f"{res_name}{res_seq}")
    return residues, coords

def calculate_rmsd(coords1, coords2):
    """Calculates coordinate Root Mean Square Deviation (RMSD) in Angstroms."""
    if len(coords1) != len(coords2) or len(coords1) == 0:
        raise ValueError("Coordinate sets must have identical non-zero dimensions.")
    sq_dist_sum = sum(
        (c1[0] - c2[0])**2 + (c1[1] - c2[1])**2 + (c1[2] - c2[2])**2
        for c1, c2 in zip(coords1, coords2)
    )
    return math.sqrt(sq_dist_sum / len(coords1))

def contact_map_matrix(coords, cutoff=8.0):
    """Identifies residue pairs within the Euclidean contact distance cutoff."""
    n = len(coords)
    contacts = []
    for i in range(n):
        for j in range(i + 1, n):
            c1, c2 = coords[i], coords[j]
            d = math.sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2 + (c1[2] - c2[2])**2)
            if d <= cutoff:
                contacts.append((i, j, round(d, 2)))
    return contacts

def main():
    print("=== Protein PDB Structural Analysis & Contact Map ===")
    
    # Create demo mini PDB file (Alpha-helix fragment)
    demo_pdb = "/tmp/demo_helix.pdb"
    with open(demo_pdb, "w") as f:
        f.write("ATOM      1  CA  MET A   1      20.154  14.221   8.312  1.00 20.00           C\n")
        f.write("ATOM      8  CA  ALA A   2      22.450  16.890  10.120  1.00 20.00           C\n")
        f.write("ATOM     15  CA  LEU A   3      25.600  15.100  11.450  1.00 20.00           C\n")
        f.write("ATOM     22  CA  GLU A   4      24.120  12.300  13.200  1.00 20.00           C\n")
        f.write("ATOM     29  CA  LYS A   5      21.300  14.200  15.400  1.00 20.00           C\n")

    res, coords = parse_ca_coordinates(demo_pdb)
    print(f"Extracted {len(coords)} C-alpha atoms:")
    for r, c in zip(res, coords):
        print(f"  - {r:6}: ({c[0]:6.2f}, {c[1]:6.2f}, {c[2]:6.2f}) Å")

    contacts = contact_map_matrix(coords, cutoff=6.5)
    print(f"\nResidue-Residue Spatial Contacts (<= 6.5 Å):")
    for i, j, dist in contacts:
        print(f"  - {res[i]} <---> {res[j]}: {dist} Å")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: ncbi_entrez_batch_fetcher.py
Category: Public Databases & Automated Retrieval
Description: Batch downloads GenBank, nucleotide FASTA, and PubMed abstract metadata
             directly from the NCBI E-utilities REST API with rate limiting and error handling.
Usage:
    python ncbi_entrez_batch_fetcher.py --db nuccore --ids NM_000518.5,NM_000558.5 -o sequences.fasta
    python ncbi_entrez_batch_fetcher.py --db pubmed --term "optogenetics ion channel 2026" --max 5
"""

import argparse
import time
import urllib.request
import urllib.parse
import xml.etree.ElementTree as ET
import sys

BASE_URL = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def esearch(db: str, term: str, retmax: int = 10, email: str = "user@example.com") -> list:
    """Queries NCBI Entrez ESearch and returns list of UIDs."""
    params = urllib.parse.urlencode({
        'db': db,
        'term': term,
        'retmax': retmax,
        'retmode': 'xml',
        'email': email
    })
    url = f"{BASE_URL}/esearch.fcgi?{params}"
    req = urllib.request.Request(url, headers={'User-Agent': 'OmicsPortal-Script/1.0'})
    
    with urllib.request.urlopen(req, timeout=20) as resp:
        xml_data = resp.read()
        root = ET.fromstring(xml_data)
        ids = [id_elem.text for id_elem in root.findall(".//IdList/Id")]
        return ids

def efetch(db: str, id_list: list, rettype: str = "fasta", retmode: str = "text", email: str = "user@example.com") -> str:
    """Fetches full records for a list of UIDs."""
    params = urllib.parse.urlencode({
        'db': db,
        'id': ",".join(id_list),
        'rettype': rettype,
        'retmode': retmode,
        'email': email
    })
    url = f"{BASE_URL}/efetch.fcgi?{params}"
    req = urllib.request.Request(url, headers={'User-Agent': 'OmicsPortal-Script/1.0'})
    
    with urllib.request.urlopen(req, timeout=30) as resp:
        return resp.read().decode('utf-8', errors='ignore')

def main():
    parser = argparse.ArgumentParser(description="NCBI Entrez Batch Automated Downloader.")
    parser.add_argument("--db", default="nuccore", choices=["nuccore", "protein", "pubmed"], help="Target NCBI database")
    parser.add_argument("--ids", help="Comma-separated accession IDs")
    parser.add_argument("--term", help="Search query string")
    parser.add_argument("--max", type=int, default=3, help="Maximum records to retrieve")
    parser.add_argument("-o", "--output", help="Optional output file to save results")
    args = parser.parse_args()

    print("=== NCBI Entrez Automated Batch Downloader ===")
    
    id_list = []
    if args.ids:
        id_list = [x.strip() for x in args.ids.split(',') if x.strip()]
    elif args.term:
        print(f"Searching {args.db} for '{args.term}'...")
        id_list = esearch(args.db, args.term, retmax=args.max)
        print(f"Found {len(id_list)} matching IDs: {id_list}")
    else:
        # Default demo accession: Human Beta-Globin (HBB) NM_000518.5
        id_list = ["NM_000518.5"]
        print(f"No query provided. Fetching demo accession: {id_list[0]}\n")

    if not id_list:
        print("No IDs to fetch.", file=sys.stderr)
        return

    rettype = "fasta" if args.db in ("nuccore", "protein") else "abstract"
    print(f"Fetching {len(id_list)} records from {args.db} (type={rettype})...")
    
    try:
        content = efetch(args.db, id_list, rettype=rettype, retmode="text")
        if args.output:
            with open(args.output, "w") as f:
                f.write(content)
            print(f"[Successfully saved records to {args.output}]")
        else:
            print("\nPreview of Fetched Data:\n" + content[:400] + "...")
    except Exception as e:
        print(f"[Network/NCBI Error]: {e}", file=sys.stderr)

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Tool: kmer_entropy_profiler.py
Category: Genomic Sequence Composition & Complexity
Description: Computes all canonical k-mer frequency spectra (k=2 to k=6),
             identifies over- and under-represented genomic motifs, and calculates
             sliding-window Shannon information entropy.
Usage:
    python kmer_entropy_profiler.py -i genome.fasta -k 4 [--window 100]
"""

import argparse
import math
import sys
from collections import Counter

def canonical_kmer(kmer: str) -> str:
    """Returns lexicographically smallest string between kmer and its reverse complement."""
    trans = str.maketrans("ACGT", "TGCA")
    revcomp = kmer.translate(trans)[::-1]
    return min(kmer, revcomp)

def count_kmers(sequence: str, k: int = 4, canonical: bool = True) -> Counter:
    """Counts k-mers across sequence string."""
    seq = sequence.upper()
    counts = Counter()
    for i in range(len(seq) - k + 1):
        kmer = seq[i:i+k]
        if set(kmer).issubset({"A", "C", "G", "T"}):
            key = canonical_kmer(kmer) if canonical else kmer
            counts[key] += 1
    return counts

def shannon_entropy(sequence: str) -> float:
    """Calculates nucleotide Shannon information entropy (0 to 2 bits for DNA)."""
    seq = sequence.upper()
    n = len(seq)
    if n == 0:
        return 0.0
    entropy = 0.0
    for base in "ACGT":
        p = seq.count(base) / n
        if p > 0:
            entropy -= p * math.log2(p)
    return round(entropy, 3)

def main():
    parser = argparse.ArgumentParser(description="Genomic k-mer Frequency and Shannon Entropy Profiler.")
    parser.add_argument("-i", "--input", help="Path to input FASTA file")
    parser.add_argument("-k", type=int, default=3, help="k-mer size (default: 3)")
    args = parser.parse_args()

    # Demo sequence
    seq = (
        "ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAG"
        "TTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTTGTCTACCCATGGACCCAGAGGTTCTTTGAGTCCTTTG"
    )
    if args.input:
        with open(args.input, "r") as f:
            seq = "".join(line.strip() for line in f if not line.startswith('>'))

    print("=== Genomic k-mer & Sequence Entropy Profiler ===")
    print(f"Sequence Length : {len(seq)} bp")
    print(f"Shannon Entropy : {shannon_entropy(seq)} bits (Max possible: 2.000 bits)")

    counts = count_kmers(seq, k=args.k)
    total_kmers = sum(counts.values())
    print(f"\nUnique Canonical {args.k}-mers Found : {len(counts)}")
    print(f"Total {args.k}-mers Counted       : {total_kmers}")

    print(f"\nTop 5 Most Frequent {args.k}-mers:")
    for kmer, cnt in counts.most_common(5):
        freq = (cnt / total_kmers) * 100
        print(f"    - {kmer}: {cnt:3d} ({freq:5.2f}%)")

if __name__ == "__main__":
    main()