100-用R语言统计人类外显子信息
刘小泽写于19.4.8
今天学习了GenomicRange，那么利用Ensembl的R包统计人类外显子信息就很方便了
首先下载包
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
BiocManager::install("EnsDb.Hsapiens.v86", version = "3.8")
可以看到最新版本是V86，那么对应的是hg19还是38呢？这里做一个补充
对应NCBI(UCSC)：Ensembl
# GRCh36 (hg18): ENSEMBL release_52.
# GRCh37 (hg19): ENSEMBL release_59/61/64/68/69/75.
# GRCh38 (hg38): ENSEMBL release_76/77/78/80/81/82...
因此这里下载的最新版是hg38版本
加载外显子信息
library(EnsDb.Hsapiens.v86)
ensembl.hg38 <- EnsDb.Hsapiens.v86
ensembl.hg38.exon <- exons(ensembl.hg38) #exon比较多，需要加载1分钟左右
# 看下有多少个exon
length(ensembl.hg38.exon) # 748400个
# 长度分布-v1
hist(width(ensembl.hg38.exon)) # 发现图中受极大值影响，大部分数据展示不出来
# （进行改进）长度分布-v2
hist(log2(width(ensembl.hg38.exon)+1)) # 大部分在2^7~2^8左右，也就是128-256bp之间
数据探索
mean(width(ensembl.hg38.exon)) #因此human exon 平均长度是310，这可以当成一个常识来记忆
median(width(ensembl.hg38.exon)) #中位数142
那么问题来了：如何求全部exon的总长度？
# 可能最直接想到的方法就是
sum(width(ensembl.hg38.exon)) #234,797,689大概234M
但是这样并不对，因为忽略了exon之间的overlap情况，需要做的是先reduce，再width
sum(width(GenomicRanges::reduce(ensembl.hg38.exon))) # 143,390,351
但是人类的外显子不应该是30M左右吗？
