1
$\begingroup$

Goal

I am trying to get a distribution of Exons and intron sizes in Three-Spined Stickleback (Gasterosteus aculeatus).

Data downloaded

I downloaded some data from Ensembl. More precisely, I went there, selected "structure" in the "attributes" and selected

- Ensembl Gene ID
- Exon Chr Start (bp)
- Exon Chr End (bp)

There are two issues in these data:

  1. Some exons end before they start
  2. Some exons overlap

The first point is I assumed due to reversion, so I just reversed the start and end position when there was such inversion. The second point is more problematic as I don't know what such overlap could mean realistically.

Can you help finding the distribution of exon and intron sizes in three-spined stickleback?


What I did with the above data

Preparation

# read data

d = read.table(file.choose(), header=TRUE)
names(d)= c("ID", "start", "end")

# Reverse exons when needed

toReverse = which(d$start > d$end)
s = d$start[toReverse]
d$start[toReverse] = d$end[toReverse]
d$end[toReverse] = s

Exon Sizes -> Looks good

# Exon sizes

ExonSizes = d$end - d$start
hist(ExonSizes) # Cool, looks good!

Introns Sizes -> Looks bad

# Intron sizes
# This is a little slow. One might have better to switch to C++.
# The idea is to look at the distance between the end of an exon and the start of the next exon within each gene.

d$ID = paste(d$ID)
IntronSizes = c()
uniquegenes = unique(d$ID)
nbgenes = length(uniquegenes)
i=0
exon = 0
for (gene in uniquegenes)
{
    i=i+1
    cat(paste0(i," / ",nbgenes,"\n"))

    while (TRUE)
    {
        exon=exon+1
        if (d$ID[exon] == gene) {
			IntronSizes = append(IntronSizes, d[exon+1,]$start - d[exon,]$end)
        } else 
        {
            break
        }
    }
}
hist(IntronSizes) # There are negative values. I don't know how to interpret them. 

hist(log(abs(IntronSizes))) # That looks good but I doubt it makes much sense!
$\endgroup$
  • 1
    $\begingroup$ When it looks like an exon ends before it begins, thats because its on the minus strand. Strand information should also be in your gtf file, so you can check for this just in case. As for overlapping, sometimes an exon has multiple splice sites. Also, genes can overlap. Both cases would cause exons to overlap. Nothing weird going on there, you can include them without worry $\endgroup$ – von Mises Jun 24 '16 at 15:54
3
$\begingroup$

At the command line:

wget ftp://ftp.ensembl.org/pub/release-84/gtf/gasterosteus_aculeatus/Gasterosteus_aculeatus.BROADS1.84.gtf.gz
gunzip Gasterosteus_aculeatus.BROADS1.84.gtf.gz

Then in R:

library(GenomicFeatures)
txdb = makeTxDbFromGFF("Gasterosteus_aculeatus.BROADS1.84.gtf")
hist(log10(width(unique(exons(txdb))))) # exons
hist(log10(width(unique(unlist(intronsByTranscript(txdb)))))) # introns

Note that some of the annotated exons (and the introns between them) are incredibly unlikely to be correct. For example, there are 1 base long exons and introns. I don't think anyone actually believes those are correct, but the distribution is probably otherwise approximately correct.

Edit: For what it's worth, there are no negative sized exons in the GTF file. My guess is that the overlapping exons are from different transcripts (or genes on opposite strands).

Edit2: If you want to get your introns from a "union gene model", then use something like reduce(exonsBy(txdb, by='gene')) and then lapply gaps() to that. The results will be correct and the process will probably take less time than what you've been trying.

$\endgroup$
  • $\begingroup$ Perfect +1 Thanks a lot! For future users, please note that GenomicFeatures is a bioconductor package and that with MAC OSX wget must be replaced by curl $\endgroup$ – Remi.b Jun 1 '16 at 21:44
  • $\begingroup$ May I ask you to further explain your second edit? What is a "union gene model"? Do you mean that I might want to get a distribution of per-gene average exon length? $\endgroup$ – Remi.b Jun 1 '16 at 21:52
  • $\begingroup$ The explanation of a union gene model ends up being a bit long for a comment. Have a look at the right part of subfigure A here: nature.com/nbt/journal/v31/n1/images/nbt.2450-F1.jpg The "exon-union model" is the same thing as a "union gene model". In this case, you're still getting the distribution of all introns in the genome, you're just redefining "intron" to not be something biologically correct (though this redefinition can be useful for some purposes). $\endgroup$ – Devon Ryan Jun 1 '16 at 22:40
  • $\begingroup$ Oh sure I got it. I did not thing of that. I will need to go with a union gene model for my purpose. I'll try to get that as soon as I'll understand what kind of objects I am dealing with! $\endgroup$ – Remi.b Jun 1 '16 at 23:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.