Genetics: From Mendel to Modern Genomics

# Load required packages
using Random, Statistics, Plots, StatsBase
using LinearAlgebra, DataFrames

# Set random seed for reproducibility
Random.seed!(42)

println(" Julia Genetics Analysis Environment Ready")
println("Version: ", VERSION)

# DNA base pairing rules
const BASE_PAIRS = Dict('A' => 'T', 'T' => 'A', 
                       'G' => 'C', 'C' => 'G')

# Generate complementary strand
function complement_strand(dna::String)
    return join([BASE_PAIRS[base] for base in uppercase(dna)])
end

# Example
template = "ATCGATCG"
comp = complement_strand(template)
println("Template:    5'-$template-3'")
println("Complement:  3'-$comp-5'")

# Analyze DNA composition
function analyze_dna(sequence::String)
    seq = uppercase(sequence)
    bases = Dict(b => count(==(b), seq) for b in "ATGC")
    total = length(seq)
    
    gc_content = 100.0 * (bases['G'] + bases['C']) / total
    
    println("Base Composition:")
    for b in "ATGC"
        pct = round(100 * bases[b] / total, digits=1)
        println("  $b: $(bases[b]) ($pct%)")
    end
    println("\nGC Content: $(round(gc_content, digits=1))%")
    
    # Verify Chargaff's rules
    println("\nChargaff's Rules:")
    println("  A/T ratio: $(round(bases['A']/bases['T'], digits=2))")
    println("  G/C ratio: $(round(bases['G']/bases['C'], digits=2))")
end

# Human beta-globin gene fragment
hbb_fragment = "ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTG"
analyze_dna(hbb_fragment)

# Transcription: DNA to RNA
function transcribe(dna::String)
    # RNA uses U instead of T
    return replace(uppercase(dna), 'T' => 'U')
end

# Genetic code (codon table)
const CODON_TABLE = Dict(
    "UUU"=>"F", "UUC"=>"F", "UUA"=>"L", "UUG"=>"L",
    "UCU"=>"S", "UCC"=>"S", "UCA"=>"S", "UCG"=>"S",
    "UAU"=>"Y", "UAC"=>"Y", "UAA"=>"*", "UAG"=>"*",
    "UGU"=>"C", "UGC"=>"C", "UGA"=>"*", "UGG"=>"W",
    "AUG"=>"M"  # Start codon
    # ... simplified for brevity
)

# Translation: RNA to Protein
function translate(rna::String)
    protein = Char[]
    
    for i in 1:3:length(rna)-2
        codon = rna[i:i+2]
        aa = get(CODON_TABLE, codon, '?')
        
        if aa == '*'  # Stop codon
            break
        else
            push!(protein, aa[1])
        end
    end
    
    return String(protein)
end

# Example: insulin gene fragment
dna_seq = "ATGGCCCTGTGGATGCGC"
rna_seq = transcribe(dna_seq)
protein_seq = translate(rna_seq)

println("DNA:     $dna_seq")
println("RNA:     $rna_seq")
println("Protein: $protein_seq")

# Simulate Mendelian inheritance
struct Genotype
    allele1::Char
    allele2::Char
end

# Determine phenotype (A = dominant)
phenotype(g::Genotype) = ('A' in [g.allele1, g.allele2]) ? "Dominant" : "Recessive"

# Gamete formation (meiosis)
function make_gamete(parent::Genotype)
    rand() < 0.5 ? parent.allele1 : parent.allele2
end

# Simulate genetic cross
function cross(parent1::Genotype, parent2::Genotype, n=1000)
    offspring = Genotype[]
    
    for _ in 1:n
        g1 = make_gamete(parent1)
        g2 = make_gamete(parent2)
        push!(offspring, Genotype(g1, g2))
    end
    
    return offspring
end

# Classic F1 × F1 cross (Aa × Aa)
p1 = Genotype('A', 'a')
p2 = Genotype('A', 'a')
f2_generation = cross(p1, p2, 10000)

# Count genotypes
genotypes = ["$(o.allele1)$(o.allele2)" for o in f2_generation]
geno_counts = countmap(genotypes)

println("F2 Generation (n=10000):")
for (g, count) in sort(collect(geno_counts))
    println("  $g: $count ($(round(100count/10000, digits=1))%)")
end

# Visualize Mendelian ratios
phenos = [phenotype(o) for o in f2_generation]
pheno_counts = countmap(phenos)

# Create visualization
bar(["Dominant", "Recessive"], 
    [pheno_counts["Dominant"], pheno_counts["Recessive"]],
    title="Phenotypic Ratio (Aa × Aa)",
    ylabel="Count",
    color=[:blue, :orange],
    legend=false)

# Add expected 3:1 ratio line
hline!([7500], label="Expected 3:1", linestyle=:dash, color=:red)

# Simulate genetic disease inheritance (recessive)
function disease_risk(parent1::Genotype, parent2::Genotype)
    offspring = cross(parent1, parent2, 10000)
    
    affected = count(o -> o.allele1 == 'a' && o.allele2 == 'a', offspring)
    carriers = count(o -> ('A' in [o.allele1, o.allele2]) && 
                          ('a' in [o.allele1, o.allele2]), offspring)
    unaffected = count(o -> o.allele1 == 'A' && o.allele2 == 'A', offspring)
    
    println("Disease Risk Analysis (10000 simulations):")
    println("  Affected (aa): $(round(100affected/10000, digits=1))%")
    println("  Carriers (Aa): $(round(100carriers/10000, digits=1))%")
    println("  Unaffected (AA): $(round(100unaffected/10000, digits=1))%")
end

# Both parents are carriers
carrier1 = Genotype('A', 'a')
carrier2 = Genotype('A', 'a')
disease_risk(carrier1, carrier2)