7  Bed tools - Análisis de regiones genómicas

En el capítulo anterior aprendimos a trabajar con un archivo GFF3 de Drosophila melanogaster: eliminamos comentarios, filtramos solo los cromosomas principales y extraímos las anotaciones de genes del cromosoma 2L. Después, convertimos esas anotaciones a formato BED, obteniendo el archivo genes_2L.bed, el cual será utilizado en este capítulo.

En este capítulo aprenderemos a combinar esas regiones anotadas con datos experimentales reales utilizando bedtools, una herramienta diseñada para operar con coordenadas genómicas. Trabajaremos con picos de ChIP-seq para la marca epigenética H3K4me3. Con estos datos podemos responder preguntas como:

Para ello descargaremos un archivo BED de picos desde GEO e integraremos con el archivo que se creó en el capítulo anterior de los genes 2L.

El Gene Expression Omnibus (GEO) es un repositorio público del NCBI dedicado a almacenar, organizar y distribuir datos de genómica funcional. Reúne resultados de técnicas como microarreglos, ChIP-seq, RNA-seq y otras plataformas de alto rendimiento. En GEO se pueden encontrar tanto los datos crudos (por ejemplo, archivos FASTQ o CEL) como los datos procesados y sus archivos suplementarios, que pueden incluir matrices de expresión, archivos BED de picos, anotaciones y metadatos experimentales.

7.1 Descarga los picos de ChiP-seq

En GEO (NCBI) es posible obtener archivos suplementarios que contienen picos ya procesados. Esto evita tener que correr herramientas de peak-calling como MACS2 y nos da un archivo listo para usar.

Para este capítulo usaremos el archivo BED suplementario del experimento GSM853481, que corresponde a ChIP-seq de H3K4me3 en Drosophila. El archivo contiene las coordenadas de picos identificados por MACS.

Ejercicio
  1. Usa el comando wget para descargar el archivo: GSM853481_2011-1333_2011-1334.merged_2011-1335_peaks.bed.gz ubicado en la ruta: https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM853nnn/GSM853481/suppl/

¿Qué comando usarías para descargar este archivo en tu directorio actual?

  1. Descomprime el archivo
  1. wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM853nnn/GSM853481/suppl/GSM853481_2011-1333_2011-1334.merged_2011-1335_peaks.bed.gz
  2. gunzip GSM853481_2011-1333_2011-1334.merged_2011-1335_peaks.bed.gz

7.2 Filtrar el cromosoma 2L

Para mantener coherencia con el archivo de genes preparado en el capítulo anterior, trabajaremos únicamente con los picos del cromosoma 2L.

Ejercicio

Usa el mismo patrón del ejercicio del capítulo anterior para seleccionar solo las líneas que comiencen con 2L seguido de un tabulador.

grep -P '^2L\t' GSM853481_2011-1333_2011-1334.merged_2011-1335_peaks.bed > peaks_2L.bed

7.3 Preparar el archivo de picos para bedtools

Algunos comandos de bedtools requieren que los archivos BED estén ordenados por cromosoma y coordenadas. Para ordenar un archivo, usamos el comando sort del sistema.

sort -k1,1 -k2,2n peaks_2L.bed > peaks_2L.sorted.bed

  • -k1,1: ordena por la primera columna (cromosoma)

  • -k2,2n: luego ordena por la segunda columna (inicio), numéricamente

El archivo peaks_2L.sorted.bed será el que usemos con bedtools.

7.4 Usar bedtools para comparar genes y picos

Hasta este punto, ya tenemos los dos archivos BED que necesitamos:

  • genes_2L.bed — archivo con anotaciones génicas de 2L (obtenido del archivo GFF3 en el capítulo previo).

  • peaks_2L.sorted.bed — generado a partir de los datos de ChIP-seq en 2L

Ahora podemos comenzar a integrar ambas fuentes de datos utilizando bedtools.

7.4.1 Intersect: identificar genes con picos

bedtools intersect permite encontrar regiones que se solapan. Usaremos el archivo de genes como referencia:

bedtools intersect -a genes_2L.bed -b peaks_2L.sorted.bed > genes_peaks_intersect.bed

El archivo resultante incluye las entradas de genes_2L.bed que tienen un pico encima o que se solapan con uno. Esto permite responder preguntas como: ¿Qué genes del cromosoma 2L muestran enriquecimiento de H3K4me3?

7.4.2 Coverage: identificar qué tan cubierto está cada gen

bedtools coverage calcula, para cada región del archivo A, qué proporción está cubierta por regiones del archivo B.

bedtools coverage -a genes_2L.bed -b peaks_2L.sorted.bed > genes_coverage.tsv

La salida incluye información como:

  • cuántos picos caen en cada gen

  • cuántas bases del gen están cubiertas

  • qué fracción del gen está cubierta

Este resultado será la base del ejercicio final en R, donde exploraremos estos valores y generaremos resúmenes visuales a partir de ellos.

7.5 Referencias