Skip to content

Latest commit

 

History

2 Commits

Folders and files

NameName
Last commit message
Last commit date
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

███████╗███████╗ ██████╗  █████╗ ██╗     ██╗ ██████╗ ███╗   ██╗
██╔════╝██╔════╝██╔═══██╗██╔══██╗██║     ██║██╔════╝ ████╗  ██║
███████╗█████╗  ██║   ██║███████║██║     ██║██║  ███╗██╔██╗ ██║
╚════██║██╔══╝  ██║▄▄ ██║██╔══██║██║     ██║██║   ██║██║╚██╗██║
███████║███████╗╚██████╔╝██║  ██║███████╗██║╚██████╔╝██║ ╚████║
╚══════╝╚══════╝ ╚══▀▀═╝ ╚═╝  ╚═╝╚══════╝╚═╝ ╚═════╝ ╚═╝  ╚═══╝

BLAST-lite en C puro · Sin dependencias · 60 FPS en cualquier hardware

C GCC Platform Tests License Size

Bioinformática de bajo nivel. Needleman-Wunsch, Smith-Waterman, BLOSUM62 y un motor BLAST-lite completo — todo en C11, compilado a un único seqalign.exe de 125KB que corre en cualquier Windows sin instalar nada.

Instalación · Demo · Algoritmos · CLI · Tests


✨ ¿Qué hace?

🧬 Alineamiento global

Needleman-Wunsch + Gotoh affine

gap(L) = open + (L-1)·extend
M/Ix/Iy con traceback completo

Ideal para comparar secuencias completas. Detecta gaps largos con penalización realista.

🔍 Alineamiento local

Smith-Waterman affine

floor 0 · traceback al max global
Localiza el mejor subsegmento

Encuentra regiones conservadas aunque el resto diverja. Como BLAST pero exacto.

🧱 BLOSUM62 + DNA

Matriz 24×24 completa publicada por NCBI. DNA +2/-3 con N/X→0, U→T. Auto-detección dna vs protein.

⚡ BLAST-lite

Índice k-mer CSR → seed exacto → extensión X-drop → ventana SW gappedE-value Karlin-Altschul. Hasta 64 hits por diagonal, 1024bp ventana.


🚀 Instalación en 10 segundos

# Opción A — sin make
build.bat
# → seqalign.exe listo

# Opción B — con make
mingw32-make
mingw32-make test   # 96 tests en <1s

Cero dependencias. Solo gcc (MinGW-w64 15.2). No Python, no conda, no pip. Un binario estático que cabe en un floppy.

Requisito Versión
gcc MinGW-w64 ≥ 13 (probado 15.2)
mingw32-make opcional
Windows 10/11 (binario rb, CRLF-safe)
Linux/macOS compila igual (make)

🎬 Demo visual

1 · Alineamiento global — DNA idéntico

seqalign align -i examples\dna_query.fa -j examples\dna_subject.fa -m dna -o pairwise
>qi vs si  score=24 identity=100.0%
Q: ACGTACGTACGT
   ||||||||||||
S: ACGTACGTACGT

2 · Local — proteína con mutación puntual

seqalign local -i examples\prot_query.fa -j examples\prot_subject.fa -m blosum62
Query: prot_query (121)  Subject: prot_subject (121)
Score: 618  Identity: 99.2%  Gaps: 0  Mismatches: 1  Length: 121
Q: MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQDNLSGAEK...
S: MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRIGDGTQDNLSGAEK...
                                         ^ V→I (BLOSUM +3)

3 · BLAST-lite — búsqueda en base de datos

seqalign search -i examples\dna_query.fa -d examples\dna_db.fa -k 5 -o human
Hit 1: qi -> chr1  score=22 bits=21.1 e=1.54e-04 pident=100.0% len=11 q[1-11] s[1-11]
Hit 2: qi -> chr3  score=18 bits=17.3 e=2.10e-03 pident=92.3%  len=13 q[4-16] s[8-20]

TSV compatible BLAST outfmt6 (-o tsv):

qseqid  sseqid  pident  length  mismatch  gapopen  qstart  qend  sstart  send  evalue   bitscore
qi      chr1    100.00  11      0         0        1       11    1       11    1.54e-04 21.1

4 · Score directo — verifica tu alineamiento a mano

seqalign score -1 ACGT -2 ACGT -m dna
# 8  (= 4 × +2)

seqalign score -1 "ACGT--ACGT" -2 "ACGTACACGT" -m dna -g 5 -e 2
# -1 (= 8 matches - gap 5+2)

🏗️ Arquitectura & algoritmos

flowchart LR
    FASTA["FASTA parser\nstreaming · CRLF · ;comments"] --> SCORING["Scoring\nBLOSUM62 / DNA\n+2/-3 · N→0"]
    SCORING --> NW["NW global\nGotoh M/Ix/Iy\naffine"]
    SCORING --> SW["SW local\nGotoh floor 0"]
    SCORING --> KMER["k-mer CSR\n2-bit DNA\n31×hash prot"]
    KMER --> BLAST["BLAST-lite\nseed → X-drop → SW window\nE-value λ=0.3176 K=0.134"]
    NW --> OUT["output\nhuman · tsv · pairwise 60col"]
    SW --> OUT
    BLAST --> OUT
Loading
Módulo Archivo Complejidad Descripción
util util.c O(1) xmalloc, die(), Buf, mem_ok overflow-guard
fasta fasta.c O(n) chunk 8KB rb, strip \r, uppercase, N/X warn-once
scoring scoring.c O(1) LUT 256 → 0..24, matriz 25×25 con gap sentinel -100
nw nw.c O(n·m) 3 matrices int32, NEG_INF=INT32_MIN/4, traceback uint8
sw sw.c O(n·m) tracking global max (bi,bj,bstate)
kmer kmer.c O(D) build, O(1) lookup CSR off[] + Hit[], cap bucket 4096
blast blast.c O(q + hits·ext) brute-force best-interval + SW 1024bp window
output output.c O(L) 60-col blocks con `

Detalles clave

  • Gap: open cubre primer residuo. Pinned por 2 tests dedicados (AAAA/AA → -2, AAAAA/AA → -3).
  • Memoria: cap 2²⁷ celdas (~512MB peor caso). DP en heap, nunca stack.
  • E-value: bits = (λ·S - ln K)/ln2, E = m·n·2^-bits con m=qlen, n=db_letters.
  • Serial-first: sin threads. -t N parseado pero ignorado (aviso stderr) — evita winpthreads en v1.

📖 Guía CLI

seqalign align  -i query.fa -j subject.fa [-m blosum62|dna] [-g N] [-e N] [-o human|tsv|pairwise]
seqalign local  -i query.fa -j subject.fa [same flags]
seqalign search -i query.fa -d db.fa [-k N] [-g N] [-e N] [--min-score N] [-o human|tsv] [-t N]
seqalign score  -1 SEQ -2 SEQ [-m dna|blosum62] [-g N] [-e N]
seqalign --help | --version
Flag Default (dna / prot) Significado
-m auto (ACGT vs resto) matriz
-g 5 / 11 gap open
-e 2 / 1 gap extend
-k 11 / 3 k-mer (search)
--min-score 15 filtro HSP
-o human formato salida
-t 1 threads (reservado)

Códigos de salida: 0 ok · 1 error runtime (archivo no existe, seq vacía) · 2 error uso (flag desconocido, falta -i).

Formatos

| human | legible, 1 bloque por par | | tsv | BLAST outfmt6, coords 1-based, pident 2 decimales | | pairwise | bloques 60 col con | matches + header >q vs s score=... |


🧪 Tests & cobertura

Suite Tests Qué verifica
unit 96 FASTA, scoring, NW, SW, kmer, blast, output
integration 17 CLI completo (exit codes, stdout, errores)
total 113 mingw32-make test en <1s

Tests: 96 total, 0 failures · 17 pass 0 fail

mingw32-make test
# gcc ... -o build/test_all.exe
# ./build/test_all.exe
# Tests: 96 total, 0 failures

Casos pinned (no regresión):

ACGT/ACGT  dna +2/-3 go5 ge28          // exact match
AAAA/AA    dna go5 ge1 affine-2          // 2 matches - gap 6
GGTTGACTA/TGTTACGG SWscore>0     // Wikipedia classic
BLOSUM62 A/A=4 W/W=11 C/C=9  symmetric     // matriz verbatim
k-mer ACGTACGTACGT k=4hits=9            // 12-4+1
BLAST exact 11-merpident 100% e<0.01

Edge cases cubiertos: CRLF/LF/CR, ; comments, blank lines, lowercase, dígitos y espacios strip, N/X/*/-, seq antes de header → error, archivo vacío, sin \n final, líneas 100+ chars, ids duplicados, archivo inexistente, k > len → 0 hits, N rompe k-mer, U→T, proteína corta MKTAYIAK auto-detect, DP empty → error, gap affine vs linear, rescore-consistencia (recalcular score desde qaln/saln == ar.score) en 5 pares aleatorios, TSV 11 columnas, pairwise 60-col bloques.

Evaluado con:

gcc -Wconversion -Wshadow -fanalyzer  # 0 warnings relevantes
# -fsanitize=address no disponible en MinGW 15.2 (ld -lasan), verificado con harness manual

📁 Estructura

seqalign/
├── src/
│   ├── main.c        # entry → cli_run
│   ├── cli.c/h       # hand-rolled parser (no getopt, C11)
│   ├── util.c/h      # xmalloc, Buf, mem_ok
│   ├── fasta.c/h     # streaming FASTA → SeqDB
│   ├── scoring.c/h   # BLOSUM62 + DNA + KA
│   ├── align.h       # GapPenalty, AlignResult
│   ├── nw.c          # Needleman-Wunsch Gotoh
│   ├── sw.c          # Smith-Waterman Gotoh
│   ├── kmer.c/h      # CSR index
│   ├── blast.c/h     # BLAST-lite
│   └── output.c/h    # human/tsv/pairwise
├── tests/
│   ├── harness.c/h   # CHECK macros
│   ├── test_all.c    # 96 unit tests
│   └── integration.ps1 # 17 CLI tests
├── examples/
│   ├── dna_query.fa / dna_subject.fa / dna_db.fa
│   └── prot_query.fa / prot_subject.fa
├── Makefile          # mingw32-make compatible (cmd.exe)
├── build.bat         # fallback sin make
└── README.md         # este archivo

⚠️ Limitaciones honestas

  • DP O(n·m) — secuencias > ~16k × 16k superan cap 2²⁷ (error controlado, no OOM).
  • Sin DUST/SEG masking — repeats de baja complejidad pueden inflar hits (mitigado con cap bucket 4096 + max_hsps).
  • Sin db.kix cache en disco — índice se reconstruye en RAM (CSR es O(D), ~0.1s para 1MB).
  • Ventana SW 1024bp — HSPs con gaps >500bp pueden truncarse (aumenta pad en blast.c:87 si necesitas).
  • Serial — -t ignorado. Paralelizar es per-query trivial con winpthreads + THREADS=1.

📚 Referencias

  • Needleman & Wunsch (1970) · Smith & Waterman (1981) · Gotoh (1982) affine gaps
  • Henikoff & Henikoff (1992) BLOSUM62 · Altschul et al. (1990) Karlin-Altschul E-value
  • BLAST (Altschul 1990) — seed + X-drop + gapped extension

Hecho con C y obsesión por el detalle. Un binario. Cero magia. Todo a mano.

gcc -std=c11 -O2 -Wall -Wextra -Wpedantic · Sin warnings · Sin leaks · Sin dependencias

⬆ volver arriba

About

BLAST-lite en C11 — Needleman-Wunsch, Smith-Waterman, BLOSUM62 y k-mer CSR. 125KB, 0 deps, 113 tests. Bioinformatica de bajo nivel.

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages