DNAnexus · Genómica en la Nube · dxpy SDK

Pipeline Somático NSCLC — GenomaBio S.L.

🧬 WES Tumor/Normal · hg38 📦 4 Applets · 1 Workflow ⚡ 200+ muestras/mes 🔒 GDPR · HIPAA-ready 💶 ~12€/muestra estimado
4
Applets del pipeline
~45min
Tiempo por muestra
200+
Muestras/mes escalable
mem1_ssd1
_v2_x8
Instancia optimizada coste
FASTQ
→ VCF
End-to-end automatizado
Arquitectura del Pipeline
1
QC & Trimming
FastQC · Trimmomatic
↓ FASTQ R1/R2 (tumor+normal)
↑ FASTQ filtrado · QC report
2
Alineamiento
BWA-MEM2 · samtools
↓ FASTQ trimado
↑ BAM sorted+indexed
3
Variant Calling
GATK Mutect2
↓ BAM tumor · BAM normal
↑ VCF somático filtrado
4
Anotación
VEP · COSMIC · ClinVar
↓ VCF filtrado
↑ VCF anotado · TSV report
5
Informe
Python · pandas
↓ VCF anotado · métricas
↑ HTML report · variantes accionables
📄 dxapp.json
JSON · Configuración del Applet
{
  "name": "somatic-variant-pipeline",
  "title": "Somatic Variant Pipeline NSCLC",
  "summary": "WES tumor/normal pipeline: QC→align→Mutect2→VEP",
  "version": "2.1.0",
  "inputSpec": [
    {
      "name": "tumor_fastq_r1",
      "label": "Tumor FASTQ R1",
      "class": "file",
      "patterns": ["*.fastq.gz"]
    },
    {
      "name": "normal_fastq_r1",
      "label": "Normal FASTQ R1",
      "class": "file",
      "patterns": ["*.fastq.gz"]
    },
    {
      "name": "sample_id",
      "label": "Sample ID",
      "class": "string"
    },
    {
      "name": "reference_genome",
      "label": "Reference (default: hg38)",
      "class": "string",
      "default": "hg38",
      "optional": true
    }
  ],
  "outputSpec": [
    {
      "name": "annotated_vcf",
      "class": "file",
      "label": "Annotated VCF (VEP+COSMIC)"
    },
    {
      "name": "qc_report",
      "class": "file",
      "label": "QC Metrics Report"
    }
  ],
  "runSpec": {
    "interpreter": "python3",
    "file": "src/somatic_pipeline.py",
    "distribution": "Ubuntu",
    "release": "22.04",
    "version": "0",
    "timeoutPolicy": {
      "*": {
        "hours": 6
      }
    },
    "execDepends": [
      {"name": "bwa-mem2"},
      {"name": "samtools"},
      {"name": "fastqc"},
      {
        "name": "trimmomatic",
        "package_manager": "apt"
      }
    ]
  },
  "regionalOptions": {
    "aws:eu-west-2": {
      "systemRequirements": {
        "*": {
          "instanceType": "mem1_ssd1_v2_x8"
        }
      }
    }
  }
}
🐍 src/somatic_pipeline.py
Python · Entrypoint principal
import dxpy
import subprocess
import os
import logging

logging.basicConfig(level=logging.INFO,
    format="%(asctime)s [%(levelname)s] %(message)s")
logger = logging.getLogger("somatic-pipeline")


@dxpy.entry_point('main')
def main(tumor_fastq_r1, tumor_fastq_r2,
          normal_fastq_r1, normal_fastq_r2,
          sample_id, reference_genome="hg38"):

    logger.info(f"Iniciando pipeline para {sample_id}")
    os.makedirs("/genomabio/tmp", exist_ok=True)

    # === PASO 1: Descargar inputs ===
    logger.info("[1/4] Descargando FASTQ tumor + normal")
    for link, path in [
        (tumor_fastq_r1, "/genomabio/tmp/tumor_R1.fastq.gz"),
        (tumor_fastq_r2, "/genomabio/tmp/tumor_R2.fastq.gz"),
        (normal_fastq_r1, "/genomabio/tmp/normal_R1.fastq.gz"),
        (normal_fastq_r2, "/genomabio/tmp/normal_R2.fastq.gz"),
    ]:
        dxpy.download_dxfile(
            link["$dnanexus_link"], path)

    # === PASO 2: QC con FastQC + Trimmomatic ===
    logger.info("[2/4] QC y trimming")
    run_qc(sample_id)
    run_trimmomatic(sample_id)

    # === PASO 3: Alineamiento BWA-MEM2 ===
    logger.info("[3/4] Alineamiento tumor y normal")
    tumor_bam = align_sample(
        "tumor", sample_id, reference_genome)
    normal_bam = align_sample(
        "normal", sample_id, reference_genome)

    # === PASO 4: Mutect2 variant calling ===
    logger.info("[4/4] GATK Mutect2 variant calling")
    raw_vcf = run_mutect2(
        tumor_bam, normal_bam,
        sample_id, reference_genome)
    filtered_vcf = filter_variants(raw_vcf, sample_id)

    # === PASO 5: Anotación VEP ===
    annotated_vcf = run_vep(filtered_vcf, sample_id)

    # === Upload resultados ===
    vcf_file = dxpy.upload_local_file(
        annotated_vcf,
        project=dxpy.PROJECT_CONTEXT_ID,
        folder=f"/results/{sample_id}",
        properties={
            "sample_id": sample_id,
            "pipeline_version": "2.1.0",
            "reference": reference_genome
        },
        tags=["somatic", "annotated", "NSCLC"]
    )
    qc_file = dxpy.upload_local_file(
        f"/genomabio/tmp/{sample_id}_qc.txt",
        folder=f"/results/{sample_id}/qc"
    )

    logger.info(f"Pipeline completado: {sample_id}")
    return {
        "annotated_vcf": dxpy.dxlink(vcf_file),
        "qc_report": dxpy.dxlink(qc_file)
    }


def align_sample(sample_type, sample_id, ref):
    out_bam = f"/genomabio/tmp/{sample_type}.sorted.bam"
    subprocess.check_call([
        "bwa-mem2", "mem", "-t", "8",
        "-R", f"@RG\\tID:{sample_id}\\tSM:{sample_type}\\tPL:ILLUMINA",
        f"/ref/{ref}.fa",
        f"/genomabio/tmp/{sample_type}_R1.trimmed.fastq.gz",
        f"/genomabio/tmp/{sample_type}_R2.trimmed.fastq.gz"
    ], stdout=open(f"/genomabio/tmp/{sample_type}.sam", "w"))
    # Sort + index
    subprocess.check_call([
        "samtools", "sort", "-@", "4", "-o", out_bam,
        f"/genomabio/tmp/{sample_type}.sam"
    ])
    subprocess.check_call(["samtools", "index", out_bam])
    return out_bam


dxpy.run()
🚀 run_batch.py
Python · Orquestación batch — 200+ muestras en paralelo
import dxpy
import json
import sys
from datetime import datetime

"""
run_batch.py — GenomaBio S.L.
Lanza el pipeline somático para todas las muestras de un CSV de manifiesto.
Uso: python run_batch.py manifest.json --project project-genomabio-prod
"""

APPLET_ID = "applet-genomabio-somatic-v210"   # ID del applet desplegado
PROJECT_ID = "project-genomabio-prod"
REF_GENOME = "hg38"


def load_manifest(path: str) -> list:
    """Carga el manifiesto de muestras desde JSON."""
    with open(path) as f:
        return json.load(f)["samples"]


def find_fastq(sample_id: str, sample_type: str, read: str) -> dict:
    """Busca el FASTQ de una muestra en el proyecto DNAnexus."""
    results = list(dxpy.find_data_objects(
        classname="file",
        name=f"{sample_id}_{sample_type}_R{read}*.fastq.gz",
        name_mode="glob",
        project=PROJECT_ID,
        folder="/raw_data",
        return_handler=False
    ))
    if not results:
        raise ValueError(f"FASTQ no encontrado: {sample_id}/{sample_type}/R{read}")
    return {"$dnanexus_link": results[0]["id"]}


def launch_sample(sample: dict, applet) -> dxpy.DXJob:
    """Lanza el pipeline para UNA muestra y devuelve el job."""
    sid = sample["sample_id"]
    inputs = {
        "tumor_fastq_r1":  find_fastq(sid, "tumor", "1"),
        "tumor_fastq_r2":  find_fastq(sid, "tumor", "2"),
        "normal_fastq_r1": find_fastq(sid, "normal", "1"),
        "normal_fastq_r2": find_fastq(sid, "normal", "2"),
        "sample_id":       sid,
        "reference_genome": REF_GENOME
    }
    job = applet.run(
        inputs,
        project=PROJECT_ID,
        folder=f"/results/{sid}",
        name=f"somatic-{sid}-{datetime.now():%Y%m%d}",
        instance_type="mem1_ssd1_v2_x8",
        priority="normal",
        depends_on=[],
        extra_args={"tags": ["batch-run", sample.get("cohort", "NSCLC")]}
    )
    print(f"  ✓ {sid} → {job.get_id()} [launched]")
    return job


def monitor_jobs(jobs: list[dxpy.DXJob]):
    """Espera a todos los jobs y muestra un resumen."""
    print(f"\nMonitorizando {len(jobs)} jobs...\n")
    results = {"done": [], "failed": []}
    for job in jobs:
        try:
            job.wait_on_done(interval=30)
            results["done"].append(job.get_id())
        except dxpy.exceptions.DXJobFailureError as e:
            results["failed"].append({"job": job.get_id(), "error": str(e)})
    return results


def main():
    manifest_path = sys.argv[1] if len(sys.argv) > 1 else "manifest.json"
    samples = load_manifest(manifest_path)
    print(f"GenomaBio S.L. — Batch Somatic Pipeline")
    print(f"Muestras a procesar: {len(samples)}\n")

    dxpy.set_workspace_id(PROJECT_ID)
    applet = dxpy.DXApplet(APPLET_ID)

    jobs = []
    for sample in samples:
        job = launch_sample(sample, applet)
        jobs.append(job)

    results = monitor_jobs(jobs)
    print(f"\n=== Resumen ===")
    print(f"Completados: {len(results['done'])}/{len(jobs)}")
    if results["failed"]:
        print(f"Fallidos: {len(results['failed'])}")
        for f in results["failed"]:
            print(f"  ✗ {f['job']}: {f['error']}")


if __name__ == "__main__":
    main()
🗄️ data_utils.py
Python · Gestión y organización de datos genómicos en DNAnexus
import dxpy
from typing import Iterator

PROJECT_PROD = "project-genomabio-prod"
PROJECT_DEV  = "project-genomabio-dev"


def upload_sequencing_run(local_dir: str, run_id: str,
                           cohort: str = "NSCLC") -> list:
    """
    Sube todos los FASTQs de un run de secuenciación con metadatos.
    Retorna lista de DXFile objects.
    """
    import pathlib
    uploaded = []
    folder = f"/raw_data/{cohort}/{run_id}"

    # Crear estructura de carpetas
    dxpy.api.project_new_folder(
        PROJECT_PROD,
        {"folder": folder, "parents": True}
    )

    for fastq_path in pathlib.Path(local_dir).glob("**/*.fastq.gz"):
        sample_id = fastq_path.stem.split("_")[0]
        dx_file = dxpy.upload_local_file(
            filename=str(fastq_path),
            project=PROJECT_PROD,
            folder=folder,
            properties={
                "run_id": run_id,
                "sample_id": sample_id,
                "cohort": cohort,
                "data_type": "WES"
            },
            tags=["raw", cohort]
        )
        uploaded.append(dx_file)
        print(f"  ↑ {fastq_path.name} → {dx_file.get_id()}")

    return uploaded


def find_samples_by_cohort(cohort: str) -> Iterator[dict]:
    """Busca todas las muestras de un cohort con sus archivos de resultado."""
    yield from dxpy.find_data_objects(
        classname="file",
        name="*.annotated.vcf.gz",
        name_mode="glob",
        properties={"cohort": cohort, "pipeline_version": "2.1.0"},
        project=PROJECT_PROD,
        describe={"fields": {"name": True, "properties": True,
                               "size": True, "created": True}}
    )


def download_results_batch(cohort: str, output_dir: str):
    """Descarga todos los VCFs anotados de un cohort a directorio local."""
    import os
    os.makedirs(output_dir, exist_ok=True)
    count = 0
    for result in find_samples_by_cohort(cohort):
        filename = result["describe"]["name"]
        dxpy.download_dxfile(
            result["id"],
            os.path.join(output_dir, filename),
            project=PROJECT_PROD
        )
        count += 1
    print(f"Descargados {count} VCFs del cohort '{cohort}'")


def archive_old_bams(days_old: int = 90):
    """Archiva BAMs de más de N días para reducir costes de almacenamiento."""
    from datetime import datetime, timedelta
    cutoff = int((datetime.now() - timedelta(days=days_old)).timestamp() * 1000)
    archived = 0
    for obj in dxpy.find_data_objects(
        classname="file", name="*.bam", name_mode="glob",
        project=PROJECT_PROD, created_before=cutoff
    ):
        dxpy.DXFile(obj["id"]).add_tags(["archived"])
        archived += 1
    print(f"Archivados {archived} BAMs > {days_old} días")
    return archived
Manifiesto de Muestras + Stack Tecnológico
📋 manifest.json
JSON · Ejemplo batch real
{
  "project": "project-genomabio-prod",
  "run_date": "2026-06-18",
  "cohort": "NSCLC-fase2",
  "samples": [
    {
      "sample_id": "GB-2026-0041",
      "patient_code": "PAT-0041-T",
      "cohort": "NSCLC",
      "tumor_purity": 0.78,
      "sequencer": "NovaSeq-6000"
    },
    {
      "sample_id": "GB-2026-0042",
      "patient_code": "PAT-0042-T",
      "cohort": "NSCLC",
      "tumor_purity": 0.65,
      "sequencer": "NovaSeq-6000"
    },
    {
      "sample_id": "GB-2026-0043",
      "patient_code": "PAT-0043-T",
      "cohort": "NSCLC",
      "tumor_purity": 0.82,
      "sequencer": "NextSeq-2000"
    }
    // ... 47 muestras más en este lote
  ]
}
📊 Estructura de resultados en DNAnexus
Árbol de carpetas · project-genomabio-prod
project-genomabio-prod/
├── raw_data/
│   └── NSCLC/
│       └── RUN-2026-0618/
│           ├── GB-2026-0041_tumor_R1.fastq.gz
│           ├── GB-2026-0041_tumor_R2.fastq.gz
│           ├── GB-2026-0041_normal_R1.fastq.gz
│           └── GB-2026-0041_normal_R2.fastq.gz
├── results/
│   └── GB-2026-0041/
│       ├── GB-2026-0041.annotated.vcf.gz    ← VEP+COSMIC
│       ├── GB-2026-0041.annotated.vcf.gz.tbi
│       └── qc/
│           ├── GB-2026-0041_fastqc.html
│           └── GB-2026-0041_qc_metrics.txt
├── references/
│   ├── hg38.fa                              ← 3.2 GB
│   ├── hg38.fa.bwt.2bit.64                  ← BWA-MEM2 index
│   ├── gnomAD_v3.1.vcf.gz                   ← filtro población
│   └── cosmic_v97_grch38.vcf.gz             ← anotación COSMIC
└── apps/
    └── somatic-variant-pipeline-v2.1.0/

Total muestras procesadas (junio 2026): 127
Tasa de éxito: 98.4% (2 muestras re-corridas)
Coste medio por muestra: 11.80 €