Introduction

Welcome to the introduction of data management with ORFik experiment. This vignette will walk you through which data and how you import and export with ORFik, focusing on NGS library import and annotation. ORFik is an R package containing various functions for analysis of RiboSeq, RNASeq, RCP-seq, TCP-seq, Chip-seq and Cage data, we advice you to read ORFikOverview vignette, before starting this one.

Motivation

There exists a myriad of formats in bioinformatics, ORFik focuses on NGS data anlysis and therefor supports many functions to import that data. We will here learn how to use ORFik effectivly in importing data (how to convert to the format you should use). That is the format with optimal information content relative to load speed.

Importing Sequencing reads

Sequencing reads starts with fastq files. In the Annotation & Alignment vignette we walk through how to acquire (download usually) and map these reads into bam files. We here presume that this is done and you have the bam files.

ORFik currently supports these formats:

bam (only format with mismatch/full read sequence information)
wig (supported by many downstream tools)
bigwig (fastest for random access on single gene)
bed5 (minimal bed format)
bed12 (supports spliced regions, similar to gtf, also color coding)
ofst (ORFik fst table format, very fast version of bam, bed or wig, contains chromosome, start, end (optional), strand, cigar (optional), all is saved with fst zstandard run length encoding, also writes/reads using multithreading)
covRle (run length encoding of coverage, split by strand, fastest format for full genome calculations, formats: (.rds and .qs, rds to share easier to others, qs using the qs package which supports better compression and reading/writing with multithreading, through the zstandard compression))
covRleList (a list covRles, split by read lengths, very fast for whole genome coverage, where you need each read length separate, formats save as covRle)

All these formats can be imported using the function fimport

ORFik also supports bam files that are collapsed, the number of duplicates is stored in the header information per read (i.e. >seq1_x150, for first read duplicated 150 times). This is a very effective speed up for short read sequencing like Ribo-seq. bigwig etc, internally supports re-collapsing (different cigar for same read length) since it uses the ‘score’ column of GAlignments object.

Loading files (general)

We advice you always to create an ORFik experiment of your data, to make it simpler to code (see Data management vignette).

But here we will show you how to do direct import of files. To load a file, usually it is sufficient to use fimport.

Bam files

# Use bam file that exists in ORFik package
library(ORFik)
bam_file <- system.file("extdata/Danio_rerio_sample", "ribo-seq.bam", package = "ORFik")
footprints <- fimport(bam_file)
# This is identical to:
footprints <- readBam(bam_file)

Bam files are very cumbersome to read, so we should only do this once, and then convert to something faster (described bellow)

Other formats

Other formats are loaded in the same way

# Use bam file that exists in ORFik package
ofst_file <- system.file("extdata/Danio_rerio_sample/ofst", "ribo-seq.ofst", package = "ORFik")
footprints <- fimport(ofst_file)
# This is identical to:
footprints <- import.ofst(ofst_file)

Exporting to new formats

Since bam files are cumbersome to load, we should convert files to more optimized formats. Which formats to convert and export to, depends on if you have the files loaded already or not.

Files are not preload into R

There are several converters in ORFik, here are some examples:

BAM to OFST (keep cigar information)

ofst_out_dir <- file.path(tempdir(), "ofst/")
convert_bam_to_ofst(NULL, bam_file, ofst_out_dir)
# Find the file again
ofst_file <- list.files(ofst_out_dir, full.names = TRUE)[1]
# Load it
fimport(ofst_file)

## GAlignments object with 16649 alignments and 0 metadata columns:
##           seqnames strand       cigar    qwidth     start       end     width
##              <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
##       [1]    chr23      -         28M        28  17599129  17599156        28
##       [2]    chr23      -         28M        28  17599129  17599156        28
##       [3]    chr23      -         28M        28  17599129  17599156        28
##       [4]    chr23      -         28M        28  17599129  17599156        28
##       [5]    chr23      -         28M        28  17599129  17599156        28
##       ...      ...    ...         ...       ...       ...       ...       ...
##   [16645]     chr8      +         29M        29  24068894  24068922        29
##   [16646]     chr8      +         28M        28  24068907  24068934        28
##   [16647]     chr8      +         30M        30  24068919  24068948        30
##   [16648]     chr8      +         30M        30  24068919  24068948        30
##   [16649]     chr8      +         30M        30  24068939  24068968        30
##               njunc
##           <integer>
##       [1]         0
##       [2]         0
##       [3]         0
##       [4]         0
##       [5]         0
##       ...       ...
##   [16645]         0
##   [16646]         0
##   [16647]         0
##   [16648]         0
##   [16649]         0
##   -------
##   seqinfo: 1133 sequences from an unspecified genome; no seqlengths

BAM to bigwig (do not keep cigar information or read lengths)

Bigwig is a bit special, in that it is very fast (and good compression), but to make it failsafe we need the information about size for all chromosomes. Bioconductor functions are very picky about correct chromosome sizes of GRanges objects etc.

bigwig_out_dir <- file.path(tempdir(), "bigwig/")
convert_to_bigWig(NULL, bam_file, bigwig_out_dir, 
                  seq_info = seqinfo(BamFile(bam_file)))
# Find the file again
bigwig_file <- list.files(bigwig_out_dir, full.names = TRUE) 
# You see we have 2 files here, 1 for forward strand, 1 for reverse
# Load it
fimport(bigwig_file)

## GRanges object with 3104 ranges and 1 metadata column:
##            seqnames     ranges strand |     score
##               <Rle>  <IRanges>  <Rle> | <numeric>
##      [1]         MT    1-16596      + |         0
##      [2]    Zv9_NA1    1-22128      + |         0
##      [3]   Zv9_NA10    1-48062      + |         0
##      [4]  Zv9_NA100    1-26026      + |         0
##      [5] Zv9_NA1000      1-942      + |         0
##      ...        ...        ...    ... .       ...
##   [3100]       chr5 1-75682077      - |         0
##   [3101]       chr6 1-59938731      - |         0
##   [3102]       chr7 1-77276063      - |         0
##   [3103]       chr8 1-56184765      - |         0
##   [3104]       chr9 1-58232459      - |         0
##   -------
##   seqinfo: 1133 sequences from an unspecified genome

Files are preload into R

For preloaded files it is better to just convert it there and then, and not convert through the filepath, because then you have just loaded the file twice.

For details, see ?convertLibs and ?convertToOneBasedRanges()

Random access

The really cool thing with bigwig is that it has very fast random access.

bigwig_file <- list.files(bigwig_out_dir, full.names = TRUE) 
# Let us access a location without loading the full file.
random_point <- GRangesList(GRanges("chr24", 22711508, "-"))
# Getting raw vector (Then you need to know which strand it is on:)
bigwig_for_random_point <- bigwig_file[2] # the reverse strand file
rtracklayer::import.bw(bigwig_for_random_point, as = "NumericList",
                              which = random_point[[1]])

## NumericList of length 1
## [[1]] 4

# 4 reads were there

ORFik also has an automatic wrapper for spliced transcript coordinates: As data.table

dt <- coveragePerTiling(random_point, bigwig_file)
dt # Position is 1, because it is relative to first

##    count genes position
##    <int> <int>    <int>
## 1:     4     1        1

Importing Annotation

Annotation consists of 2 primary files:

The full genome sequence (fasta genome + index, 2 individual files)
The gtf (The location of genes on the genome etc.)

For your own analysis, you can download a matching genome and annotation from Ensembl or RefSeq. This yeast example requires an internet connection and is not run when building the vignette:

library(ORFik)
organism <- "Saccharomyces cerevisiae" # Baker's yeast
# This is where you should usually store you annotation references ->
#output_dir <- file.path(ORFik::config()["ref"], gsub(" ", "_", organism))
output_dir <- tempdir()
annotation <- getGenomeAndAnnotation(
                    organism = organism,
                    output.dir = output_dir,
                    assembly_type = "toplevel"
                    )
genome <- annotation["genome"]
gtf <- annotation["gtf"]

The runnable examples below use the small synthetic reference included in ORFik. Copy it to a temporary directory so that indexing and annotation optimization do not write into the installed package. This reference is for learning the API, not for biological analysis.

reference_dir <- system.file("extdata/references/homo_sapiens", package = "ORFik")
output_dir <- tempfile("orfik-reference-")
dir.create(output_dir)
reference_files <- c("Homo_sapiens_dummy.fasta", "Homo_sapiens_dummy.gtf")
file.copy(file.path(reference_dir, reference_files), output_dir)

## [1] TRUE TRUE

genome <- file.path(output_dir, reference_files[1])
gtf <- file.path(output_dir, reference_files[2])
organism <- "Homo sapiens"

The gtf is a very slow format and ORFik will almost never use this directly. We therefor use a much faster format called TxDb (transcript database) The nice thing with using getGenomeAndAnnotation is that it will do a lot of important fixes for you related to import. It will make the fasta index and txdb, and a lot more optimizations that for large species like human make import time go from minutes to less than a second.

If you don’t want to use getGenomeAndAnnotation, you can do it for your own annotation like this:

# Index fasta genome
indexFa(genome)

## [1] "/tmp/Rtmp0yielu/orfik-reference-19fd79fc84cf/Homo_sapiens_dummy.fasta.fai"

# Make TxDb object for large speedup:
txdb <- makeTxdbFromGenome(gtf, genome, organism, optimize = TRUE, return = TRUE)

## Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
## is not available for this TxDb object

The txdb is saved in the same as gtf with a “.db” extension.

Loading Genome index (fasta index)

# Access a FaFile
fa <- findFa(genome)
# Get chromosome information
seqinfo(findFa(genome))

## Seqinfo object with 6 sequences from an unspecified genome:
##   seqnames seqnames seqlengths isCircular genome
##   1            chr1       1161         NA   <NA>
##   2            chr2       1161         NA   <NA>
##   3            chr3       1161         NA   <NA>
##   4            chr4       1161         NA   <NA>
##   5            chr5       1161         NA   <NA>
##   6            chr6       1161         NA   <NA>

# Load a 10 first bases from chromosome 1
txSeqsFromFa(GRangesList(GRanges("chr1", IRanges(1, 10), "+")), fa, TRUE)

## DNAStringSet object of length 1:
##     width sequence
## [1]    10 AGTACACTCC

# Load a 10 first bases from chromosome 1 and 2.
txSeqsFromFa(GRangesList(GRanges("chr1", IRanges(1, 10), "+"),
                       GRanges("chr2", IRanges(1, 10), "+")), fa, TRUE)

## DNAStringSet object of length 2:
##     width sequence
## [1]    10 AGTACACTCC
## [2]    10 GCAAACTCGC

Loading Gene annotation (Txdb)

ORFik makes it very easy to load specific regions from a txdb. We already have the txdb loaded, but let us load it again

  txdb_path <- paste0(gtf, ".db")
  txdb <- loadTxdb(txdb_path)

Loading Transcript regions

Lets get these regions of the transcripts:

All transcripts
All mRNAs (all transcripts with CDS)
Leaders (5’ UTRs)
CDS (Coding sequences, main ORFs)
Trailers (3’ UTRs)

  loadRegion(txdb, "tx")[1]

## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   306-445      + |         1          NA         1
##   [2]     chr1   446-751      + |         2          NA         2
##   [3]     chr1   752-856      + |         3          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome

  loadRegion(txdb, "mrna")[1]

## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   306-445      + |         1          NA         1
##   [2]     chr1   446-751      + |         2          NA         2
##   [3]     chr1   752-856      + |         3          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome

  loadRegion(txdb, "leaders")[1]

## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   306-445      + |         1          NA         1
##   -------
##   seqinfo: 6 sequences from an unspecified genome

  loadRegion(txdb, "cds")[1]

## GRangesList object of length 1:
## $ENSTTEST10001
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |    cds_id    cds_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   446-751      + |         1        <NA>         2
##   -------
##   seqinfo: 6 sequences from an unspecified genome

  loadRegion(txdb, "trailers")

## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   752-856      + |         3          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10002
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr2   306-445      - |         4          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10003
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr3   752-856      + |         9          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## ...
## <3 more elements>

These are output as GRangesList, which are list elements of GRanges (1 list elements per transcript, which can contain multiple GRanges). If the gene region is not spliced, it has only 1 GRanges object.

Loading Transcript regions (filtering)

Your annotations contains many transcripts, the 2022 human GTFs usually contain around 200K transcripts, at least half of these are from non coding transcripts (they have no CDS). So how to filter out what you do not need?

Loading Transcript regions (filtering by length requirements)

Let’s say you only want mRNAs, which have these properties: - Leaders >= 1nt - CDS >= 150nt - Trailers >= 0nt

We can in ORFik get this with:

  tx_to_keep <- filterTranscripts(txdb, minFiveUTR = 1, minCDS = 150, minThreeUTR = 0)
  loadRegion(txdb, "mrna", names.keep = tx_to_keep)

## GRangesList object of length 6:
## $ENSTTEST10001
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1   306-445      + |         1          NA         1
##   [2]     chr1   446-751      + |         2          NA         2
##   [3]     chr1   752-856      + |         3          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10002
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr2   752-856      - |         6          NA         1
##   [2]     chr2   446-751      - |         5          NA         2
##   [3]     chr2   306-445      - |         4          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10003
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr3   306-445      + |         7          NA         1
##   [2]     chr3   446-751      + |         8          NA         2
##   [3]     chr3   752-856      + |         9          NA         3
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## ...
## <3 more elements>

Increasing the minimum leader length makes the filter more restrictive:

  filterTranscripts(txdb, minFiveUTR = 10, minCDS = 150, minThreeUTR = 0)

## [1] "ENSTTEST10001" "ENSTTEST10002" "ENSTTEST10003" "ENSTTEST10004"
## [5] "ENSTTEST10005" "ENSTTEST10006"

Some references, including yeast annotations, have missing or very short leaders. ORFik can add pseudo 5’ UTRs to the TxDb when needed. These are analysis windows, not experimentally supported transcript boundaries; the source GTF is not changed. The bundled reference already has leaders.

  txdb <- makeTxdbFromGenome(gtf, genome, organism, optimize = TRUE, return = TRUE, 
                             pseudo_5UTRS_if_needed = 100)

## Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
## is not available for this TxDb object

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   105.0   113.8   140.0   128.3   140.0   140.0

  filterTranscripts(txdb, minFiveUTR = 10, minCDS = 150, minThreeUTR = 0)[1:3]

## [1] "ENSTTEST10001" "ENSTTEST10002" "ENSTTEST10003"

For detailed access to single points such as start sites or regions, see the ‘working with transcripts vignette’.

Loading Transcript regions (filtering by canonical isoform)

To get canonical mRNA isoform (primary isoform defined e.g. by Ensembl)

  filterTranscripts(txdb, minFiveUTR = 0, minCDS = 1, minThreeUTR = 0, 
                    longestPerGene = TRUE)[1:3]

## [1] "ENSTTEST10001" "ENSTTEST10002" "ENSTTEST10003"

You here get all canonical isoform mRNAs (cds exists since minCDS >= 1)

Loading uORF annotation

You can also add in a uORF annotation defined by ORFik, like this: uORF discovery searches the annotated or inferred leaders. Interpret results from pseudo leaders with the limitations of those inferred boundaries in mind.

  findUORFs_exp(txdb, findFa(genome), startCodon = "ATG|CTG", save_optimized = TRUE)

## GRangesList object of length 9:
## $ENSTTEST10001_1
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10001     chr1   313-432      + | ENSTTEST10001_1
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10001_2
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10001     chr1   380-388      + | ENSTTEST10001_2
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10002_1
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10002     chr2   824-838      - | ENSTTEST10002_1
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## ...
## <6 more elements>

  loadRegion(txdb, "uorfs") # For later use, output like this

## GRangesList object of length 9:
## $ENSTTEST10001_1
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10001     chr1   313-432      + | ENSTTEST10001_1
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10001_2
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10001     chr1   380-388      + | ENSTTEST10001_2
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## $ENSTTEST10002_1
## GRanges object with 1 range and 1 metadata column:
##                 seqnames    ranges strand |           names
##                    <Rle> <IRanges>  <Rle> |     <character>
##   ENSTTEST10002     chr2   824-838      - | ENSTTEST10002_1
##   -------
##   seqinfo: 6 sequences from an unspecified genome
## 
## ...
## <6 more elements>

Session information

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] GenomicFeatures_1.65.0   AnnotationDbi_1.75.2     Biobase_2.73.2          
##  [4] fstcore_0.10.0           ORFik_1.33.1             GenomicAlignments_1.49.2
##  [7] Rsamtools_2.29.0         Biostrings_2.81.9        XVector_0.53.0          
## [10] GenomicRanges_1.65.4     Seqinfo_1.3.2            IRanges_2.47.5          
## [13] S4Vectors_0.51.10        BiocGenerics_0.59.12     generics_0.1.4          
## [16] BiocStyle_2.41.0        
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.3.0                   bitops_1.1-0               
##   [3] gridExtra_2.3.1             httr2_1.3.0                
##   [5] biomaRt_2.69.3              rlang_1.3.0                
##   [7] magrittr_2.0.5              otel_0.2.0                 
##   [9] matrixStats_1.5.0           compiler_4.6.1             
##  [11] RSQLite_3.53.3              txdbmaker_1.9.0            
##  [13] png_0.1-9                   vctrs_0.7.3                
##  [15] stringr_1.6.0               pkgconfig_2.0.3            
##  [17] crayon_1.5.3                fastmap_1.2.0              
##  [19] dbplyr_2.6.0                biomartr_1.0.7             
##  [21] rmarkdown_2.32              UCSC.utils_1.9.0           
##  [23] bit_4.6.0                   xfun_0.61                  
##  [25] cachem_1.1.0                cigarillo_1.3.1            
##  [27] GenomeInfoDb_1.49.1         jsonlite_2.0.0             
##  [29] progress_1.2.3              blob_1.3.0                 
##  [31] DelayedArray_0.39.6         BiocParallel_1.47.0        
##  [33] parallel_4.6.1              prettyunits_1.2.0          
##  [35] R6_2.6.1                    stringi_1.8.9              
##  [37] bslib_0.12.0                RColorBrewer_1.1-3         
##  [39] rtracklayer_1.73.0          jquerylib_0.1.4            
##  [41] Rcpp_1.1.2                  SummarizedExperiment_1.43.0
##  [43] knitr_1.52                  R.utils_2.13.0             
##  [45] BiocBaseUtils_1.15.1        Matrix_1.7-6               
##  [47] tidyselect_1.2.1            abind_1.4-8                
##  [49] yaml_2.3.12                 stringfish_0.19.2          
##  [51] codetools_0.2-20            curl_8.0.0                 
##  [53] lattice_0.23-1              tibble_3.3.1               
##  [55] withr_3.0.3                 KEGGREST_1.53.6            
##  [57] S7_0.2.2                    evaluate_1.0.5             
##  [59] RcppParallel_6.2.1          BiocFileCache_3.3.0        
##  [61] xml2_1.6.0                  pillar_1.11.1              
##  [63] BiocManager_1.30.27         filelock_1.0.3             
##  [65] MatrixGenerics_1.25.0       RCurl_1.98-1.20            
##  [67] hms_1.1.4                   ggplot2_4.0.3              
##  [69] scales_1.4.0                qs2_0.3.1                  
##  [71] glue_1.8.1                  maketools_1.3.2            
##  [73] tools_4.6.1                 BiocIO_1.23.3              
##  [75] data.table_1.18.6.1         sys_3.4.3                  
##  [77] BSgenome_1.81.1             locfit_1.5-9.12            
##  [79] buildtools_1.0.0            XML_3.99-0.24              
##  [81] cowplot_1.2.0               grid_4.6.1                 
##  [83] GenomeInfoDbData_1.2.15     restfulr_0.0.17            
##  [85] cli_3.6.6                   fst_0.9.8                  
##  [87] S4Arrays_1.13.0             dplyr_1.2.1                
##  [89] gtable_0.3.6                R.methodsS3_1.8.2          
##  [91] DESeq2_1.53.3               sass_0.4.10                
##  [93] digest_0.6.39               SparseArray_1.13.2         
##  [95] rjson_0.2.23                farver_2.1.2               
##  [97] memoise_2.0.1               htmltools_0.5.9            
##  [99] R.oo_1.27.1                 lifecycle_1.0.5            
## [101] httr_1.4.9                  bit64_4.8.6

Importing data