omics.co.in logo omics.co.in

7: Bioconductor Core & GenomicRanges

Tested Stacks: GATK 4.6.1 Samtools 1.21 Nextflow 24.10 Seurat v5.1 Python 3.12
Prerequisites: Basic Linux CLI & Sequencing Formats
📚 Chapter 7 of 13 0% Course Completed
📊 Downloadable Lesson R Script 5.1 KB

r07_genomicranges_s4.R

S4 object models, IRanges/GRanges intervals arithmetic, TSS promoter extraction, finding ChIP-seq peak overlaps, and SummarizedExperiment container overview.

$ Rscript r07_genomicranges_s4.R

Chapter 7: Bioconductor Core: S4 Classes & GenomicRanges

The GenomicRanges package represents the foundational data structure of modern computational genomics. By encapsulating chromosomes, interval coordinates, and strandedness inside validated S4 classes, GRanges provides fast interval arithmetic and overlap queries across millions of features.

1. Constructing GRanges Objects

library(GenomicRanges)

# Define gene models
genes_gr <- GRanges(
    seqnames = c("chr1", "chr1", "chr2", "chr17"),
    ranges   = IRanges(start = c(10000, 50000, 25000, 7668402),
                       end   = c(15000, 58000, 32000, 7687550)),
    strand   = c("+", "-", "+", "-"),
    gene_id  = c("GENE_A", "GENE_B", "GENE_C", "TP53")
)
print(genes_gr)

2. Genomic Arithmetic: Promoters & Flanks

Extracting transcription start site (TSS) promoter regions (e.g. 2000bp upstream, 200bp downstream):

promoters_gr <- promoters(genes_gr, upstream = 2000, downstream = 200)
print(promoters_gr)

3. Finding Interval Overlaps (ChIP-Seq Peaks vs Promoters)

Identifying which transcription factor peaks overlap promoter regions:

overlaps <- findOverlaps(tf_peaks_gr, promoters_gr)

# Extract matching peak-gene pairs
query_idx <- queryHits(overlaps)
subj_idx  <- subjectHits(overlaps)
cat("Overlapping Hits:\n")
print(data.frame(Peak = tf_peaks_gr$peak_id[query_idx], Gene = promoters_gr$gene_id[subj_idx]))
⚠️

Common Pitfalls & Troubleshooting Checklist

  • Java Heap Exhaustion (GATK / Picard): If processes terminate with java.lang.OutOfMemoryError, set your environment variable explicitly: export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp".
  • Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with samtools index -@ 4 sample.bam.
  • Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority: conda config --set channel_priority strict.
  • Chromosome Naming Contamination: Always verify whether reference FASTA headers use chr1 or 1; divergent naming causes silent alignment drops.
Was this guide helpful to your research?
Help us maintain 100% accurate, peer-reviewed bioinformatics standards.
📝 Propose Edit on GitHub ↗