Similar as with differential gene expression, we need to make sure that observed differences of exon usage values between conditions are statistically significant, i.e., are sufficiently unlikely to be just due to random fluctuations such as those seen even between samples from the same condition, i.e., between replicates. To this end, DEXSeq assesses the strength of these fluctuations (quantified by the so-called dispersion) by comparing replicates before comparing the averages between the sample groups.
The preceding description is somewhat simplified (and perhaps over-simplified), and we recommend that users consult the publication for a more complete description, as well as Appendix 10.6 of this vignette, which describes how the current implementation of DEXSeq differs from the original approach described in the paper. Nevertheless, two important aspects should be mentioned already here: First, DEXSeq does not actually work on the exon usage ratios, but on the counts in the numerator and denominator, to be able to make use of the information that is contained in the magnitude of count values. For example, 3000 reads versus 1000 reads is the same ratio as 3 reads versus 1 read, but the latter is a far less reliable estimate of the underlying true value, because of statistical sampling. Second, DEXSeq is not limited to simple two-group comparisons; rather, it uses so-called generalized linear models (GLMs) to permit ANOVA-like analysis of potentially complex experimental designs.
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.
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).
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)
##  "dexseq_count.py" "dexseq_prepare_annotation.py"
system.file( "python_scripts", package="DEXSeq", mustWork=TRUE )
##  "/tmp/RtmpQMiG9x/Rinst452e12f3b250/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.
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.
For each SAM file, we next count the number of reads that overlap with each of the exon counting bins defined in the flattened GFF file. This is done with the script python_count.py:
python /path/to/library/DEXSeq/python_scripts/dexseq_count.py Dmel.BDGP5.25.62.DEXSeq.chr.gff untreated1.sam untreated1fb.txt
This command (again, to be entered as a single line) expects two files in the current working directory, namely the GFF file produced in the previous step (here Dmel.BDGP5.25.62.DEXSeq.chr.gff) and a SAM file with the aligned reads from a sample (here the file untreated1.sam with the aligned reads from the first control sample). The command generates an output file, here called untreated1fb.txt, with one line for each exon counting bin defined in the flattened GFF file. The lines contain the exon counting bin IDs (which are composed of gene IDs and exon bin numbers), followed by a integer number which indicates the number of reads that were aligned such that they overlap with the counting bin.
Use the script multiple times to produce a count file from each of your SAM files.
There are a number of crucial points to pay attention to when using the python_count.py script:
Paired-end data: If your data is from a paired-end sequencing run, you need to add the option
-p yes to
the command to call the script. (As usual, options have to be placed before the file names, surrounded by spaces.) In addition, the SAM file needs to be sorted, either by read name or by position. Most aligners produce sorted SAM files; if your SAM file is not sorted, use
samtools sort -n to sort by read name (or
samtools sort) to sort by position. (See Anders et al. (2013), if you need further explanations on how to sort SAM files.) Use the option
-r pos or
-r name to indicate whether your paired-end data is sorted by alignment position or by read name.
Strandedness: By default, the counting script assumes your library to be strand-specific, i.e., reads are aligned to the same strand as the gene they originate from. If you have used a library preparation protocol that does not preserve strand information (i.e., reads from a given gene can appear equally likely on either strand), you need to inform the script by specifying the option
-s no. If your library preparation protocol reverses the strand (i.e., reads appear on the strand opposite to their gene of origin), use
-s reverse. In case of paired-end data, the default (
-s yes) means that the read from the first sequence pass is on the same strand as the gene and the read from the second pass on the opposite strand (“forward-reverse” or “fr” order in the parlance of the Bowtie/TopHat manual) and the options
-s reverse specifies the opposite case.
SAM and BAM files: By default, the script expects its input to be in plain-text SAM format.
However, it can also read BAM files, i.e., files in the the compressed binary variant of the SAM format.
If you wish to do so, use the option
-f bam. This works only if you have installed the Python
Alignment quality: The scripts takes a further option,
-a to specify the minimum alignment quality (as given in the fifth column of the SAM file). All reads with a lower quality than specified (with default
-a 10) are skipped.
Help pages: Calling either script without arguments displays a help page with an overview of all options and arguments.
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)
##  "treated1fb.txt" "treated2fb.txt" "treated3fb.txt" "untreated1fb.txt" ##  "untreated2fb.txt" "untreated3fb.txt" "untreated4fb.txt"
flattenedFile = list.files(inDir, pattern="gff$", full.names=TRUE) basename(flattenedFile)
##  "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.
## 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
library( "DEXSeq" ) dxd = DEXSeqDataSetFromHTSeq( countFiles, sampleData=sampleTable, design= ~ sample + exon + condition:exon, flattenedfile=flattenedFile )
## Warning in DESeqDataSet(rse, design, ignoreRank = TRUE): some variables in ## design formula are characters, converting to factors
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.
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)[] 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
## DataFrame with 14 rows and 4 columns ## sample condition libType exon ## <factor> <factor> <character> <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] [,12] ## FBgn0000256:E001 92 28 43 54 131 51 49 1390 829 923 1115 2495 ## FBgn0000256:E002 124 80 91 76 224 82 95 1358 777 875 1093 2402 ## FBgn0000256:E003 340 241 262 347 670 260 297 1142 616 704 822 1956 ## FBgn0000256:E004 250 189 201 219 507 242 250 1232 668 765 950 2119 ## FBgn0000256:E005 96 38 39 71 76 57 62 1386 819 927 1098 2550 ## [,13] [,14] ## FBgn0000256:E001 1054 1052 ## FBgn0000256:E002 1023 1006 ## FBgn0000256:E003 845 804 ## FBgn0000256:E004 863 851 ## FBgn0000256:E005 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 )
## $others ##  8 9 10 11 12 13 14 ## ## $this ##  1 2 3 4 5 6 7
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.000 1250.67 ## FBgn0000256:E002 110.286 2769.57 ## FBgn0000256:E003 345.286 22147.90 ## transcripts ## <list> ## FBgn0000256:E001 FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E002 FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E003 FBtr0077511,FBtr0077513,FBtr0077512,... ## ------- ## 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
sampleAnnotation( dxd )
## DataFrame with 7 rows and 3 columns ## sample condition libType ## <factor> <factor> <character> ## 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.
dxd = estimateSizeFactors( dxd )
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 )
Having the dispersion estimates and the size factors, we can now test for differential exon usage.
For each gene, DEXSeq fits a generalized linear model with the formula
~sample + exon + condition:exon and compare it to the smaller model (the null model)
~ sample + exon.
In these formulae (which use the standard notation for linear model formulae in R; consult a text book
on R if you are unfamiliar with it),
sample is a factor with different levels for each sample,
condition is the factor of experimental conditions that we defined when constructing the
DEXSeqDataSet object at the beginning of the analysis, and
exon is a factor with
others, that were specified when we generated our DEXSeqDataSet
object. The two models described by these formulae are fit for each counting bin, where the data
supplied to the fit comprise two read count values for each sample, corresponding to the
two levels of the
exon factor: the number of reads mapping to the bin in question
this), and the sum of the read counts from all other bins of the same gene
others). Note that this approach differs from the approach described in
the paper (Anders, Reyes, and Huber (2012)) and used in older versions of DEXSeq;
see Appendix 10.6
for further discussion.
We have an interaction term
condition:exon, but denote no main effect for
Note, however, that all observations from the same sample are also from the same condition, i.e., the
condition main effects are absorbed in the
sample main effects, because the
factor is nested within the
The deviances of both fits are compared using a \(\chi^2\)-distribution, providing a \(p\) value. Based on this \(p\)-value, we can decide whether the null model is sufficient to explain the data, or whether it may be rejected in favour of the alternative model, which contains an interaction coefficient for
condition:exon. The latter means that the fraction of the gene’s reads that fall onto the exon under the test differs significantly between the experimental conditions.
testForDEU() performs these tests for each exon in each gene.
dxd = testForDEU( dxd )
The resulting DEXSeqDataSet object contains slots with information regarding the test.
For some uses, we may also want to estimate relative exon usage fold changes. To this end,
estimateExonFoldChanges(). Exon usage fold changes are
calculated based on the coefficients of a GLM fit that uses the formula
count ~ condition + exon + condition:exon
condition can be replaced with any of the column names of
colData() (see man pages for details). The resulting coefficients allow the estimation of changes on the usage of exons across different conditions. Note that the differences on exon usage are distinguished from gene expression differences across conditions. For example, consider the case of a gene that is differentially expressed between conditions and has one exon that is differentially used between conditions. From the coefficients of the fitted model, it is possible to distinguish overall gene expression effects, that alter the counts from all the exons, from exon usage effects, that are captured by the interaction term
condition:exon and that affect the each exon individually.
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")
So far in the pipeline, the intermediate and final results have been stored in the meta data of a DEXSeqDataSet object, they can be accessed using the function
mcols(). In order to summarize the results without showing the values from intermediate steps, we call the function
DEXSeqResults(). The result is a DEXSeqResults object, which is a subclass of a DataFrame object.
dxr1 = DEXSeqResults( dxd ) dxr1
## ## LRT p-value: full vs reduced ## ## DataFrame with 498 rows and 13 columns ## groupID featureID exonBaseMean dispersion stat ## <character> <character> <numeric> <numeric> <numeric> ## FBgn0000256:E001 FBgn0000256 E001 58.3431 0.01719168 8.23378e-06 ## FBgn0000256:E002 FBgn0000256 E002 103.3328 0.00734733 1.57578e+00 ## FBgn0000256:E003 FBgn0000256 E003 326.4763 0.01048062 3.57359e-02 ## FBgn0000256:E004 FBgn0000256 E004 253.6548 0.01100418 1.69541e-01 ## FBgn0000256:E005 FBgn0000256 E005 60.6377 0.04402967 2.91053e-02 ## ... ... ... ... ... ... ## FBgn0261573:E012 FBgn0261573 E012 23.08422 0.0219839 8.400797 ## FBgn0261573:E013 FBgn0261573 E013 9.79723 0.2475513 1.155232 ## FBgn0261573:E014 FBgn0261573 E014 87.49678 0.0330407 1.118567 ## FBgn0261573:E015 FBgn0261573 E015 268.25566 0.0120126 2.598314 ## FBgn0261573:E016 FBgn0261573 E016 304.15134 0.0998527 0.146583 ## pvalue padj control knockdown ## <numeric> <numeric> <numeric> <numeric> ## FBgn0000256:E001 0.997711 1 11.0275 10.8820 ## FBgn0000256:E002 0.209370 1 13.2500 13.5722 ## FBgn0000256:E003 0.850062 1 18.9103 18.6029 ## FBgn0000256:E004 0.680520 1 17.7029 17.2692 ## FBgn0000256:E005 0.864536 1 10.9937 11.1006 ## ... ... ... ... ... ## FBgn0261573:E012 0.00375057 0.0977353 6.61692 8.47412 ## FBgn0261573:E013 0.28245660 1.0000000 5.97080 3.79997 ## FBgn0261573:E014 0.29022722 1.0000000 13.19130 11.94851 ## FBgn0261573:E015 0.10697745 0.9656443 18.42453 17.26506 ## FBgn0261573:E016 0.70182202 1.0000000 18.91296 18.14410 ## log2fold_knockdown_control genomicData ## <numeric> <GRanges> ## FBgn0000256:E001 -0.0513597 chr2L:3872658-3872947:- ## FBgn0000256:E002 0.1036537 chr2L:3873019-3873322:- ## FBgn0000256:E003 -0.0895554 chr2L:3873385-3874395:- ## FBgn0000256:E004 -0.1283219 chr2L:3874450-3875302:- ## FBgn0000256:E005 0.0375698 chr2L:3878895-3879067:- ## ... ... ... ## FBgn0261573:E012 0.832682 chrX:19421654-19421867:+ ## FBgn0261573:E013 -1.395561 chrX:19422668-19422673:+ ## FBgn0261573:E014 -0.410999 chrX:19422674-19422856:+ ## FBgn0261573:E015 -0.341490 chrX:19422927-19423634:+ ## FBgn0261573:E016 -0.224595 chrX:19423707-19424937:+ ## countData transcripts ## <matrix> <list> ## FBgn0000256:E001 92:28:43:... FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E002 124:80:91:... FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E003 340:241:262:... FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E004 250:189:201:... FBtr0077511,FBtr0077513,FBtr0077512,... ## FBgn0000256:E005 96:38:39:... FBtr0077511,FBtr0077513,FBtr0077512,... ## ... ... ... ## FBgn0261573:E012 37:23:38:... FBtr0302863 ## FBgn0261573:E013 8:3:6:... FBtr0302863,FBtr0302864 ## FBgn0261573:E014 75:66:92:... FBtr0302862,FBtr0302863,FBtr0302864 ## FBgn0261573:E015 264:234:245:... FBtr0302862,FBtr0302863,FBtr0302864 ## FBgn0261573:E016 611:187:188:... FBtr0302862,FBtr0302863,FBtr0302864
The description of each of the column of the object DEXSeqResults can be found in the metadata columns.
##  "group/gene identifier" ##  "feature/exon identifier" ##  "mean of the counts across samples in each feature/exon" ##  "exon dispersion estimate" ##  "LRT statistic: full vs reduced" ##  "LRT p-value: full vs reduced" ##  "BH adjusted p-values" ##  "exon usage coefficient" ##  "exon usage coefficient" ##  "relative exon usage fold change" ##  "GRanges object of the coordinates of the exon/feature" ##  "matrix of integer counts, of each column containing a sample" ##  "list of transcripts overlapping with the exon"
From this object, we can ask how many exonic regions are significant with a false discovery rate of 10%:
table ( dxr1$padj < 0.1 )
## ## FALSE TRUE ## 426 17
We may also ask how many genes are affected
table ( tapply( dxr1$padj < 0.1, dxr1$groupID, any ) )
## ## FALSE TRUE ## 20 9
Remember that our example data set contains only a selection of genes. We have chosen these to contain interesting cases; so the fact that such a large fraction of genes is significantly affected here is not typical.
To see how the power to detect differential exon usage depends on the number of reads that
map to an exon, a so-called MA plot is useful, which plots the logarithm of fold change versus average
normalized count per exon and marks by red colour the exons which are considered significant; here,
the exons with an adjusted \(p\) values of less than \(0.1\) (Figure 2. There is of course nothing special about the number \(0.1\), and you can specify other thresholds in the call to
plotMA( dxr1, cex=0.8 )