13  Introducción a BEDtools en R 🧬

13.1 💻

Este documento proporciona una guía básica para trabajar con R, Bedtools, y archivos genómicos.
Incluye:

  • Preparación del ambiente
  • Comandos esenciales de Bedtools
  • Lectura y manipulación de BED/GFF en R
  • Ejemplos prácticos

13.2 📦 Cargar paquetes

Note

{r} # Instalar paquetes si es la primera vez install.packages(c("readr", "dplyr", "ggplot2"))

#Cargar de librerías necesarias:

library(readr) # Lectura de archivos (TSV, CSV, texto) library(dplyr) # Manipulación de datos (pipes, mutate, filter…) library(ggplot2) # Gráficos (histogramas, dispersión, boxplots)

13.2.1 Recordemos lo siguiente:

Note

📌 Descripción de las columnas del resumen de genes

A continuación se explica el significado de cada columna generada en el análisis:

  • chrom
    Cromosoma donde se encuentra el gen (por ejemplo: chr1, chr2, chrX).

  • start
    Posición inicial del gen en el cromosoma (coordenadas genómicas en formato BED).

  • end
    Posición final del gen en el cromosoma.

  • gene_id
    Identificador único del gen (por ejemplo: ENSG00000123456 o el nombre del gen).

  • number_of_peaks
    Número total de picos que intersectan el gen.
    Se refiere a cuántas regiones (por ejemplo, ATAC-seq o ChIP-seq) se solapan con el gen.

  • bases_overlapped
    Cantidad de bases (nucleótidos) dentro del gen que están cubiertas por picos.
    Resume cuánto del gen está siendo “tocado” por regiones de interés.

13.2.2 Lectura de archivos generados por Bedtools

Archivo de coverage

Note

#Leer archivo obtenido con bedtools coverage

coverage <- read.delim(file = "Archivos_Bash/genes_coverage.tsv", col.names = c("chrom", # Cromosoma"start", # Inicio del gen"end", # Fin del gen"gene_id", # ID del gen"number_of_peaks", # Número de picos que intersectan el gen"bases_overlapped", # Número total de bases cubiertas por picos"gene_length", # Longitud total del gen"fraction_covered" # Proporción cubierta por picos))

#Vista general del objeto cargado

glimpse(coverage)head(coverage)

13.2.3 Archivo de intersect (bedtools intersect -wa)

Note

#Leer archivo de intersecciones (solo genes con al menos un pico)intersect_genes <- read.delim("Archivos_Bash/genes_peaks_intersect.bed", col.names = c("chrom", "start", "end", "gene_id"))

#vista general del objeto cargado

glimpse(intersect_genes) head(intersect_genes)

13.2.4 🧪 Chequeos de consistencia

Validar la fracción cubierta

Note

coverage_check <- coverage %>%mutate(fraction_calc = bases_overlapped / gene_length, # Recalcular fraccióndiff_fraction = abs(fraction_calc - fraction_covered) # Comparar con la reportada)

#Si todo está bien, los valores deberían ser cercanos a 0summary(coverage_check$diff_fraction)

13.2.5 Genes con y sin picos (según coverage)

Note

coverage %>% summarise(total_genes = n(), # Número total de genesgenes_con_picos = sum(number_of_peaks > 0), # Genes con ≥1 picogenes_sin_picos = sum(number_of_peaks == 0) # Genes sin picos)

13.2.6 Verificar consistencia del archivo intersect

Note

intersect_genes %>% summarise( genes_intersect_unicos = n_distinct(gene_id), # Cuántos genes aparecen al menos una vezfilas_intersect = n() # Número total de filas (varios picos por gen))

13.2.7 Cruzar coverage con intersect para validar

Añadir columna indicando si un gen aparece en intersect

Note

coverage_validado <- coverage %>%mutate(en_intersect = gene_id %in% intersect_genes$gene_id)

13.2.8 Comparar consistencia entre ambos archivos

Note

coverage_validado %>%count(number_of_peaks > 0, en_intersect) %>%rename(tiene_picos_según_coverage = number_of_peaks > 0)

(Idealmente, genes con picos según coverage deben estar presentes en intersect).

13.2.9 Exploración visual y estadística básica

Note

#Histograma de fracción cubiertaggplot(coverage, aes(x = fraction_covered)) +geom_histogram(bins = 50) +labs( title = "Distribución de la fracción de cobertura por gen",x = "Fracción cubierta",y = "Número de genes")

13.2.10 Resumen de cobertura por grupo (con vs sin picos)

Note

coverage %>%
mutate(tiene_picos = number_of_peaks > 0) %>%
group_by(tiene_picos) %>%
summarise(
n_genes = n(),
mediana_fraccion = median(fraction_covered),
promedio_fraccion = mean(fraction_covered),
max_fraccion = max(fraction_covered)
)

13.2.11 Boxplot de fracción cubierta

Note

coverage %>%mutate(tiene_picos = ifelse(number_of_peaks > 0, "Con picos", "Sin picos")) %>%ggplot(aes(x = tiene_picos, y = fraction_covered)) +geom_boxplot() +labs(title = "Cobertura de picos por gen",x = "Grupo de genes",y = "Fracción cubierta")

13.3 Referencias