需求
tpm、fpkm这样的数据，计算公式里都有“基因有效长度”，有些地方直接叫基因长度。之前计算这些，都是采用最快的方式—找豆豆。今天想要把这部分逐渐加到课程里，我就必须要研究一下基因长度这个东东了。翻遍全网，找到了关于基因长度最优秀的描述，当属曾老板写的：基因长度知多少（给曾老板手动点赞👍）。
其中的基因长度依据的参考基因组是hg19，我想要的是从gtf文件开始计算，在曾老板的代码基础上改吧一下就行。因为可能要涉及到其他物种，或者其他版本，这个需求还是比较普遍的。
1.从gtf文件开始
gtf文件的下载，主要是GENECODE和ENSEMBL这两个网站，搜一哈，点进去就能找到最新版本。（目前最新版本是103，以后肯定是继续更新的）如果是其他模式生物，可以自己查查参考基因组注释文件从哪里下载。
library(rtracklayer)
gtf = rtracklayer::import("Homo_sapiens.GRCh38.103.chr.gtf.gz")
class(gtf)
## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"
gtf = as.data.frame(gtf);dim(gtf)
## [1] 3073646 27
table(gtf$type)
##
## gene transcript exon CDS start_codon
## 60607 234326 1460128 806304 91602
## stop_codon five_prime_utr three_prime_utr Selenocysteine
## 84603 160663 175296 117
首先把gtf读取进R语言，会成为一个Grange对象。不了解Grange无所谓的，反正可以转换成数据框。这里面有约300万行（不同版本行数不一样，逐渐进化中）。
其中基因占了6w，转录本、外显子、CDS都被单独罗列出来了。
引用曾老板，几种主流的基因长度计算方法：
挑选基因的最长转录本
选取多个转录本长度的平均值
非冗余外显子(EXON)长度之和
非冗余 CDS（Coding DNA Sequence） 长度之和
今天就玩这个啦！挨个试试。
1.挑选基因的最长转录本
library(tidyverse)
tra = gtf[gtf$type=="transcript",
c("start","end","gene_name")]
length(unique(tra$gene_name))
## [1] 59409
glt = mutate(tra, length = end - start) %>%
arrange(desc(length)) %>%
filter(!duplicated(gene_name)) %>%
select(c(3,4))
head(glt)
## gene_name length
## 1 RBFOX1 2471656
## 2 CNTNAP2 2304197
## 3 DLG2 2172171
## 4 DMD 2092327
## 5 PTPRD 2084571
## 6 CSMD1 2059553
2.选取多个转录本长度的平均值
gltm = mutate(tra, length = end - start) %>%
group_by(gene_name) %>%
summarise(length = mean(length))
head(gltm)
## # A tibble: 6 x 2
## gene_name length
## <chr> <dbl>
## 1 5_8S_rRNA 151.
## 2 5S_rRNA 112.
## 3 7SK 266.
## 4 A1BG 4146
## 5 A1BG-AS1 6124.
## 6 A1CF 68834.
3.非冗余外显子(EXON)长度之和
exon = gtf[gtf$type=="exon",
c("start","end","gene_name")]
length(unique(exon$gene_name))
## [1] 59409
if(F){
gle = lapply(split(exon,exon$gene_name),function(x){
tmp=apply(x,1,function(y){
y[1]:y[2]
})
length(unique(unlist(tmp)))
})
gle=data.frame(gene_name=names(gle),
length=as.numeric(gle))
save(gle,file = "gle.Rdata")
}
load("gle.Rdata")
head(gle)
## gene_name length
## 1 5_8S_rRNA 911
## 2 5S_rRNA 1021
## 3 7SK 1870
## 4 A1BG 3999
## 5 A1BG-AS1 3374
## 6 A1CF 9603
4.非冗余 CDS长度之和
CDS = gtf[gtf$type=="CDS",
c("start","end","gene_name")]
length(unique(CDS$gene_name))
## [1] 20382
if(F){
glc = lapply(split(CDS,CDS$gene_name),function(x){
tmp=apply(x,1,function(y){
y[1]:y[2]
})
length(unique(unlist(tmp)))
})
glc=data.frame(gene_name=names(glc),
length=as.numeric(glc))
save(glc,file = "glc.Rdata")
}
load("glc.Rdata")
head(glc)
## gene_name length
## 1 A1BG 1572
## 2 A1CF 1905
## 3 A2M 4422
## 4 A2ML1 4502
## 5 A3GALT2 1020
## 6 A4GALT 4236
比较一下几种方法得到的基因长度
a = gle %>% inner_join(glc,"gene_name") %>%
inner_join(glt,"gene_name") %>%
inner_join(gltm,"gene_name")
colnames(a) = c("gene_name","exon","CDS","tra_max","tra_mean")
a[1:4,]
## gene_name exon CDS tra_max tra_mean
## 1 A1BG 3999 1572 8309 4146.00
## 2 A1CF 9603 1905 86266 68834.40
## 3 A2M 6318 4422 48211 11661.31
## 4 A2ML1 7156 4502 54166 17511.89
boxplot(log2(a[,-1]),outline = F)
接下来是一点补充讨论
其实上面列出来的几种方法，用哪个都可以的。我问豆豆,你说这里面哪个最好呢。豆豆说，他用非冗余外显子。好的，以后没什么特殊要求，默认的长度就它了。为什么呢？因为默认大(lao)佬(gong)说的都是对的。
有选择困难症的时候，最简单的办法是看看大佬是怎么做的。查查文献或者翻翻网上的帖子都行呗。
豆豆还找到了一篇关于基因长度的文章，还挺新的，这一篇里面选择了最长转录本长度作为基因长度。
https://www.frontiersin.org/articles/10.3389/fgene.2021.559998/full
对接count与fpkm、tpm的转换
豆豆写过了，我选择直接复制粘贴：https://www.jianshu.com/p/9dfb65e405e8
网友评论