I have downloaded the latest GTF file (Homo_sapiens.GRCh37.75.gtf) from Ensembl，I was trying to find all protein coding genes and calculte the exome size, I did as follows:
awk '{if($3=="gene" && $2=="protein_coding"){print $0}}' Homo_sapiens.GRCh37.75.gtf
I found 22810 protein_coding genes, with a total length of 1,395,684,274 bp ( | awk '{sum+=$5-$4}END{print sum}'); 1.3G seems too large for me.
Similarly when I search for all exons
awk '{if($3=="exon" && $2=="protein_coding"){print $0}}' Homo_sapiens.GRCh37.75.gtf
I found 809933 exons, with a total length of 191,357,777 bp
As far as I remember, in a coure I attended before, we have tried R to calculate the total exome size of human, which is ~30M, why my calculation here is so big? did I made some mistakes?
Thanks for your help.
2 answers
If a gene has multiple transcripts then you'll end up double/triple/etc. counting some of its exons with that method. In R, you can import the GTF, split() by gene_id and then reduce() that to get non-overlapping exons for this purpose. That'd be a real pain in awk.
Log in to answer this question.
-
Custom reference sequence length discrepancy with STAR/Salmon/gffread
written by ki •
-
Drosophila melanogaster dm3 intron percentage too high
written by Christian •
-
how to calculate gene length for TPM
written by adi.gershon1 •
-
Qualimap whole exome sequencing depth of coverage - use whole genome or exome as reference?
written by rebeliscu •
-
samtools - Finding Average Fragment size for pair?
written by oars
-
Tool to compare whole genebody enrichment between different samples (MeDip-seq)
written by florian.noack •
-
Scaffolding with SSPACE returned scaffold the same as initial contigs set
written by pbigbig
-
Calculating Exome Coverage
written by Coryza •
-
Tutorial: Extract Total Non-Overlapping Exon Length Per Gene With Bioconductor
written by Irsan
-
Total Contig Length From Velvet Output
written by deepthithomaskannan
I have another naive questions, since intron is not included in the gtf (is it used to?), how I infer the intron regions, regions in the gene but without annotion should be the intron region? am I right? or these are other easy method?
Introns aren't explicitly included since they're just the regions between the exons. You can add them in and then use the same methods to calculate their size. Example scripts are provided here in R (from me) and perl (from Alejandro Reyes).
thanks a lot for your help