Question: Question: Calculating each gene's mean of intron length
gravatar for fabregas23213
14 months ago by
fabregas232130 wrote:

Hi, I want to calculate the mean of intron length of each gene in mouse(mm10) or human(hg19 or hg38). The rationale of doing this is trying to categorize genes into 4 or 5 bins based on their intron length. As the previously published paper "", in Table 1, they classified genes into 4 groups depending on their intron length(0-1kb, 1-10kb,10-100kb,>100kb) and count the gene number in each class. Could anyone give me any suggestion to do this?

sequencing genome • 486 views
ADD COMMENTlink modified 14 months ago by shawn.w.foley1.2k • written 14 months ago by fabregas232130

This question appears to be identical to: Calculating each gene's mean of intron length. It is not a good practice to post the same question using two different accounts.

ADD REPLYlink written 14 months ago by genomax92k
gravatar for shawn.w.foley
14 months ago by
shawn.w.foley1.2k wrote:

You can obtain a .bed file with gene-body coordinates and a .bed file with exon coordinates from UCSC or Ensembl. After that you can use bedtools subtract to subtract the exons from the gene bodies - this will give you an output file with all of the intron coordinates.

bedtools subtract -a gene_body.bed -b exons.bed -s > introns.bed

After that you can calculate the length and average the size across genes using whichever language you prefer. Personally, I'd perform this in R:

Input format (introns.bed):

chr1    11873   12226   DDX11L1.1   .   +
chr1    12612   12720   DDX11L1.1   .   +
chr1    13220   14408   DDX11L1.1   .   +

I'd run the code:

R --vanilla
inFile <- read.table('introns.bed',header=F,,sep='\t')
colnames(inFile) <- c('chr','start','stop','gene','score','strand')
inFile$intronLength <- inFile$stop - inFile$start
meanIntronLength <- aggregate(inFile$intronLength,list(inFile$gene),mean)

This will read in the file, calculate the intron length, then aggregate via gene name by taking the mean of the intron lengths. The output file will be a two-column file with gene name and intron length.

ADD COMMENTlink written 14 months ago by shawn.w.foley1.2k
gravatar for chouaib23213
14 months ago by
chouaib232130 wrote:

using BioAlcidaeJdk

I'm treating each intron of each transcript as a distinct entity.

$ wget -O - -q "" | gunzip -c |\
java -jar dist/bioalcidaejdk.jar -F GTF -e 'stream().forEach(G->println(G.getId()+"\t"+G.getGeneName()+"\t"+G.getTranscripts().stream().flatMap(T->T.getIntrons().stream()).mapToInt(I->I.getLengthOnReference()).average().orElse(-99) ));' 

ENSMUSG00000029836  Cbx3    2392.5
ENSMUSG00000029837  1700111E14Rik   11957.25
ENSMUSG00000029830  Svopl   3663.7948717948716
ENSMUSG00000005864  Cnga2   3509.5714285714284
ENSMUSG00000029831  Npvf    1478.0
ENSMUSG00000029832  Nfe2l3  8307.166666666666
ENSMUSG00000029833  Trim24  5817.84
ENSMUSG00000029838  Ptn 23120.14285714286
ENSMUSG00000030827  Fgf21   327.0
ENSMUSG00000030825  Hsd17b14    1523.090909090909
ENSMUSG00000030826  Bcat2   1856.8392857142858
ENSMUSG00000030823  9130019O22Rik   533.0
ENSMUSG00000030824  Nucb1   1569.7358490566037
ENSMUSG00000030822  Prr14   622.8695652173913
ENSMUSG00000005871  Apc 6496.068493150685
ENSMUSG00000029823  Luc7l2  5344.425
ENSMUSG00000029826  Zc3hav1 3994.5384615384614
ENSMUSG00000029821  Gsdme   6530.684210526316
ENSMUSG00000029822  Osbpl3  5717.4789915966385
ENSMUSG00000005873  Reep5   8005.071428571428
ADD COMMENTlink modified 14 months ago by genomax92k • written 14 months ago by chouaib232130

It looks like you are re-posting @Pierre's answer here. You should provide proper attribution if you are going to use someone else's code.

ADD REPLYlink modified 14 months ago • written 14 months ago by genomax92k

ADD REPLYlink written 14 months ago by Pierre Lindenbaum131k

The question is asked separately by two users. Possible to merge ?

ADD REPLYlink written 14 months ago by geek_y11k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2016 users visited in the last hour