本文包括内容：
1. 基因长度定义
2. 常用注释包简介
3. 基因长度计算:用GenomicFeatures包 计算TxDb.Hsapiens.UCSC.hg38.knownGene中四种不同定义的基因长度。
参考：
https://blog.csdn.net/Candle_light/article/details/90598308
https://mp.weixin.qq.com/s?__biz=MzAxMDkxODM1Ng==&mid=2247488656&idx=1&sn=c2472cca80a0b0f8b831fbf5c69f2a34&chksm=9b48542bac3fdd3dd99474c42b2af092e462842007e21b28724936a3c9fc9f9c14625d355fd1&scene=21#wechat_redirect
https://mp.weixin.qq.com/s?__biz=MzU4NjU4ODQ2MQ==&mid=2247490627&idx=1&sn=f23242af5baa6cd6c07ff3558d65c97b&chksm=fdf85401ca8fdd17017107413706c4c259c91385ddc17e0778539d327f2f1d49f0f5072ea688&scene=21#wechat_redirect
1. 基因长度定义
非冗余外显子(EXON)长度之和;
挑选基因的最长转录本 ;
选取多个转录本长度的平均值;
非冗余 CDS 长度之和
2. 常用注释包简介
2.1 TxDb包
包含位置信息。包里的注释内容，可以通过包名来判断。包名的格式是：TxDb.Species.Source.Build.Table
例如，TxDb.Hsapiens.UCSC.hg19.knownGene，代表：
物种是 Homo sapiens
来源于 UCSC genome browser
基因组版本是 hg19 (their version of GRCh37)
包含的注释内容是 knownGene table
1.2 EnsDb包
EnsDb和TxDb类似，只不过它是基于Ensembl数据库的
例如：EnsDb.Hsapiens.v79
3. GRanges，Ranges
Ranges 对象功能十分强大，可以让我们基于基因组位置信息轻松选取出相关数据。GRanges和GRangesLists对象，就类似于data.frame和List数据结构，我们可以使用[符号去选取子集。
3.1 GRanges
注释包比较常用的一个功能，是将注释包中的位置信息提取出来，添加到GRanges或者GRangesList对象中。
- 将GRanges对象写成txt
exon_txdb=exons(txdb)
tmp=as.data.frame(exon_txdb)
write.table(tmp,"exon.pos",row.names=F)
genes_txdb=genes(txdb)
tmp=as.data.frame(genes_txdb)
write.table(tmp,"gene.pos",row.names=F)
- 其他操作
seqnames(exon_txdb)
ranges(exon_txdb) #返回外显子的起始终止位点，长度，以及其它信息
strand(exon_txdb)#返回外显子的正负链信息，要么在正链要么在负链
mcols(exon_txdb)#返回exon的id编号，1到724366个
seqlengths(exon_txdb)#返回每条染色体的长度信息 names,length
GRanges对象还有很多其它类型的操作，非常好玩的，split,shift,resize,flank,reduce,gaps,disjoin,coverage
其它求交集并集和都可以用,union,intersect,setdiff,pintersect,psetdiff
3.2 GRangesLists
以List形式呈现，每个基因所对应的转录本，在基因组上的坐标
3.3 其他
transcripts(), genes()函数选取位置信息， 编码区可以使用cds(),promoters(),exons()来选取
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
genes(txdb)#29794个基因
exons(txdb)#724366 个exon
transcripts(txdb) #272352个转录本
cds(txdb)#306522个编码区
transcriptsBy(txdb,by="gene")#对象按照gene来对转录本分组
exonsBy,cdsBy来对它进行处理，都会返回Granges对象
4. 根据四种定义计算基因长度
目前了解到两种方法:1) 使用R包GenomicFeatures 2)使用R包 rtracklayer；
首先，载入txdb文件
##创建txdb文件
if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",format="gtf")
##或者下载
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
#### R包GenomicFeatures####
library(GenomicFeatures)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
##① 非冗余外显子(EXON)长度之和
exons.list.per.gene <- exonsBy(txdb,by="gene")
head(as.data.frame(exons.list.per.gene))
group group_name seqnames start end width strand exon_id
1 1 1 chr19 58345178 58347029 1852 - 603914
2 1 1 chr19 58345183 58347029 1847 - 603915
3 1 1 chr19 58346854 58347029 176 - 603916
4 1 1 chr19 58346858 58347029 172 - 603917
5 1 1 chr19 58346860 58347029 170 - 603918
6 1 1 chr19 58347353 58347634 282 - 603919
exon_name
1 <NA>
2 <NA>
3 <NA>
4 <NA>
5 <NA>
6 <NA>
# 然后计算每个基因的非冗余外显子长度(widths)并加和。
exonic.gene.sizes <- sum(width(GenomicRanges::reduce(exons.list.per.gene)))
exonic.gene.sizes
head(as.data.frame(exonic.gene.sizes))
1 6528
10 1834
100 12177
1000 8583
10000 22459
100008586 693
gle <- data.frame(gene_id=names(exonic.gene.sizes),
length=exonic.gene.sizes) ##31456个基因
gene_id length
1 1 6528
10 10 1834
100 100 12177
1000 1000 8583
10000 10000 22459
100008586 100008586 693
##②非冗余cds长度之和
cds.list.per.gene <- cdsBy(txdb,by="gene")
head(as.data.frame(cds.list.per.gene))
group group_name seqnames start end width strand cds_id
1 1 1 chr19 58347022 58347029 8 - 252822
2 1 1 chr19 58347353 58347640 288 - 252823
3 1 1 chr19 58350370 58350651 282 - 252824
4 1 1 chr19 58350594 58350651 58 - 252825
5 1 1 chr19 58351391 58351687 297 - 252826
6 1 1 chr19 58352098 58352184 87 - 252827
cds_name
1 <NA>
2 <NA>
3 <NA>
4 <NA>
5 <NA>
6 <NA>
# 然后计算每个基因的非冗余外显子长度(widths)并加和。
cds.gene.sizes <- sum(width(GenomicRanges::reduce(cds.list.per.gene)))
head(as.data.frame(cds.gene.sizes))
cds.gene.sizes
1 1575
10 873
100 1979
1000 2978
10000 3296
100008586 354
glc <- data.frame(gene_id=names(cds.gene.sizes),
length=cds.gene.sizes)#19626个基因
gene_id length
1 1 1575
10 10 873
100 100 1979
1000 1000 2978
10000 10000 3296
100008586 100008586 354
##③挑选基因的最长转录本
t_l=transcriptLengths(txdb) #272352
t_l=na.omit(t_l)#229531
t_l=t_l[order(t_l$tx_len,decreasing = T),]# 这里把同样的基因的多个转录本长度排序了
head(t_l)
tx_id tx_name gene_id nexon tx_len
38851 38851 ENST00000589042.5 7273 363 109224
38852 38852 ENST00000591111.5 7273 313 104301
38849 38849 ENST00000342992.11 7273 312 101520
138913 138913 ENST00000597346.1 10984 1 91667
38850 38850 ENST00000460472.6 7273 191 82029
38853 38853 ENST00000342175.11 7273 191 81357
t_l=t_l[!duplicated(t_l$gene_id),]#去重，直接保留转录本最长的
t_l=t_l[order(t_l$gene_id,decreasing = F),]
glm=data.frame(gene_id=t_l$gene_id,length=t_l$tx_len)
head(glm)
gene_id length
1 1 3382
2 10 1285
3 100 5965
4 1000 4016
5 10000 7950
6 100008586 585
##④ 选取多个转录本长度的平均值
t_l=transcriptLengths(txdb) #272352
t_l=na.omit(t_l)#229531
t_l=t_l[order(t_l$gene_id,decreasing = T),]
gla <- t_l %>% group_by(gene_id) %>%
summarise(tx_ave_len = round(mean(tx_len),0))#31456个基因
gla <-as.data.frame(gla)
head(gla)
gene_id tx_ave_len
1 1 1842
2 10 850
3 100 1708
4 1000 2744
5 10000 2436
6 100008586 573
比较不同方法得到的基因长度:
a = gle %>% inner_join(glc,"gene_id") %>%
inner_join(glm,"gene_id") %>%
inner_join(gla,"gene_id")
colnames(a) = c("gene_id","exon","CDS","max","average") # 19626个基因
head(a)
gene_id exon CDS max average
1 1 6528 1575 3382 1842
2 10 1834 873 1285 850
3 100 12177 1979 5965 1708
4 1000 8583 2978 4016 2744
5 10000 22459 3296 7950 2436
6 100008586 693 354 585 573
boxplot(log2(a[,-1]),outline = F)
不同计算方法基因长度
gene_id 转换到symbol看下
library(org.Hs.eg.db)
s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
aa =merge(s2g,a,by='gene_id') #19601个基因，比id转换前少了25个基因
head(aa)
gene_id symbol exon CDS max average
1 1 A1BG 6528 1575 3382 1842
2 10 NAT2 1834 873 1285 850
3 100 ADA 12177 1979 5965 1708
4 1000 CDH2 8583 2978 4016 2744
5 10000 AKT3 22459 3296 7950 2436
6 100008586 GAGE12F 693 354 585 573
1.不同版本的注释文件，得到的长度不太一样，我的结果是来自TxDb.Hsapiens.UCSC.hg38.knownGene，小洁的结果是Homo_sapiens.GRCh38.103.chr.gtf.gz。 可以看到exon 和cds 的计算方式差别更小，最长转录本和平均转录本的计算方式差的比较多。
我的结果
小洁的结果2. boxplot的结果和小洁的也不太一样？ 我的看起来exon计算结果基因长度最长，小洁的是按最长转录本的方式，基因长度最长.
