python 学习之 featureCounts 软件的基因长度是怎么算的?
你需要一颗沉静的心
1引言
featureCounts 软件问题
RNA-seq 流程上游 最后一步需要对基因或者转录本进行定量, featureCounts 这款软件在对基因进行定量时,我们需要指定 -t exon 参数, 在定量的结果文件中会有 基因长度 这一列:
但是没有仔细去研究这个 基因长度是怎么计算出来的, 去看这个软件官方> 文档也没有仔细说明,那就只有自己去探索了。 是不是选的最长的转录本?
提取 GTF 文件转录本长度
想要解答疑问,我们首先需要 获取 gtf 注释文件每个基因的每个转录本的长度信息。
此外不同注释文件的记录信息也会有所差别,尤其是不同数据库的注释文件, UCSC 数据库的结果比较保守, ENSEMBL 和 GENCODE 数据库的信息是最全的,但是准确性可能会低一些。
注意:
UCSC和ENSEMBL,GENCODE的 gtf 注释文件有一些不同。
所以我会写两个 python 脚本来提取不同数据库的 gtf 注释文件的基因转录本长度信息。这里我下载了 UCSC 的 hg38 和 gencode GRCh38 的 gtf 注释文件
hg38.ncbiRefSeq.gtf,gencode.v38.annotation_human.gtf。
2提取不同数据库 gtf 的基因转录本长度
计算的基本思路就是 累加每个基因的每个转录本 exon 长度并保存信息到文件里。
提取 UCSC 数据库转录本长度
outinfo = open('ucsc_hg38_gtf_info.txt','w')
# 打开测试 gtf 文件
with open('hg38.ncbiRefSeq.gtf','r') as gtf:
# 信息保存在字典里
info = {}
for line in gtf:
# 跳过注释行
if line.startswith('#'):
continue
# 分割
fields = line.split()
# 类型
type = fields[2]
if type == 'exon':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
trans_id = fields[11].replace('"','').replace(';','')
# biotype = fields[23].replace('"','').replace(';','')
# 连接名称
key = gene_name + 't' + gene_id + 't' + trans_id
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4])
length = end - start + 1
# 累计求和
info.setdefault(key,0)
info[key] += length
# 导出保存
for key,val in info.items():
newline = key + 't' + str(val) + 'n'
outinfo.write(newline)
# 关闭文件
outinfo.close()
查看结果,分别为 gene_name,gene_id,trans_id 和 长度:
上图可以看到有些基因只有一个转录本,有些有多个转录本。
提取 GENCODE 数据库转录本长度
outinfo = open('gencode_hg38_gtf_info.txt','w')
# 打开测试 gtf 文件
with open('gencode.v38.annotation_human.gtf','r') as gtf:
# 信息保存在字典里
info = {}
for line in gtf:
# 跳过注释行
if line.startswith('#'):
continue
# 分割
fields = line.split()
# 类型
type = fields[2]
if type == 'exon':
# 名称
gene_name = fields[15].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
trans_id = fields[11].replace('"','').replace(';','')
gene_type = fields[13].replace('"','').replace(';','')
# 连接名称
key = gene_name + 't' + gene_id + 't' + trans_id + 't' + gene_type
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4])
length = end - start + 1
# 累计求和
info.setdefault(key,0)
info[key] += length
# 导出保存
for key,val in info.items():
newline = key + 't' + str(val) + 'n'
outinfo.write(newline)
# 关闭文件
outinfo.close()
查看结果,分别为 gene_name,gene_id,trans_id,gene_type 和 长度:
3结果对比
我们随便找一个 编码基因 ACTB 比较:
# UCSC
ACTB ACTB NM_001101.5 1812
# GENCODE
ACTB ENSG00000075624.17 ENST00000674681.1 protein_coding 2554
ACTB ENSG00000075624.17 ENST00000642480.2 protein_coding 2021
ACTB ENSG00000075624.17 ENST00000676397.1 protein_coding 1804
ACTB ENSG00000075624.17 ENST00000676319.1 protein_coding 811
ACTB ENSG00000075624.17 ENST00000676189.1 protein_coding 1815
ACTB ENSG00000075624.17 ENST00000473257.3 protein_coding 1687
# 相同长度的
ACTB ENSG00000075624.17 ENST00000646664.1 protein_coding 1812
#
ACTB ENSG00000075624.17 ENST00000464611.1 protein_coding 715
ACTB ENSG00000075624.17 ENST00000477812.2 protein_coding 2271
ACTB ENSG00000075624.17 ENST00000675515.1 protein_coding 1919
ACTB ENSG00000075624.17 ENST00000425660.5 protein_coding 1845
ACTB ENSG00000075624.17 ENST00000462494.5 protein_coding 2245
ACTB ENSG00000075624.17 ENST00000493945.6 protein_coding 1495
ACTB ENSG00000075624.17 ENST00000432588.6 protein_coding 1300
ACTB ENSG00000075624.17 ENST00000645576.1 protein_coding 1009
ACTB ENSG00000075624.17 ENST00000647275.1 protein_coding 556
ACTB ENSG00000075624.17 ENST00000484841.6 protein_coding 561
ACTB ENSG00000075624.17 ENST00000645025.1 protein_coding 572
ACTB ENSG00000075624.17 ENST00000480301.1 protein_coding 657
ACTB ENSG00000075624.17 ENST00000443528.5 protein_coding 569
ACTB ENSG00000075624.17 ENST00000417101.2 protein_coding 472
ACTB ENSG00000075624.17 ENST00000414620.1 protein_coding 561
ACTB ENSG00000075624.17 ENST00000646584.1 protein_coding 425
可以看到 UCSC 只有一个转录本,而 GENCODE 多达 22 个转录本!
我们看看 featureCounts 软件计算出来的长度却并不一样,是 4405:
软件定量使用的 GTF 注释文件就是 gencode 的, 如果是选择最长转录本的话应该是 2554!所以猜想错误。
那么 featureCounts 软件是如何选择计算一个基因的长度的呢?
4计算 featureCounts 软件基因长度
我们仔细观察定量结果文件会发现每个基因都有其 start 和 end 两列信息:
代表了 GTF 文件里这个 基因所有外显子的起始位置和终止位置的坐标, 经过摸索终于搞清了长度是如何计算出来的,示意图如下:
原理就是 取所有非冗余的外显子 或者是 所有外显子的并集 !
我们可以提取 featureCounts 结果文件的信息自己计算验证一下:
library(tidyverse)
# 读取结果
gene_info <- read.delim('featureCounts_gene_Length.txt',header = T)
# 查看内容
head(gene_info,3)
Geneid
1 ENSG00000223972.5
2 ENSG00000227232.5
3 ENSG00000278267.1
Start
1 11869;12010;12179;12613;12613;12975;13221;13221;13453
2 14404;15005;15796;16607;16858;17233;17606;17915;18268;24738;29534
3 17369
End Length
1 12227;12057;12227;12721;12697;13052;14409;13374;13670 1735
2 14501;15038;15947;16765;17055;17368;17742;18061;18366;24891;29570 1351
3 17436 68
gene_name
1 DDX11L1
2 WASH7P
3 MIR6859-1
写个循环计算每个基因的所有外显子的并集及长度, 注意: 基因名不能作为循环的元素:
# 基因名不唯一
gene_name = gene_info$gene_name
length(gene_name)
[1] 60651
length(unique(gene_info$gene_name))
[1] 59386
# 基因id唯一
gene_id = gene_info$Geneid
length(gene_id)
[1] 60651
length(unique(gene_info$Geneid))
[1] 60651
循环计算:
lapply(gene_id[1:50], function(x){
# 筛选基因
tmp <- gene_info %>% filter(Geneid == x)
# 去除exon起始和终止位置
start <- unlist(strsplit(tmp$Start,split = ';')) %>% as.numeric()
end <- unlist(strsplit(tmp$End,split = ';')) %>% as.numeric()
# 转为数据框
tmp1 <- data.frame(start,end)
# 计算所有外显子的并集
fields = c()
# loop
for (j in 1:nrow(tmp1)) {
fields <- union(fields,tmp1[j,1]:tmp1[j,2])
}
GneneLenth = length(fields)
# 合并结果
gene_name = tmp$gene_name
res <- data.frame(gene_id = x,gene_name = gene_name,Length = GneneLenth)
return(res)
}) %>% Reduce('rbind',.) %>%
as.data.frame() -> final_res
我们把自己计算的结果和 featureCounts 软件的 Length 列 对比:
可以看到结果是一模一样的。
5结尾
-
定量的话,对于 GTF 文件的选择其实大家也不必纠结,基于以上探索结果的话, 如果用
UCSC的注释文件来定量的话, featureCounts 软件计算的基因长度是会比ENSEMBL/GENCODE数据库的小的。那么后面计算出来的fpkm/rpkm/tpm也会有差异。 -
用哪个都没有对错, 根据自己选择, 我觉得其实都没必要对基因长度进行标准化矫正,直接用 rpm 即可, 因为我们往往是比较同一个基因在不同样本之间的变化,而同一个基因在不同样本之间的长度肯定是一样的。
-
如果你想 比较同一个样本里不同基因的表达相对高低, 那么这就需要对基因长度进行标准化了。
-
有时候我们经常只有 count 矩阵文件, 但没有基因长度文件, 后面需要计算
fpkm/rpkm/tpm的话就比较麻烦,后续我会按照 featureCounts 的计算原理写一个脚本来从 GTF 文件获取基因的长度信息,敬请期待…
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
欢迎小伙伴留言评论！
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀python 学习之 python 里也能用 dplyr?
◀python 学习之 提取 Ensembl,Gencode 和 Ucsc 基因 TSS 位点
◀python 学习之 R and Python: 循环函数
◀...
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
