2.1 Example data

To demonstrate the use of DEXSeq, we use the pasilla dataset, an RNA-Seq dataset generated by (Brooks et al. 2011). They investigated the effect of siRNA knock-down of the gene pasilla on the transcriptome of fly S2-DRSC cells. The RNA-binding protein pasilla protein is thought to be involved in the regulation of splicing. (Its mammalian orthologs, NOVA1 and NOVA2, are well-studied examples of splicing factors.) Brooks et al. prepared seven cell cultures, treated three with siRNA to knock-down pasilla and left four as untreated controls, and performed RNA-Seq on all samples. They deposited the raw sequencing reads with the NCBI Gene Expression Omnibus (GEO) under the accession number GSE18508.

##Executability of the code

Usually, Bioconductor vignettes contain automatically executable code, i.e., you can follow the vignette by directly running the code shown, using functionality and data provided with the package. However, it would not be practical to include the voluminous raw data of the pasilla experiment here. Therefore, the code in this section is not automatically executable. You may download the raw data yourself from GEO, as well as the required extra tools, and follow the work flow shown here and in the pasilla vignette (Reyes 2013). From Section 4 on, code is directly executable, as usual. Therefore, we recommend that you just read this section, and try following our analysis in R only from the next section onwards. Once you work with your own data, you will want to come back and adapt the workflow shown here to your data.

2.2 Alignment

The first step of the analysis is to align the reads to a reference genome. If you are using classic alignment methods (i.e. not pseudo-alignment approaches) it is important to align them to the genome, not to the transcriptome, and to use a splice-aware aligner (i.e., a short-read alignment tool that can deal with reads that span across introns) such as TopHat2 (Kim et al. 2013), GSNAP (Wu and Nacu 2010), or STAR (Dobin et al. 2013). The explanation of the analysis workflow presented here starts with the aligned reads in the SAM format. If you are unfamiliar with the process of aligning reads to obtain SAM files, you can find a summary how we proceeded in preparing the pasilla data in the vignette for the pasilla data package (Reyes 2013) and a more extensive explanation, using the same data set, in our protocol article on differential expression calling (Anders et al. 2013).

2.3 HTSeq

The initial steps of a DEXSeq analysis are done using two Python scripts that we provide with DEXSeq. Importantly, these preprocessing steps can also be done using tools equivalent to these Python scripts, for example, using GenomicRanges infrastructure (10.2) or Rsubread (10.3). The following two steps describe how to do this steps using the Python scripts that we provide within DEXSeq.

You do not need to know how to use Python; however you have to install the Python package HTSeq, following the explanations given on the HTSeq web page.

Once you have installed HTSeq, you can use the two Python scripts, dexseq_prepare_annotation.py (described in Section 2.4) and dexseq_count.py (see Section 3), that come with the DEXSeq package. If you have trouble finding them, start R and ask for the installation directory with:

pythonScriptsDir = system.file( "python_scripts", package="DEXSeq" )
list.files(pythonScriptsDir)

## [1] "dexseq_count.py"              "dexseq_prepare_annotation.py"

system.file( "python_scripts", package="DEXSeq", mustWork=TRUE )

## [1] "C:/Users/biocbuild/bbs-3.9-bioc/tmpdir/RtmpYx55ZJ/Rinst1eb030e55ae1/DEXSeq/python_scripts"

The displayed path should contain the two files. If it does not, try to re-install DEXSeq (as usual, with the function BiocManager::install()). An alternative work flow, which replaces the two Python-based steps with R based code, is also available in the Appendix section of this vignette.

2.4 Preparing the annotation

The Python scripts expect a GTF file with gene models for your species. We have tested our tools chiefly with GTF files from Ensembl and hence recommend to use these, as files from other providers sometimes do not adhere fully to the GTF standard and cause the preprocessing to fail. Ensembl GTF files can be found in the “FTP Download” sections of the Ensembl web sites (i.e., Ensembl, EnsemblPlants, EnsemblFungi, etc.). Make sure that your GTF file uses a coordinate system that matches the reference genome that you have used for aligning your reads. The safest way to ensure this is to download the reference genome from Ensembl, too. If you cannot use an Ensembl GTF file, see Appendix ?? for advice on converting GFF files from other sources to make them suitable as input for the dexseq_prepare_annotation.py script.

In a GTF file, many exons appear multiple times, once for each transcript that contains them. We need to “collapse” this information to define exon counting bins, i.e., a list of intervals, each corresponding to one exon or part of an exon. Counting bins for parts of exons arise when an exonic region appears with different boundaries in different transcripts. See Figure 1 of the DEXSeq paper (Anders, Reyes, and Huber (2012)) for an illustration. The Python script dexseq_prepare_annotation.py takes an Ensembl GTF file and translates it into a GFF file with collapsed exon counting bins.

Make sure that your current working directory contains the GTF file and call the script (from the command line shell, not from within R) with

python /path/to/library/DEXSeq/python_scripts/dexseq_prepare_annotation.py Drosophila_melanogaster.BDGP5.72.gtf Dmel.BDGP5.25.62.DEXSeq.chr.gff

In this command, which should be entered as a single line, replace /path/to…/python_scripts with the correct path to the Python scripts, which you have found with the call to system.file() shown above. Drosophila_melanogaster.BDGP5.72.gtf is the Ensembl GTF file (here the one for fruit fly, already de-compressed) and Dmel.BDGP5.25.62.DEXSeq.chr.gff is the name of the output file.

In the process of forming the counting bins, the script might come across overlapping genes. If two genes on the same strand are found with an exon of the first gene overlapping with an exon of the second gene, the script’s default behaviour is to combine the genes into a single “aggregate gene” which is subsequently referred to with the IDs of the individual genes, joined by a plus (‘+’) sign. If you do not like this behaviour, you can disable aggregation with the option -r no. Without aggregation, exons that overlap with other exons from different genes are simply skipped.

3.1 Reading the data in to R

The remainder of the analysis is now done in R. We will use the output of the python scripts for the pasilla experiment, that can be found in the package pasilla. Open an R session and type:

inDir = system.file("extdata", package="pasilla")
countFiles = list.files(inDir, pattern="fb.txt$", full.names=TRUE)
basename(countFiles)

## [1] "treated1fb.txt"   "treated2fb.txt"   "treated3fb.txt"  
## [4] "untreated1fb.txt" "untreated2fb.txt" "untreated3fb.txt"
## [7] "untreated4fb.txt"

flattenedFile = list.files(inDir, pattern="gff$", full.names=TRUE)
basename(flattenedFile)

## [1] "Dmel.BDGP5.25.62.DEXSeq.chr.gff"

Now, we need to prepare a sample table. This table should contain one row for each library, and columns for all relevant information such as name of the file with the read counts, experimental conditions, technical information and further covariates. To keep this vignette simple, we construct the table on the fly.

sampleTable = data.frame(
   row.names = c( "treated1", "treated2", "treated3", 
      "untreated1", "untreated2", "untreated3", "untreated4" ),
   condition = c("knockdown", "knockdown", "knockdown",  
      "control", "control", "control", "control" ),
   libType = c( "single-end", "paired-end", "paired-end", 
      "single-end", "single-end", "paired-end", "paired-end" ) )

While this is a simple way to prepare the table, it may be less error-prone and more prudent to used an existing table that had already been prepared when the experiments were done, save it in CSV format and use the R function read.csv() to load it.

In any case, it is vital to check the table carefully for correctness.

sampleTable

##            condition    libType
## treated1   knockdown single-end
## treated2   knockdown paired-end
## treated3   knockdown paired-end
## untreated1   control single-end
## untreated2   control single-end
## untreated3   control paired-end
## untreated4   control paired-end

Our table contains the sample names as row names and the two covariates that vary between samples: first the experimental condition (factor condition with levels control and treatment) and the library type (factor libType), which we included because the samples in this particular experiment were sequenced partly in single-end runs and partly in paired-end runs. For now, we will ignore this latter piece of information, and postpone the discussion of how to include such additional covariates to Section 6. If you have only a single covariate and want to perform a simple analysis, the column with this covariate should be named condition.

Now, we construct an DEXSeqDataSet object from this data. This object holds all the input data and will be passed along the stages of the subsequent analysis. We construct the object with the function DEXSeqDataSetFromHTSeq(), as follows:

library( "DEXSeq" )

dxd = DEXSeqDataSetFromHTSeq(
   countFiles,
   sampleData=sampleTable,
   design= ~ sample + exon + condition:exon,
   flattenedfile=flattenedFile )

The function takes four arguments. First, a vector with names of count files, i.e., of files that have been generated with the dexseq_count.py script. The function will read these files and arrange the count data in a matrix, which is stored in the DEXSeqDataSet object dxd. The second argument is our sample table, with one row for each of the files listed in the first argument. This information is simply stored as is in the object’s meta-data slot (see below). The third argument is a formula of the form ~ sample + exon + condition:exon that specifies the contrast with of a variable from the sample table columns and the ‘exon’ variable. Using this formula, we are interested in differences in exon usage due to the “condition” variable changes. Later in this vignette, we will how to add additional variables for complex designs. The fourth argument is a file name, now of the flattened GFF file that was generated with dexseq_prepare_annotation.py and used as input to dexseq_count.py when creating the count file.

There are other ways to get a DEXSeq analysis started. See Appendix 10.2 and Ref. (Love (2013)) for details.

4.1 Loading and inspecting the example data

To demonstrate the this work flow, we will use the DEXSeqDataSet constructed in the previous section. However, in order to keep the run-time of this vignette small, we will subset the object to only a few genes.

genesForSubset = read.table( 
  file.path(inDir, "geneIDsinsubset.txt"), 
  stringsAsFactors=FALSE)[[1]]

dxd = dxd[geneIDs( dxd ) %in% genesForSubset,]

The DEXSeqDataSet class is derived from the DEXSeqDataSet. As such, it contains the usual accessor functions for the column data, row data, and some specific ones. The core data in an DEXSeqDataSet object are the counts per exon. Each row of the DEXSeqDataSet contains in each column the count data from a given exon (‘this’) as well as the count data from the sum of the other exons belonging to the same gene (‘others’). This annotation, as well as all the information regarding each column of the DEXSeqDataSet, is specified in the colData().

colData(dxd)

## DataFrame with 14 rows and 4 columns
##         sample condition    libType     exon
##       <factor>  <factor>   <factor> <factor>
## 1     treated1 knockdown single-end     this
## 2     treated2 knockdown paired-end     this
## 3     treated3 knockdown paired-end     this
## 4   untreated1   control single-end     this
## 5   untreated2   control single-end     this
## ...        ...       ...        ...      ...
## 10    treated3 knockdown paired-end   others
## 11  untreated1   control single-end   others
## 12  untreated2   control single-end   others
## 13  untreated3   control paired-end   others
## 14  untreated4   control paired-end   others

We can access the first 5 rows from the count data by doing,

head( counts(dxd), 5 )

##                  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## FBgn0000256:E001   92   28   43   54  131   51   49 1390  829   923  1115
## FBgn0000256:E002  124   80   91   76  224   82   95 1358  777   875  1093
## FBgn0000256:E003  340  241  262  347  670  260  297 1142  616   704   822
## FBgn0000256:E004  250  189  201  219  507  242  250 1232  668   765   950
## FBgn0000256:E005   96   38   39   71   76   57   62 1386  819   927  1098
##                  [,12] [,13] [,14]
## FBgn0000256:E001  2495  1054  1052
## FBgn0000256:E002  2402  1023  1006
## FBgn0000256:E003  1956   845   804
## FBgn0000256:E004  2119   863   851
## FBgn0000256:E005  2550  1048  1039

Note that the number of columns is 14, the first seven (we have seven samples) corresponding to the number of reads mapping to out exonic regions and the last seven correspond to the sum of the counts mapping to the rest of the exons from the same gene on each sample.

split( seq_len(ncol(dxd)), colData(dxd)$exon )

## $this
## [1] 1 2 3 4 5 6 7
## 
## $others
## [1]  8  9 10 11 12 13 14

We can also access only the first five rows from the count belonging to the exonic regions (‘this’) (without showing the sum of counts from the rest of the exons from the same gene) by doing,

head( featureCounts(dxd), 5 )

##                  treated1 treated2 treated3 untreated1 untreated2 untreated3
## FBgn0000256:E001       92       28       43         54        131         51
## FBgn0000256:E002      124       80       91         76        224         82
## FBgn0000256:E003      340      241      262        347        670        260
## FBgn0000256:E004      250      189      201        219        507        242
## FBgn0000256:E005       96       38       39         71         76         57
##                  untreated4
## FBgn0000256:E001         49
## FBgn0000256:E002         95
## FBgn0000256:E003        297
## FBgn0000256:E004        250
## FBgn0000256:E005         62

In both cases, the rows are labelled with gene IDs (here Flybase IDs), followed by a colon and the counting bin number. Since a counting bin corresponds to an exon or part of an exon, this ID is called the feature ID or exon ID within DEXSeq. The table content indicates the number of reads that have been mapped to each counting bin in the respective sample.

To see details on the counting bins, we also print the first 3 lines of the feature data annotation:

head( rowRanges(dxd), 3 )

## GRanges object with 3 ranges and 5 metadata columns:
##                    seqnames          ranges strand |   featureID     groupID
##                       <Rle>       <IRanges>  <Rle> | <character> <character>
##   FBgn0000256:E001    chr2L 3872658-3872947      - |        E001 FBgn0000256
##   FBgn0000256:E002    chr2L 3873019-3873322      - |        E002 FBgn0000256
##   FBgn0000256:E003    chr2L 3873385-3874395      - |        E003 FBgn0000256
##                        exonBaseMean      exonBaseVar
##                           <numeric>        <numeric>
##   FBgn0000256:E001               64 1250.66666666667
##   FBgn0000256:E002 110.285714285714 2769.57142857143
##   FBgn0000256:E003 345.285714285714 22147.9047619048
##                                                                                                                                                 transcripts
##                                                                                                                                                      <list>
##   FBgn0000256:E001 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
##   FBgn0000256:E002 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
##   FBgn0000256:E003 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
##   -------
##   seqinfo: 14 sequences from an unspecified genome; no seqlengths

So far, this table contains information on the annotation data, such as gene and exon IDs, genomic coordinates of the exon, and the list of transcripts that contain an exon.

The accessor function annotationData() shows the design table with the sample annotation (which was passed as the second argument to DEXSeqDataSetFromHTSeq()):

sampleAnnotation( dxd )

## DataFrame with 7 rows and 3 columns
##       sample condition    libType
##     <factor>  <factor>   <factor>
## 1   treated1 knockdown single-end
## 2   treated2 knockdown paired-end
## 3   treated3 knockdown paired-end
## 4 untreated1   control single-end
## 5 untreated2   control single-end
## 6 untreated3   control paired-end
## 7 untreated4   control paired-end

In the following sections, we will update the object by calling a number of analysis functions, always using the idiom dxd = someFunction( dxd ), which takes the dxd object, fills in the results of the performed computation and writes the returned and updated object back into the variable dxd.

4.2 Normalisation

DEXSeq uses the same method as DESeq and DESeq2, which is provided in the function estimateSizeFactors().

dxd = estimateSizeFactors( dxd )

4.3 Dispersion estimation

To test for differential exon usage, we need to estimate the variability of the data. This is necessary to be able to distinguish technical and biological variation (noise) from real effects on exon usage due to the different conditions. The information on the strength of the noise is inferred from the biological replicates in the data set and characterized by the so-called dispersion. In RNA-Seq experiments the number of replicates is typically too small to reliably estimate variance or dispersion parameters individually exon by exon, and therefore, variance information is shared across exons and genes, in an intensity-dependent manner.

In this section, we discuss simple one-way designs: In this setting, samples with the same experimental condition, as indicated in the condition factor of the design table (see above), are considered as replicates – and therefore, the design table needs to contain a column with the name condition. In Section 6, we discuss how to treat more complicated experimental designs which are not accommodated by a single condition factor.

To estimate the dispersion estimates, DEXSeq uses the approach of the package DESeq2. Internally, the functions from DEXSeq are called, adapting the parameters of the functions for the specific case of the DEXSeq model. Briefly, per-exon dispersions are calculated using a Cox-Reid adjusted profile likelihood estimation, then a dispersion-mean relation is fitted to this individual dispersion values andfinally, the fitted values are taken as a prior in order to shrink the per-exon estimates towards the fitted values. See the DESeq2 paper for the rational behind the shrinkage approach (Love, Huber, and Anders (2014)).

dxd = estimateDispersions( dxd )

As a shrinkage diagnostic, the DEXSeqDataSet inherits the method plotDispEsts() that plots the per-exon dispersion estimates versus the mean normalised count, the resulting fitted valuesand the a posteriori (shrinked) dispersion estimates (Figure 1).

plotDispEsts( dxd )

Figure 1: Fit Diagnostics
The initial per-exon dispersion estimates (shown by black points), the fitted mean-dispersion values function (red line), and the shrinked values in blue

10.1 Controlling FDR at the gene level

DEXSeq returns one \(p\)-value for each exonic region. Using these p-values, a rejection threshold is established according to a multiple testing correction method. Consider a setting with \(M\) exonic regions, where \(p_{k}\) is the \(p\)-value for the \(k\) exonic region and \(\theta \in {p_{1}, ..., p_{M}} \subset [0, 1]\) is a rejection threshold, leading to \(V=\left|\left\{k: p_k\le\theta\right\}\right|\) rejections (i.e. number of exonic regions being DEU). The Benjamini-Hochberg method establishes that FDR is controlled at the level \(q^{*}\) where

\(q^* = \frac{M \theta}{\left|\left\{k: p_k\le\theta\right\}\right|}\).

Note that \(\theta\) in the denominator of this expression is the probability of rejecting a single hypothesis if it is true. \(q^*\) allows to estimate the number of exons that are differentially used at a specific FDR.

However, sometimes it is informative to know what is the FDR at the gene level, i.e. knowing the numer of genes with at least one differentially used exon while keeping control of FDR. For this scenario, let \(p_{i,l}\) be the \(p\)-value from the DEU test for exon \(l\) in gene \(i\), let \(n_i\) the number of exons in gene \(i\), and \(M\) the number of genes. A gene \(i\) has at least one DEU exon if \(\exists l\) such that \(p_{i,l}\le\theta\). Then, FDR is controlled at the level

\(\frac{\sum_{i=1}^{M} 1-(1-\theta)^{n_i}}{\left|\left\{i: \exists l:p_{il}\le\theta\right\}\right|}\).

The implementation of the formula above is provided by the function perGeneQValue(). For the pasilla example, the code below calculates the number of genes (at a FDR of 10%) with at least one differentially used exon.

numbOfGenes <- sum( perGeneQValue(dxr) < 0.1)
numbOfGenes

## [1] 9

10.2 Preprocessing within R

As an alternative to the approach described in the introductory section, users can also create DEXSeqDataSeq objects from other R data objects. The code for implementationg these functions was kindly contributed by Michael I. Love. For details, see the parathyroidSE package vignette (Love 2013). The work flow is similar to the one using the python scripts.

Note: The code in this section is not run when the vignette is built, as some of the commands have long run time. Therefore, no output is given.

We use functionality from the following Bioconductor packages

library(GenomicRanges)
library(GenomicFeatures)
library(GenomicAlignments)

We demonstrate the workflow briefly (for more details, see (Love 2013)) on the data set by Haglund et al. (2012), which is provided as example data in the parathyroidSE data package.

First, we download the current human gene model annotation from Ensembl via Biomart and create a transcript data base from these. Note that this step takes some time.

hse = makeTxDbFromBiomart( biomart="ensembl", dataset="hsapiens_gene_ensembl", host="grch37.ensembl.org" )

Next, we collapse the gene models into counting bins.

exonicParts = exonicParts( hse, linked.to.single.gene.only = TRUE )

As before, we have to choose how to handle genes with overlapping exons. The linked.to.single.gene.only = option here plays the same role as the -r option to dexseq_prepare_anotation.py. The exonicParts object contains a GRanges object with our counting bins. We use it to count the number of read fragments that overlap with the bins by means of the function summarizeOverlaps(). To demonstrate this, we first determine the paths to the example BAM files in the parathyroidSE data package.

bamDir = system.file( "extdata", package="parathyroidSE", mustWork=TRUE )
fls = list.files( bamDir, pattern="bam$", full=TRUE )

Then, use the following code to count the reads overlapping the bins.

bamlst = BamFileList( fls, index=character(), yieldSize=100000, obeyQname=TRUE )
SE = summarizeOverlaps( exonicParts, bamlst, mode="Union", singleEnd=FALSE, ignore.strand=TRUE, inter.feature=FALSE, fragments=TRUE )

We can now call the function DEXSeqDataSetFromSE() to build an DEXSeqDataSet object. We modify the colData() slot in order to specify the sample annotation, indicating that the first two BAM files form one experimental condition. Then we create our DEXSeqDataSet object.

colData(SE)$condition = c("A", "A", "B")
DEXSeqDataSetFromSE( SE, design= ~ sample + exon + condition:exon )

10.3 Preprocessing using featureCounts

As an alternative to HTSeq, one could use the function featureCounts() to count the number of reads overlapping with each exonic region. The output files from featureCounts() can then be used as input to DEXSeq. In this Github repository, Vivek Bhardwaj provides scripts and documentation that integrate featureCounts() with DEXSeq. He also provides a function to import the count files into R and generate a DEXSeqDataSet object.

10.4 Further accessors

The function geneIDs() returns the gene ID column of the feature data as a character vector, and the function exonIDs() return the exon ID column as a factor.

head( geneIDs(dxd) )

## [1] "FBgn0000256" "FBgn0000256" "FBgn0000256" "FBgn0000256" "FBgn0000256"
## [6] "FBgn0000256"

head( exonIDs(dxd) )

## [1] "E001" "E002" "E003" "E004" "E005" "E006"

These functions are useful for subsetting an DEXSeqDataSet object.

10.5 Overlap operations

The methods subsetByOverlaps() and findOverlaps() have been implemented for the DEXSeqResults object, the query= argument must be a DEXSeqResults object.

interestingRegion = GRanges( "chr2L", 
   IRanges(start=3872658, end=3875302) )
subsetByOverlaps( x=dxr, ranges=interestingRegion )

## 
## LRT p-value: full vs reduced
## 
## DataFrame with 4 rows and 16 columns
##                      groupID   featureID     exonBaseMean
##                  <character> <character>        <numeric>
## FBgn0000256:E001 FBgn0000256        E001 58.3431235035746
## FBgn0000256:E002 FBgn0000256        E002 103.332775062259
## FBgn0000256:E003 FBgn0000256        E003  326.47632106146
## FBgn0000256:E004 FBgn0000256        E004  253.65479841846
##                           dispersion               stat            pvalue
##                            <numeric>          <numeric>         <numeric>
## FBgn0000256:E001  0.0171917182971432 8.235054792749e-06 0.997710330959157
## FBgn0000256:E002 0.00734736582550583    1.5757743019185 0.209370422702931
## FBgn0000256:E003  0.0104806430117485 0.0357357766926327  0.85006218148039
## FBgn0000256:E004  0.0110041990800429  0.169540730538131 0.680520297875874
##                       padj          control        knockdown
##                  <numeric>        <numeric>        <numeric>
## FBgn0000256:E001         1 11.0274690389273 10.8820481621958
## FBgn0000256:E002         1 13.2499830426846 13.5721774023876
## FBgn0000256:E003         1 18.9102859004382 18.6029373357064
## FBgn0000256:E004         1 17.7029156673542 17.2692066316926
##                  log2fold_knockdown_control        control.1
##                                   <numeric>        <numeric>
## FBgn0000256:E001         -0.051359704837787 11.0274690389273
## FBgn0000256:E002          0.103653700142338 13.2499830426846
## FBgn0000256:E003        -0.0895554394887519 18.9102859004382
## FBgn0000256:E004         -0.128321909158976 17.7029156673542
##                       knockdown.1 log2fold_knockdown_control.1
##                         <numeric>                    <numeric>
## FBgn0000256:E001 10.8820481621958           -0.051359704837787
## FBgn0000256:E002 13.5721774023876            0.103653700142338
## FBgn0000256:E003 18.6029373357064          -0.0895554394887519
## FBgn0000256:E004 17.2692066316926           -0.128321909158976
##                              genomicData       countData
##                                <GRanges>        <matrix>
## FBgn0000256:E001 chr2L:3872658-3872947:-    92:28:43:...
## FBgn0000256:E002 chr2L:3873019-3873322:-   124:80:91:...
## FBgn0000256:E003 chr2L:3873385-3874395:- 340:241:262:...
## FBgn0000256:E004 chr2L:3874450-3875302:- 250:189:201:...
##                                                                                                                                               transcripts
##                                                                                                                                                    <list>
## FBgn0000256:E001 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
## FBgn0000256:E002 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
## FBgn0000256:E003 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")
## FBgn0000256:E004 c("FBtr0077511", "FBtr0077513", "FBtr0077512", "FBtr0290077", "FBtr0290079", "FBtr0290078", "FBtr0290082", "FBtr0290080", "FBtr0290081")

findOverlaps( query=dxr, subject=interestingRegion )

## Hits object with 4 hits and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         1           1
##   [2]         2           1
##   [3]         3           1
##   [4]         4           1
##   -------
##   queryLength: 498 / subjectLength: 1

This functions could be useful for further downstream analysis.

10.6 Methodological changes since publication of the paper

In our paper (Anders, Reyes, and Huber 2012), we suggested to fit for each exon a model that includes separately the counts for all the gene’s exons. However, this turned out to be computationally inefficient for genes with many exons, because the many exons required large model matrices, which are computationally expensive to deal with. We have therefore modified the approach: when fitting a model for an exon, we now sum up the counts from all the other exon and use only the total, rather than the individual counts in the model. Now, computation time per exon is independent of the number of other exons in the gene, which improved DEXSeq’s scalability. While the \(p\) values returned by the two approaches are not exactly equal, the differences were very minor in the tests that we performed.

Deviating from the paper’s notation, we now use the index \(i\) to indicate a specific counting bin, with \(i\) running through all counting bins of all genes. The samples are indexed with \(j\), as in the paper. We write \(K_{ij0}\) for the count or reads mapped to counting bin \(i\) in sample \(j\) and \(K_{ij1}\) for the sum of the read counts from all other counting bins in the same gene. Hence, when we write \(K_{ijl}\), the third index \(l\) indicates whether we mean the read count for bin \(i\) (\(l=0\)) or the sum of counts for all other bins of the same gene (\(l=1\)). As before, we fit a GLM of the negative binomial (NB) family $ K_{ijl} (=s_j_{ijl},=_i)$ and we write \(\log_2 \mu_{ijl} = \beta^\text{S}_{ij} + l \beta^\text{E}_{i} + \beta^\text{EC}_{i\rho_j}\).

This model is fit separately for each counting bin \(i\). The coefficient \(\beta^\text{S}_{ij}\) accounts for the sample-specific contribution (factor sample), the term \(\beta^\text{E}_{i}\) is only included if \(l=1\) and hence estimates the logarithm of the ratio \(K_{ij1}/K_{ij0}\) between the counts for all other exons and the counts for the tested exon. As this coefficient is estimated from data from all samples, it can be considered as a measure of “average exon usage”. In the R model formula, it is represented by the term exon with its two levels this (\(l=0\)) and others (\(l=1\)). Finally, the last term, \(\beta^\text{EC}_{i,\rho_j}\), captures the interaction condition:exon, i.e., the change in exon usage if sample \(j\) is from experimental condition group \(\rho(j)\). Here, the first condition, \(\rho=0\), is absorbed in the sample coefficients, i.e., \(\beta^\text{EC}_{i0}\) is fixed to zero and does not appear in the model matrix.

For the dispersion estimation, one dispersion value \(\alpha_i\) is estimated with Cox-Reid-adjusted maximum likelihood using the full model given above. A mean-variance relation is fitted using the individual dispersion values. Finally, the individual values are shrinked towards the fitted values. For more details about this shrinkage approach look at the DESeq2 vignette and/or its manuscript (Love, Huber, and Anders 2014). For the likelihood ratio test, this full model is fit and compared with the fit of the reduced model, which lacks the interaction term \(\beta^\text{EC}_{i\rho_j}\). As described in Section @ref{complexdesigns}, alternative model formulae can be specified.

#Requirements on GTF files

In the initial preprocessing step described in Section @ref{preparing-the-annotation}, the Python script dexseq_prepare_annotation.py is used to convert a GTF file with gene models into a GFF file with collapsed gene models. We recommend to use GTF files downloaded from Ensembl as input for this script, as files from other sources may deviate from the format expected by the script. Hence, if you need to use a GTF or GFF file from another source, you may need to convert it to the expected format. To help with this task, we here give details on how the dexseq_prepare_annotation.py script interprets a GFF file.

• The script only looks at exon lines, i.e., at lines which contain the term exon in the third (“type”) column. All other lines are ignored.
• Of the data in these lines, the information about chromosome, start, end, and strand (1st, 4th, 5th, and 7th column) are used, and, from the last column, the attributes gene_id and transcript_id. The rest is ignored.
• The gene_id attribute is used to see which exons belong to the same gene. It must be called gene_id (and not Parent as in GFF3 files, or GeneID as in some older GFF files), and it must give the same identifier to all exons from the same gene, even if they are from different transcripts of this gene. (This last requirement is not met by GTF files generated by the Table Browser function of the UCSC Genome Browser.)
• The transcript_id attribute is used to build the transcripts attribute in the flattened GFF file, which indicates which transcripts contain the described counting bin. This information is needed only to draw the transcript model at the bottom of the plots when the displayTranscripts= option to plotDEXSeq() is used.

Therefore, converting a GFF file to make it suitable as input to dexseq_prepare_annotation.py amounts to making sure that the exon lines have type exon and that the atributes giving gene ID (or gene symbol) and transcript ID are called gene_id and transcript_id, with this exact spelling. Remember to also take care that the chromosome names match those in your SAM files, and that the coordinates refer to the reference assembly that you used when aligning your reads.