Advanced Sequence Analysis

BioLang's advanced sequence analysis functions operate directly on DNA sequences and support IUPAC ambiguity codes throughout. They cover pattern searching with mismatch tolerance, restriction digest simulation, sliding window extraction, GC skew profiling, Shannon entropy scanning, and single-sequence IUPAC matching. All functions follow BioLang's pipe-first design and return structured lists that compose cleanly with filter, map, and to_table.

IUPAC Ambiguity Codes

Several functions accept IUPAC-coded patterns in place of exact nucleotide sequences. The full code table is:

find_pattern

Search for an IUPAC pattern on both strands of a DNA sequence, with optional mismatch tolerance. Each hit reports the position, strand, the actual matched subsequence, and the number of mismatches relative to the pattern.

# find_pattern(seq, pattern, max_mismatches=0)
# Returns List<{pos, strand, matched, mismatches}>
#
# pos     — 0-based start position on the forward strand
# strand  — "+" (forward) or "-" (reverse complement)
# matched — the actual nucleotides at that position
# mismatches — number of positions that differ from the pattern

Exact IUPAC search

# Search for the DnaA box consensus (TTATCCACA) exactly
# The DnaA replication initiator protein recognises this motif near oriC
let chrom = read_fasta("ecoli_K12.fasta") |> first()

let boxes = find_pattern(chrom.seq, "TTATCCACA")
print(len(boxes), "exact DnaA boxes found")
boxes |> map(|h| print(h.strand, h.pos, h.matched))

Search with mismatches

# Allow up to 1 mismatch — typical for biological DnaA box surveys
let boxes = find_pattern(chrom.seq, "TTATCCACA", 1)

# Filter to hits on the forward strand only
let fwd_boxes = boxes |> filter(|h| h.strand == "+")
print(len(fwd_boxes), "forward-strand DnaA boxes (0–1 mm)")

# Inspect mismatching hits
boxes
  |> filter(|h| h.mismatches == 1)
  |> map(|h| {
    pos:      h.pos,
    strand:   h.strand,
    matched:  h.matched,
    mm:       h.mismatches
  })
  |> to_table()
  |> print()

IUPAC degenerate pattern

# Sigma-70 -10 element consensus: TATAAT (W=A|T, so TATAWT is more sensitive)
let promoter_region = dna"GCTTTATAATAGCAAACTTTTCAATCGCTTTATAATAGCGCTTTATAATAGCAAAC"

let hits = find_pattern(promoter_region, "TATAWT")
hits |> map(|h| print("pos:", h.pos, "strand:", h.strand, "seq:", h.matched))

# N wildcard: find all 8-mers starting with GG and ending with CC
let seq = read_fasta("data/genome.fasta") |> first()
find_pattern(seq.seq, "GGNNNNNCC")
  |> to_table()
  |> print()

restriction_sites

Locate restriction enzyme recognition sequences and their cut positions within a DNA sequence. Pass a common enzyme name, "all" to scan for all built-in enzymes at once, or a custom IUPAC pattern string.

# restriction_sites(seq, enzyme_or_pattern)
# Returns List<{pos, end, strand, enzyme, site}>
#
# pos    — 0-based start of the recognition sequence
# end    — 0-based end (exclusive) of the recognition sequence
# strand — "+" or "-"
# enzyme — enzyme name (or "custom" for bare IUPAC patterns)
# site   — the matched recognition sequence in the molecule

Built-in enzyme names

The following common enzymes are recognised by name: EcoRI, BamHI, HindIII, NotI, XhoI, XbaI, SalI, ClaI, SphI, KpnI, SmaI, PstI, SacI, AvaI, NcoI, NheI, MluI, SpeI, ApaI, EcoRV.

In silico single digest

# Find all EcoRI sites in a cloning vector
# EcoRI recognises GAATTC and cuts between G and AATTC
let vector = read_fasta("pUC19.fasta") |> first()

let ecori_sites = restriction_sites(vector.seq, "EcoRI")
print(len(ecori_sites), "EcoRI sites in pUC19")
ecori_sites |> map(|s| print("pos:", s.pos, "site:", s.site))

# Predict fragment sizes after digestion
let positions = ecori_sites |> map(|s| s.pos) |> collect()
let vec_len   = len(vector.seq)

# Append the circular wrap-around position then compute gaps
let cuts = [positions, [positions[0] + vec_len]] |> flatten()
zip(cuts, cuts |> drop(1))
  |> map(|pair| pair[1] - pair[0])
  |> map(|frag_size| print("Fragment:", frag_size, "bp"))

Double digest and multi-enzyme survey

# Simulate an EcoRI + BamHI double digest of a PCR product
let insert = dna"GAATTCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGGATCC"

let ecori = restriction_sites(insert, "EcoRI")
let bamhi  = restriction_sites(insert, "BamHI")

print("EcoRI sites:", len(ecori))
print("BamHI sites:", len(bamhi))

# Survey all built-in enzymes to find unique cutters
let all_sites = restriction_sites(vector.seq, "all")

let enzymes = all_sites |> map(|s| s.enzyme) |> unique()
enzymes
  |> map(|enzyme| {
    enzyme: enzyme,
    sites: all_sites |> filter(|s| s.enzyme == enzyme)
  })
  |> filter(|group| len(group.sites) == 1)
  |> map(|group| {
    enzyme: group.enzyme,
    pos:    group.sites[0].pos,
    site:   group.sites[0].site
  })
  |> to_table()
  |> print()

Custom IUPAC recognition pattern

# Find all sites matching a degenerate recognition sequence
# e.g. a hypothetical enzyme that cuts GRCGYC (R=AG, Y=CT)
let custom_sites = restriction_sites(vector.seq, "GRCGYC")
custom_sites |> map(|s| print(s.pos, s.site))

windows

Extract overlapping or non-overlapping subsequences from a sequence using a sliding window. The default step of 1 produces fully overlapping windows; set step equal to size for non-overlapping tiles.

# windows(seq, size, step=1)
# Returns List<{pos, end, seq, length}>
#
# pos    — 0-based start of the window
# end    — 0-based end (exclusive)
# seq    — the subsequence in that window
# length — window length (equals size for interior windows)

AT-content scan with a sliding window

# Scan a sequence for AT-rich regions using 200 bp windows, 50 bp step
let chrom = read_fasta("ecoli_K12.fasta") |> first()

chrom.seq
  |> windows(200, 50)
  |> map(|w| {
    let bases = base_counts(w.seq)
    {
      pos:    w.pos,
      end:    w.end,
      at_pct: (bases.A + bases.T) / w.length * 100.0
    }
  })
  |> filter(|w| w.at_pct > 70.0)
  |> to_table()
  |> print()

Non-overlapping tiles (genome binning)

# Tile a chromosome into 1 kb non-overlapping bins
let bins = chrom.seq |> windows(1000, 1000)
print(len(bins), "1 kb bins")

# Compute GC content per bin
bins
  |> map(|w| {
    bin:    w.pos / 1000,
    pos:    w.pos,
    gc_pct: gc_content(w.seq) * 100.0
  })
  |> to_table()
  |> write_csv("gc_per_bin.csv")

K-mer frequency in each window

# Count a specific k-mer in every 500 bp window
let seq = read_fasta("data/genome.fasta") |> first()

seq.seq
  |> windows(500, 100)
  |> map(|w| {
    pos:       w.pos,
    cg_count:  motif_count(w.seq, "CG")
  })
  |> to_table()
  |> print()

gc_skew

Compute per-window GC skew and its cumulative sum along a sequence. GC skew is defined as (G - C) / (G + C) for each window. The cumulative sum changes sign at the origin of replication (oriC) and terminus (ter) of bacterial chromosomes, making this a reliable method for replication origin prediction.

# gc_skew(seq, window_size=1000)
# Returns List<{pos, end, g, c, skew, cumulative_skew}>
#
# pos              — window start (0-based)
# end              — window end (exclusive)
# g                — G count in this window
# c                — C count in this window
# skew             — (G - C) / (G + C) for this window
# cumulative_skew  — running sum of skew values up to and including this window

OriC prediction in a bacterial chromosome

# Predict the origin of replication in E. coli K-12
# The minimum of cumulative GC skew marks oriC;
# the maximum marks the terminus of replication.
let chrom = read_fasta("ecoli_K12.fasta") |> first()

let skew = gc_skew(chrom.seq, 10000)

# Find the window with the minimum cumulative skew (oriC)
let oric_window = skew |> sort_by(|w| w.cumulative_skew) |> first()
print("Predicted oriC region:", oric_window.pos, "-", oric_window.end)

# Find the window with the maximum cumulative skew (ter)
let ter_window = skew |> sort_by(|w| w.cumulative_skew) |> last()
print("Predicted ter region: ", ter_window.pos, "-", ter_window.end)

Export skew profile for visualisation

# Write the full GC skew profile to CSV for plotting
gc_skew(chrom.seq, 5000)
  |> map(|w| {
    pos:             w.pos,
    end:             w.end,
    skew:            w.skew,
    cumulative_skew: w.cumulative_skew
  })
  |> to_table()
  |> write_csv("ecoli_gc_skew.csv")

print("Wrote GC skew profile")

Compare multiple chromosomes

# Predict oriC for every sequence in a multi-FASTA file
read_fasta("bacteria_genomes.fasta")
  |> map(|rec| {
    genome: rec.id,
    oric:   (gc_skew(rec.seq, 10000) |> sort_by(|w| w.cumulative_skew) |> first()).pos,
    ter:    (gc_skew(rec.seq, 10000) |> sort_by(|w| w.cumulative_skew) |> last()).pos
  })
  |> to_table()
  |> print()

entropy

Compute Shannon entropy in bits across a sliding window of a sequence. The entropy at each window is calculated as:

# H = -sum(p_i * log2(p_i))   for i in {A, C, G, T}
# where p_i is the frequency of base i in the window.
#
# Maximum entropy is 2.0 bits (all four bases equally frequent).
# Low entropy regions contain biased base composition — repeats, low-complexity, or homopolymers.

# entropy(seq, window=100, step=1)
# Returns List<{pos, end, entropy}>
#
# pos     — 0-based window start
# end     — 0-based window end (exclusive)
# entropy — Shannon entropy in bits (0.0 – 2.0)

Find low-complexity and repeat regions

# Identify candidate repeat/low-complexity regions for masking
# Threshold of 1.5 bits catches most simple repeats and homopolymers
let seq = read_fasta("data/genome.fasta") |> first()

let low_complexity = seq.seq
  |> entropy(100, 10)
  |> filter(|w| w.entropy < 1.5)

print(len(low_complexity), "low-complexity windows (< 1.5 bits)")

low_complexity
  |> map(|w| { pos: w.pos, end: w.end, entropy: w.entropy })
  |> to_table()
  |> write_csv("low_complexity_regions.csv")

Soft-mask a sequence based on entropy

# Collect all low-entropy window coordinates, then report coverage
let windows_to_mask = seq.seq
  |> entropy(100, 1)
  |> filter(|w| w.entropy < 1.2)

let masked_bases = windows_to_mask |> map(|w| w.end - w.pos) |> sum()
let total_bases  = len(seq.seq)
print("Low-complexity coverage:", masked_bases / total_bases * 100.0, "%")

Compare entropy profiles across sequences

# Average entropy per sequence — useful for quality filtering of assemblies
read_fasta("assembly.fasta")
  |> map(|rec| {
    contig:      rec.id,
    length:      len(rec.seq),
    mean_entropy: entropy(rec.seq, 100, 100)
                    |> map(|w| w.entropy)
                    |> mean()
  })
  |> filter(|r| r.mean_entropy < 1.6)    # flag suspicious contigs
  |> to_table()
  |> print()

iupac_match

Test whether a fixed-length sequence matches an IUPAC pattern of the same length. Both the sequence and pattern must have equal length. Returns true if every position in the sequence is consistent with the corresponding IUPAC code in the pattern.

# iupac_match(seq, pattern)
# Returns Bool
#
# seq     — a DNA sequence (no ambiguity codes)
# pattern — an IUPAC pattern of the same length
#
# Sequence and pattern MUST be the same length.

Basic matching

# Check individual sequences against an IUPAC pattern
print(iupac_match(dna"GAATTC", "GAATTC"))   # true  — exact
print(iupac_match(dna"AAATTC", "RAATTC"))   # true  — R matches A or G
print(iupac_match(dna"CAATTC", "RAATTC"))   # false — C is not a purine
print(iupac_match(dna"ATCGAT", "NNNNNN"))   # true  — N matches any base
print(iupac_match(dna"ATCG",   "ATCGAT"))   # false — lengths differ

Filter a list of sequences against a degenerate motif

# Keep only sequences that match the CpG island core motif SSSSCGSSSS
# S = C or G — so this picks GC-rich regions containing a CpG
let candidates = read_fasta("promoter_candidates.fasta") |> collect()

candidates
  |> filter(|r| iupac_match(r.seq, "SSSSCGSSSS"))
  |> map(|r| r.id)
  |> print()

Validate restriction site sequences

# Verify that a library of short oligos all conform to a degenerate design
# Oligos should match the pattern NNRRNNYYNN (degenerate probe design)
let oligos = [
  dna"AAGGAATTAT",
  dna"TTGGCCTTAG",
  dna"CCAACCTTAT",
  dna"GGAAGATTAT",
]

let pattern = "NNRRNNYYNN"

oligos
  |> map(|oligo| {
    seq:   str(oligo),
    valid: iupac_match(oligo, pattern)
  })
  |> to_table()
  |> print()

Full Pipeline Example: Restriction Cloning Site Survey

Combine several advanced functions to characterise a cloning vector, find suitable restriction sites, check insert compatibility, and flag repeat regions to avoid.

# ── 1. Load vector and target insert ────────────────────────────────────────
let vector = read_fasta("pUC19.fasta") |> first()
let insert  = read_fasta("target_gene.fasta") |> first()

# ── 2. Find unique cutters in the Multiple Cloning Site (MCS) ───────────────
let mcs_start = 396
let mcs_end   = 452
let mcs       = vector.seq |> windows(mcs_end - mcs_start, mcs_end - mcs_start)
                            |> filter(|w| w.pos == mcs_start)
                            |> first()

let all_in_vector = restriction_sites(vector.seq, "all")
let unique_cutters = all_in_vector
  |> map(|s| s.enzyme)
  |> unique()
  |> map(|enzyme| {
    enzyme: enzyme,
    sites: all_in_vector |> filter(|s| s.enzyme == enzyme)
  })
  |> filter(|group| len(group.sites) == 1)

print("Unique cutters in pUC19:")
unique_cutters |> map(|group| print(" ", group.enzyme, "at pos", group.sites[0].pos))

# ── 3. Check that selected enzymes do NOT cut the insert ────────────────────
let chosen = ["EcoRI", "HindIII"]

chosen |> map(|enzyme| {
  insert_sites: restriction_sites(insert.seq, enzyme),
  enzyme:       enzyme
})
|> map(|r| print(r.enzyme, "cuts insert", len(r.insert_sites), "time(s)"))

# ── 4. Scan the insert for low-complexity regions before cloning ─────────────
let lc_regions = entropy(insert.seq, 80, 10)
  |> filter(|w| w.entropy < 1.4)

print(len(lc_regions), "low-complexity windows in insert")
lc_regions |> map(|w| print("  pos:", w.pos, "entropy:", w.entropy))

# ── 5. Confirm restriction-site sequences match expected IUPAC patterns ──────
print("Verifying EcoRI site pattern GAATTC:")
restriction_sites(vector.seq, "EcoRI")
  |> map(|s| print(" ", s.site, "valid:", iupac_match(s.site, "GAATTC")))

Full Pipeline Example: Bacterial oriC Prediction + DnaA Box Survey

# ── 1. Load genome ───────────────────────────────────────────────────────────
let genome = read_fasta("mycoplasma_genitalium.fasta") |> first()

# ── 2. Predict oriC by GC skew minimum ──────────────────────────────────────
let skew       = gc_skew(genome.seq, 5000)
let oric_pred  = skew |> sort_by(|w| w.cumulative_skew) |> first()
print("Predicted oriC:", oric_pred.pos, "–", oric_pred.end)

# ── 3. Extract a 50 kb window around the predicted oriC ─────────────────────
let flank    = 25000
let oric_pos = oric_pred.pos
let oric_seq = genome.seq |> windows(50000, 1)
  |> filter(|w| w.pos == (oric_pos - flank))
  |> first()

# ── 4. Find DnaA boxes (with 1 mismatch) near the predicted oriC ─────────────
let dnaa_boxes = find_pattern(oric_seq.seq, "TTATCCACA", 1)
print(len(dnaa_boxes), "DnaA boxes near predicted oriC")
dnaa_boxes
  |> map(|h| {
    genome_pos: oric_pos - flank + h.pos,
    strand:     h.strand,
    matched:    h.matched,
    mismatches: h.mismatches
  })
  |> to_table()
  |> print()

# ── 5. Shannon entropy scan of the oriC window ───────────────────────────────
entropy(oric_seq.seq, 100, 50)
  |> filter(|w| w.entropy < 1.5)
  |> map(|w| {
    genome_pos: oric_pos - flank + w.pos,
    entropy:    w.entropy
  })
  |> to_table()
  |> write_csv("oric_low_complexity.csv")