如何优雅的统计基因外显子长度

小编最近在统计基因组内每个基因的外显子长度，原以为非常简单，直接查找外显子的位置计算就可以，但写脚本的时候才发现非常麻烦。因为基因组中很多外显子区域是重合的，粗暴的将每个外显子的长度加在一起是不对的，这时我们可以使用R包"GenomicFeatures "去除外显子重叠的部分，优雅的统计每个基因的外显子长度。

## 安装R包"GenomicFeatures"
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GenomicFeatures")
## 加载R包
library("GenomicFeatures")
## 导入gff3文件
txdb <- makeTxDbFromGFF("genome.gff3",format="gff3")
## 获取外显子位置
exons_gene <- exonsBy(txdb, by = "gene")
## 去除外显子重叠部分，计算外显子长度
exons_gene_len <- lapply(exons_gene,function(x){sum(width(reduce(x)))})

简单几步，我们就可以准确地计算出每个基因的外显子长度啦！

"GenomicFeatures "还有非常多的功能，如提取基因的启动子序列、获得基因内含子的长度等等，大家可以查看它的manual，学习更多的用法！

本文R语言代码引自：

http://www.bioinfo-scrounger.com/archives/342/ 作者: Kai

参考资料：

http://www.bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf