DEXseq: How to do differential exon usage of genes with overlapping exons
Entering edit mode
Last seen 6 months ago

I followed the DEXSeq tutorial described here for performing differential exon usage analysis. When looking at the DEU plot for my gene of interest many of the exons were not shown. It turned out that a lncRNA having overlapping exons with my gene of interest and that many exons were removed during exon binning when I set the parameter to TRUE.

R code:

# Create transcript data base
txdb = GenomicFeatures::makeTxDbFromGFF(gtf_file)

# Collapse the gene models into counting bins
exonic_parts = GenomicFeatures::exonicParts(txdb, = TRUE)

# Create references to BAM files
bam_filelist = Rsamtools::BamFileList(bam_files, index=character(), yieldSize=100000, obeyQname=TRUE)

# Count the reads overlapping the bins
ese = GenomicAlignments::summarizeOverlaps(exonic_parts, bam_filelist, mode="Union", singleEnd=isSingleEnd, ignore.strand=TRUE, inter.feature=FALSE, fragments=!isSingleEnd)

# Convert SummarizedExperiment to DEXSeqDataSet
dxd = DEXSeq::DEXSeqDataSetFromSE(ese, design = design)

When I set the parameter to FALSE an error occurs during the build of a DEXSeq object with DEXSeqDataSetFromSE:

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘NSBS’ for signature ‘"CompressedIntegerList"’

I suppose it's because there is no exon_part column (reflects the exon bin order in a gene) included in rowRanges(ese) because a sequence can map to multiple genes. I could find a workaround by splitting exon bins with multiple comma-seperated genes in to multiple rows.

However, the problem that I see now is that DEXSeq can't differentiate between read counts that belong to lncRNA and the read counts that belong to my gene of interest. I.e. for my gene of interest DEXseq would always include the read counts of lncRNA for the overlapping exons and the same for the lncRNA.

The only solution I can think of is to use a transcript quantification tool like Salmon instead of a genome mapper and to build a DEXSeq object from it. Is there something I'm missing? I am happy to hear your opinion.

GenomicFeatures DEXSeq • 803 views

Login before adding your answer.

Traffic: 283 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6