NN methodology - danifilho/Evo2_BASF GitHub Wiki

#!/usr/bin/env python
# prepare_dataset_list.py – requer RAM ≥ 120 GB

from pathlib import Path
import tempfile, subprocess, torch, numpy as np, pandas as pd
from tqdm import tqdm

WIN        = 128
EMB_DIR    = Path("/mnt/gs21/scratch/dasilvaf/evo2/embeddings_149952_full")
CRM_BED    = Path("/mnt/gs21/scratch/dasilvaf/evo2/ReMap/remap2022_crm_macs2_TAIR10_v1_0.bed.gz")
OUT_X      = Path("X.npy")
OUT_Y      = Path("y.npy")
OUT_META   = Path("meta.tsv")

def get_labels(bed_rows):
    """Retorna 0/1 para cada linha BED via bedtools intersect."""
    with tempfile.NamedTemporaryFile("w+", delete=False) as tmp:
        tmp.writelines("\t".join(map(str, r)) + "\n" for r in bed_rows)
        tmp.flush()
        found = subprocess.check_output(
            ["bedtools", "intersect", "-u", "-a", tmp.name, "-b", str(CRM_BED)]
        ).decode().splitlines()
    hits = {l.split("\t")[3] for l in found}
    return [1 if r[3] in hits else 0 for r in bed_rows]

rows, labels, meta = [], [], []

for shard in tqdm(sorted(EMB_DIR.glob("*.pt")), desc="shards"):
    chrom, s, e = shard.stem.rsplit("_", 2)
    if int(e) - int(s) != 149_952:
        continue  # descarta shards curtos

    tensor = torch.load(shard, map_location="cpu").float()
    bed_batch = []
    for idx in range(tensor.shape[0]):
        w_start = int(s) + idx * WIN
        w_end   = w_start + WIN
        uid     = f"{shard.stem}::{idx}"
        bed_batch.append([chrom, w_start, w_end, uid])

    batch_labels = get_labels(bed_batch)

    rows.extend(tensor.numpy())
    labels.extend(batch_labels)
    meta.extend({"id": b[3], "label": lab}
                for b, lab in zip(bed_batch, batch_labels))

np.save(OUT_X, np.asarray(rows, dtype="float32"))
np.save(OUT_Y, np.asarray(labels, dtype="uint8"))
pd.DataFrame(meta).to_csv(OUT_META, sep="\t", index=False)
print("Feito!  X:", OUT_X.stat().st_size/1e9, "GB   N =", len(labels))

salloc -p intel18 -N1 -w skl-141 -c 32 --mem=0 -t 02:00:00 -J prep_evo2

module purge
module load foss/2023a
module load BEDTools/2.31.0-GCC-12.3.0

python -m pip install --quiet --no-cache-dir \
    "numpy<2" pandas tqdm "torch==2.2.2+cpu" \
    -f https://download.pytorch.org/whl/torch_stable.html

O que exatamente fizemos com este script (rotulagem ReMap)?

1. Objetivo geral Criar um dataset (X.npy, y.npy, meta.tsv) onde cada vetor Evo2 (128 bp) recebe rótulo 1 se sobrepõe a um pico ChIP-seq do ReMap 2022 (proxy de CRE), caso contrário 0.

2. Mapeamento de cromossomos (CHR_MAP) Converte IDs NCBI dos shards (NC_003070.9, …) em rótulos do ReMap/TAIR (1–5, C, M). Isso elimina o aviso inconsistent naming convention no bedtools.

3. Contagem global

   • len(SHARDS) → 1 584 shards longos (149 952 bp).
   • Cada shard contém 149 952 embeddings (1 bp). Dividindo por 128 bp → 1 171 janelas.
   • N = 1 854 864 janelas totais; d = 4 096 dimensões.

O que é um shard?

4. Passo 1 – windows.bed

   • Para cada shard cria 1 171 linhas: chrom start end id onde id = ::<índice>
   • 100 % do genoma representado em janelas 128 bp, prontas p/ intersect.

5. Passo 2 – Intersect único

bedtools intersect -u -a windows.bed -b ReMap.bed.gz > hits.bed

-u devolve só as linhas de A que tocam ≥1 pico. - Resultado é lido em HIT_IDS (set de ~1,46 M ids). - Só 1 chamada → segundos em vez de horas.

6. Passo 3 – Memmap de saída

   • X_mm (float32) com shape (N, d) abre em modo “w+” → grava direto no disco, não na RAM.
   • Y_mm (uint8) idem. Memória de trabalho < 4 GB.

7. Passo 4 – Processamento de cada shard

   1. Carrega o tensor original (149 952, d).
   2. Trim: descarta as 64 linhas finais para múltiplo de 128 bp.
   3. Agrupamento.
   4. Gera lista de ids; marca lbl = 1 se id ∈ HIT_IDS.
   5. Copia blocos direto para X_mm e Y_mm; acumula meta.

tens = raw_trim.view(n_win, 128, d).mean(dim=1)  # (1171, d)
# Cada vetor final = média dos 128 embeddings de 1 bp.

8. Passo 5 – Metadados
   • meta.tsv guarda id + label (facilita debug).
   • Impressão final confirma: ✔️ Feito! N = 1 854 864 X ≈ 30 GB.
9. O que obtivemos
   • X.npy 30 GB · (1 854 864 × 4 096) vetores 128 bp.
   • y.npy 1.9 MB · rótulos 0/1 (≈ 79 % positivos).
   • meta.tsv ID + label (47 MB).
   • Processo inteiro ~40 min (32 CPU), sem OOM (memmap).
10. Por que essa abordagem
    • Usa rótulos experimentais (picos ChIP-seq) ⇒ alta confiança.
    • Um único intersect garante rapidez e simplicidade.
    • Memmap evita 120 GB de RAM (lista) e permite escalar.
    • Serve de baseline “biológico” para comparar outras heurísticas (p.ex. janelas com gene).

ReMap 2022 oferece o catálogo mais amplo e controlado de cis‑regulatory modules (CRMs) já reunido para Arabidopsis thaliana (228 k – 233 k regiões). Ele combina curadoria de qualidade ENCODE‑like com clustering de picos ChIP/DAP para entregar hot‑spots de co‑ligação de múltiplos fatores – excelentes rótulos fracos, porém em larga escala, para redes profundas. Trabalhos recentes treinam CNNs, Transformers, GCNs e regressões nesses candidatos e usam os poucos milhares de CREs in vivo validados (VISTA, REDfly, STARR‑seq, EPD‑new) apenas para teste externo.
Essa separação maximiza poder estatístico, evita overfitting e segue o padrão de publicações de 2021‑2025.

• 795 experimentos TAIR10 (464 ChIP‑seq, 331 DAP‑seq) → 4,8 M picos de TF; o clustering gera 228 624 CRMs (TF‑only) ou 233 380 CRMs (TF + histona).
• Arabidopsis dispõe de menos epigenomas públicos que humano; ReMap é um dos poucos bancos que padroniza todas as amostras vegetais em um pipeline único.

• Cada pico passa filtros NSC > 1.05, RSC > 0.8 e FRiP definidos pelo ENCODE.
• CRMs usam somente picos não‑redundantes e exigem ≥ 2 picos por módulo, reduzindo falsos positivos.

• Arquivo BED9 traz (chr, start, end, TF‑list, score = nº picos); basta cortar para BED6 ou converter em TSV/GFF. O score pode entrar como feature numérica.
• Licença CC‑BY 4.0 permite redistribuir o arquivo (citando os autores).

• 10² – 10⁶ exemplos → evita underfitting nos embeddings Evo 2.
• Cobrem múltiplos tecidos, condições e TFs, capturando a diversidade de gramática regulatória.
• Pipeline idêntico facilita expansão para outras espécies.

Muito confiáveis, mas pequenos ou enviesados – treinar só neles leva a superespecialização.

1. Treinar Evo 2 + logistic/MLP sobre ReMap 2022 CRM TF‑only (y = 1) vs. fundo GC‑casado (y = 0).
2. Validar AUROC/AUPRC em VISTA, REDfly, EPD‑new (mantidos fora de treino).
3. Ajustar hyper‑params se o desempenho cair; opcional PU‑learning para refinamento.

• ReMap 2022 – Nucleic Acids Res. 2022
• ANANSE – Nucleic Acids Res. 2021
• C.Origami – Nat. Biotechnol. 2022
• ExplaiNN – Genome Biol. 2023
• CREaTor – Genome Biol. 2023
• (…demais estudos listados na tabela…)