We used the RNA-Seq data from the publication by Brooks et al. (Brooks et al. 2010). The experiment investigated the effect of siRNA knock-down of pasilla, a gene that is known to bind to mRNA in the spliceosome, and which is thought to be involved in the regulation of splicing. The data set contains 3 biological replicates of the knockdown as well as 4 biological replicates for the untreated control. The data files are available from the NCBI Gene Expression Omnibus under the accession GSE18508. The read sequences in FASTQ format were extracted from the NCBI short read archive file (.sra files) using the SRA toolkit.
The reads in the FASTQ files were aligned using tophat
version 1.2.0 with default parameters against the reference Drosophila
melanogaster genome. The following table summarizes the read number and
alignment statistics. The column exon.counts refers to the
number of reads that could be uniquely aligned to an exon.
dataFilesDir = system.file("extdata", package = "pasilla", mustWork=TRUE)
pasillaSampleAnno = read.csv(file.path(dataFilesDir, "pasilla_sample_annotation.csv"))
pasillaSampleAnno## file condition type number.of.lanes total.number.of.reads
## 1 treated1fb treated single-read 5 35158667
## 2 treated2fb treated paired-end 2 12242535 (x2)
## 3 treated3fb treated paired-end 2 12443664 (x2)
## 4 untreated1fb untreated single-read 2 17812866
## 5 untreated2fb untreated single-read 6 34284521
## 6 untreated3fb untreated paired-end 2 10542625 (x2)
## 7 untreated4fb untreated paired-end 2 12214974 (x2)
## exon.counts
## 1 15679615
## 2 15620018
## 3 12733865
## 4 14924838
## 5 20764558
## 6 10283129
## 7 11653031The reference genome fasta files were obtained from the Ensembl ftp
server. We ran bowtie-build to index the fasta file.
For more information on this procedure see the bowtie
webpage. The indexed form is required by bowtie, and
thus tophat.
wget ftp://ftp.ensembl.org/pub/release-62/fasta/drosophila_melanogaster/ \
dna/Drosophila_melanogaster.BDGP5.25.62.dna_rm.toplevel.fa.gz
gunzip Drosophila_melanogaster.BDGP5.25.62.dna_rm.toplevel.fa.gz
bowtie-build Drosophila_melanogaster.BDGP5.25.62.dna_rm.toplevel.fa \
d_melanogaster_BDGP5.25.62We generated the alignment BAM file using tophat. For
the single-reads data:
For the paired-end data:
tophat -r inner-fragment-size bowtie_index \
<filename>1_1.fastq,<filename>2_1.fastq,...,<filename>N_1.fastq \
<filename>1_2.fastq,<filename>2_2.fastq,...,<filename>N_2.fastqThe SAM alignment files from which pasilla was generated are available at this URL.
To generate the per-exon read counts, we first needed to define the
exonic regions.
To this end, we downloaded the file
Drosophila_melanogaster.BDGP5.25.62.gtf.gz from Ensembl.
The script dexseq_prepare_annotation.py contained in the
DEXSeq
package was used to extract the exons of the transcripts from the file,
define new non-overlapping exonic regions and reformat it to create the
file Dmel.BDGP5.25.62.DEXSeq.chr.gff contained in
pasilla/extdata. For example, for this file we ran:
wget ftp://ftp.ensembl.org/pub/release-62/gtf/ \
drosophila_melanogaster/Drosophila_melanogaster.BDGP5.25.62.gtf.gz
gunzip Drosophila_melanogaster.BDGP5.25.62.gtf.gz
python dexseq_prepare_annotation.py Drosophila_melanogaster.BDGP5.25.62.gtf \
Dmel.BDGP5.25.62.DEXSeq.chr.gffTo count the reads that fell into each non-overlapping exonic part,
the script dexseq_count.py, which is also contained in the
DEXSeq
package, was used. It took the alignment results in the form of a SAM
file (sorted by position in the case of a paired end data) and the
gtf file Dmel.BDGP5.25.62.DEXSeq.chr.gff and
returned one file for each biological replicate with the exon counts.
For example, for the file treated1.bam, which contained single-end
alignments, we ran:
samtools index treated1.bam
samtools view treated1.bam > treated1.sam
python dexseq_count.py Dmel.BDGP5.25.62.DEXSeq.chr.gff \
treated1.sam treated1fb.txtFor the file treated2.bam, which contained paired-end
alignments:
samtools index treated2.bam
samtools view treated2.bam > treated2.sam
sort -k1,1 -k2,2n treated2.sam > treated2_sorted.sam
python dexseq_count.py -p yes Dmel.BDGP5.25.62.DEXSeq.chr.gff \
treated2_sorted.sam treated2fb.txtThe output of the two HTSeq python scripts is provided in the pasilla package:
## [1] "geneIDsinsubset.txt" "treated1fb.txt" "treated2fb.txt"
## [4] "treated3fb.txt" "untreated1fb.txt" "untreated2fb.txt"
## [7] "untreated3fb.txt" "untreated4fb.txt"The Python scripts are built upon the HTSeq library.
To create a DEXSeqDataSet object, we need the count files, the sample annotations, and the GFF file with the per exon annotation.
With these, we could call the function DEXSeqDataSet to
construct the object dxd.
library("DEXSeq")
dxd = DEXSeqDataSetFromHTSeq(
countfiles = file.path(dataFilesDir, paste(pasillaSampleAnno$file, "txt", sep=".")),
sampleData = pasillaSampleAnno,
design= ~ sample + exon + condition:exon,
flattenedfile = gffFile)
dxd ## class: DEXSeqDataSet
## dim: 70463 14
## metadata(1): version
## assays(1): counts
## rownames(70463): FBgn0000003:E001 FBgn0000008:E001 ... FBgn0261575:E001
## FBgn0261575:E002
## rowData names(5): featureID groupID exonBaseMean exonBaseVar
## transcripts
## colnames: NULL
## colData names(8): sample file ... exon.counts exonWe only wanted to work with data from a subset of genes, which was defined in the following file.
genesforsubset = readLines(file.path(dataFilesDir, "geneIDsinsubset.txt"))
dxd = dxd[geneIDs( dxd ) %in% genesforsubset,]We saved the R objects in the data directory of the package:
## 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] DEXSeq_1.58.0 AnnotationDbi_1.74.0
## [3] DESeq2_1.52.0 SummarizedExperiment_1.42.0
## [5] GenomicRanges_1.64.0 Seqinfo_1.2.0
## [7] IRanges_2.46.0 S4Vectors_0.50.3
## [9] MatrixGenerics_1.24.0 matrixStats_1.5.0
## [11] Biobase_2.72.0 BiocGenerics_0.58.1
## [13] generics_0.1.4 BiocParallel_1.46.0
## [15] BiocStyle_2.40.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 dplyr_1.2.1 farver_2.1.2
## [4] blob_1.3.0 filelock_1.0.3 Biostrings_2.80.2
## [7] S7_0.2.2 bitops_1.1-0 fastmap_1.2.0
## [10] BiocFileCache_3.2.0 XML_3.99-0.25 digest_0.6.39
## [13] lifecycle_1.0.5 statmod_1.5.2 survival_3.8-12
## [16] KEGGREST_1.52.2 RSQLite_3.53.3 genefilter_1.94.0
## [19] magrittr_2.0.5 compiler_4.6.1 rlang_1.3.0
## [22] sass_0.4.10 progress_1.2.3 tools_4.6.1
## [25] yaml_2.3.12 knitr_1.52 prettyunits_1.2.0
## [28] S4Arrays_1.12.1 bit_4.6.0 curl_8.0.0
## [31] DelayedArray_0.38.2 RColorBrewer_1.1-3 abind_1.4-8
## [34] hwriter_1.3.2.1 sys_3.4.3 grid_4.6.1
## [37] xtable_1.8-8 ggplot2_4.0.3 scales_1.4.0
## [40] biomaRt_2.68.1 cli_3.6.6 rmarkdown_2.32
## [43] crayon_1.5.3 otel_0.2.0 httr_1.4.9
## [46] DBI_1.3.0 cachem_1.1.0 stringr_1.6.0
## [49] splines_4.6.1 parallel_4.6.1 BiocManager_1.30.27
## [52] XVector_0.52.0 vctrs_0.7.3 Matrix_1.7-6
## [55] jsonlite_2.0.0 geneplotter_1.90.0 hms_1.1.4
## [58] bit64_4.8.6 maketools_1.3.2 locfit_1.5-9.12
## [61] jquerylib_0.1.4 annotate_1.90.0 glue_1.8.1
## [64] codetools_0.2-20 stringi_1.8.9 gtable_0.3.6
## [67] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
## [70] R6_2.6.1 dbplyr_2.6.0 httr2_1.3.0
## [73] evaluate_1.0.5 lattice_0.23-1 png_0.1-9
## [76] Rsamtools_2.28.0 memoise_2.0.1 bslib_0.12.0
## [79] Rcpp_1.1.2 SparseArray_1.12.3 xfun_0.61
## [82] buildtools_1.0.0 pkgconfig_2.0.3