Metodologia "contém gene" - danifilho/Evo2_BASF GitHub Wiki

to_std.awk

BEGIN{FS=OFS="\t";
  map["NC_003070.9"]="1"; map["NC_003071.7"]="2"; map["NC_003074.8"]="3";
  map["NC_003075.7"]="4"; map["NC_003076.8"]="5"; map["NC_000932.1"]="C";
  map["NC_037304.1"]="M";  # organelos
  map["Chr1"]="1"; map["Chr2"]="2"; map["Chr3"]="3"; map["Chr4"]="4";
  map["Chr5"]="5"; map["ChrC"]="C"; map["ChrM"]="M";
}
{if($1 in map) $1=map[$1]; print}

awk -f to_std.awk INFILE > OUT_STD.bed

Carregando o modulo

module load BEDTools/2.31.0-GCC-12.3.0

# 3.1  149 952 bp (≈150 kb)
bedtools makewindows -g TAIR10.fna.fai -w 149952 \
  | awk -f to_std.awk > windows149k_std.bed      # 1…5 C M

# 3.2  128 bp
bedtools makewindows -g TAIR10.fna.fai -w 128 \
  | awk -f to_std.awk > tmp128.bed
awk 'BEGIN{FS=OFS="\t"}{print $0,"win128_"NR}' tmp128.bed > windows128_std.bed
rm tmp128.bed

grep -P '\tgene\t' TAIR10.gff3 \
  | awk -f to_std.awk > genes_only_std.gff3

bedtools intersect -u \
  -a windows149k_std.bed \
  -b genes_only_std.gff3 \
  > gene_win.bed

bedtools intersect -u \
  -a windows128_std.bed \
  -b gene_win.bed \
  > hits_gene_embed.bed

Criando o X.npy

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 --no-cache-dir \
    "numpy<2" pandas tqdm "torch==2.2.2+cpu" \
    -f https://download.pytorch.org/whl/torch_stable.html

#!/usr/bin/env python
"""
make_X.py
Cria X.npy (float32 memmap) e meta.tsv usando embeddings Evo2-7B
de shards de 149 952 bp.  RAM (~3 GB) • tempo ~30 min (32 CPU).
"""

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

# ── parâmetros ─────────────────────────────────────────────────────────────
WIN     = 128
EMB_DIR = Path("/mnt/gs21/scratch/dasilvaf/evo2/create_embeddings/embeddings_149952_full")     # ajuste se necessário
OUT_X   = Path("X.npy")
OUT_META= Path("meta.tsv")

# ── descobrir shape global ─────────────────────────────────────────────────
shards = sorted([p for p in EMB_DIR.glob("*.pt")
                 if int(p.stem.rsplit("_",2)[2]) - int(p.stem.rsplit("_",2)[1]) == 149_952])
d = torch.load(shards[0], map_location="cpu").shape[1]
N = len(shards) * (149_952 // WIN)           # 1171 janelas por shard

print(f"{len(shards)} shards  ·  {N:,} janelas  ·  d = {d}")

# ── memmap de saída ────────────────────────────────────────────────────────
X_mm = np.memmap(OUT_X, dtype="float32", mode="w+", shape=(N, d))
meta_rows, ptr = [], 0

# ── processar shards ───────────────────────────────────────────────────────
for sh in tqdm(shards, desc="embeddings"):
    raw = torch.load(sh, map_location="cpu").float()        # (149 952, d)
    n_win = raw.shape[0] // WIN                            # 1171
    tens = raw[: n_win*WIN].view(n_win, WIN, d).mean(dim=1).numpy()
    ids  = [f"{sh.stem}::{i}" for i in range(n_win)]

    X_mm[ptr:ptr+n_win] = tens
    meta_rows.extend({"id": uid} for uid in ids)
    ptr += n_win

# ── salvar metadados ───────────────────────────────────────────────────────
pd.DataFrame(meta_rows).to_csv(OUT_META, sep="\t", index=False)
print("✔️  X.npy e meta.tsv criados")



dasilvaf@skl-141:/mnt/gs21/scratch/dasilvaf/evo2/create_dataset/gene_methodology$ python make_X.py 
1584 shards  ·  1,854,864 janelas  ·  d = 4096
embeddings: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1584/1584 [50:12<00:00,  1.90s/it]
✔️  X.npy e meta.tsv criados

Criando o y_gene.py

#!/usr/bin/env python
# make_y_gene.py  ------------------------------------------------------------
"""
Gera y_gene.npy (rótulo 0/1) usando:
  • windows_std.bed  (janelas 128 bp, IDs shard::idx, cromossomos 1-5 C M)
  • windows149k_std.bed  (janelas 150 kb)
  • genes_only_std.gff3  (genes TAIR padronizados)
Saída: y_gene.npy – vetor uint8 alinhado a X.npy / meta.tsv
"""

from pathlib import Path
import subprocess, numpy as np, pandas as pd

WIN128   = Path("windows_std.bed")
WIN150   = Path("windows149k_std.bed")
GENES    = Path("genes_only_std.gff3")

GENE_WIN = Path("gene_win.bed")           # 150 kb com ≥1 gene
HITS_BED = Path("hits_gene_embed.bed")    # 128 bp positivos
META     = Path("meta.tsv")
OUT_YG   = Path("y_gene.npy")

# 1) janelas 150 kb que contêm gene -----------------------------------------
subprocess.check_call(
    ["bedtools","intersect","-u",
     "-a", str(WIN150), "-b", str(GENES)],
    stdout=GENE_WIN.open("w")
)

# 2) janelas 128 bp dentro dessas 150 kb ------------------------------------
subprocess.check_call(
    ["bedtools","intersect","-u",
     "-a", str(WIN128), "-b", str(GENE_WIN)],
    stdout=HITS_BED.open("w")
)
hit_ids = {l.split("\t")[3].strip() for l in HITS_BED.open()}

# 3) vetor de rótulos --------------------------------------------------------
meta = pd.read_csv(META, sep="\t")
y_gene = np.fromiter((uid in hit_ids for uid in meta.id),
                     dtype="uint8", count=len(meta))
np.save(OUT_YG, y_gene)
print("✔️  y_gene.npy salvo  –  % positivos =", y_gene.mean()*100)

y_gene.npy salvo  –  % positivos = 99.24269380396622

  1. Motivação Se uma janela de 150 kb contém pelo menos um gene, há grande chance de haver algum elemento regulatório ativo (CRE) dentro dela. Resultado: rótulo positivo = “janela 128 bp está dentro de uma janela 150 kb que contém gene”.

  2. Conversão das

bedtools makewindows -g ../evo2_arabidopsis/GCF_000001735.4_TAIR10.1_genomic.fna.fai -w 149952 > windows149k.bed

149 952 bp equivale a 1 171 tokens de 128 bp no Evo2.

  1. Padronização de cromossomos
    • Da mesma forma, criei um awk que converteu IDs NCBI → TAIR10 (Chr1…Chr5 ChrC ChrM).
    • Arquivo resultante: windows149k_chr.bed.
awk 'BEGIN{
        FS = OFS = "\t";                     # TAB como separador
        # ── mapa NCBI → TAIR10 ──────────────────────────
        map["NC_003070.9"] = "Chr1";         # cromossomos 1–5
        map["NC_003071.7"] = "Chr2";
        map["NC_003074.8"] = "Chr3";
        map["NC_003075.7"] = "Chr4";
        map["NC_003076.8"] = "Chr5";
        map["NC_000932.1"] = "ChrC";         # cloroplasto
        map["NC_037304.1"] = "ChrM";         # mitocôndria
     }
     {
        if ($1 in map) $1 = map[$1];         # converte se estiver no mapa
        print                                # devolve a linha TAB-delimitada
     }'  windows149k.bed  >  windows149k_chr.bed
  1. Seleção de janelas que contêm genes
grep -P '\tgene\t' TAIR10_GFF3_genes.gff > genes_only.gff3
bedtools intersect -u -a windows149k_chr.bed -b genes_only.gff3 > gene_win.bed

-	gene_win.bed = subconjunto de janelas 150 kb com ≥1 gene
  1. Janelas pequenas de 128 bp
    • Geramos de novo o BED de 128 bp:
bedtools makewindows -g TAIR10.fna.fai -w 128 > windows128.bed
  • Mesma conversão de cromossomos + coluna-ID → windows128_chr.bed.
  1. Propagação do rótulo
bedtools intersect -u \
    -a windows128_chr.bed \
    -b gene_win.bed \
    > hits128_gene.bed
  • hits128_gene.bed contém só as linhas 128 bp que caem em janelas 150 kb com gene.
  • IDs positivos ≈ 90 % das janelas (genoma de Arabidopsis é gene-denso).
  1. Checagens
    • Conferimos contigs no .fai: 5 nucleares + cloroplasto (NC_000932.1) + mitocôndria (NC_037304.1).
    • Verificamos que todos os BEDs usam Chr1…Chr5 ChrC ChrM.
    • cut -f4 hits128_gene.bed | sort | uniq -d → nenhum ID duplicado.

Eu tive problemas para gerar o y_genes.npy por conta do arquivo windows.bed, que está com os cromossomos identificados por 1, 2, 3, 4, 5, C e M, ao invés da anotação que estou usando. Agora, cirei outro arquivo hits_gene_embed.bed com a anotação correta com

# usa o BED de embeddings (windows.bed) vs. janelas-gene
bedtools intersect -u \
    -a windows.bed \
    -b gene_win.bed \
    > hits_gene_embed.bed