查看: 16889|回复: 10

[其他] 从gtf文件中计算每个基因exon总长python脚本,给有需要的人

[复制链接]

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2018.1.3
在线时间
338 小时

热心会员活跃会员


发表于 2019.6.3 23:32:14 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
代码如下:
[Python] 纯文本查看 复制代码
import re
import sys

patG=re.compile(r'gene_id\s+"(.+?)";')
Gene=dict()

d=open(sys.argv[1],'r')
for each in d:
        if each.startswith("#"):
                continue
        line=each.strip().split()
        typ=line[2]
        if typ != "exon":
                continue
        start=line[3]
        end=line[4]
        each_len=int(end)-int(start)+1
        mgene=patG.search(each)
        geneName=mgene.groups()[0]
        Gene.setdefault(geneName,[]).append(each_len)


for each in Gene:
        print(f'{each}\t{sum(Gene[each])}')

评分

参与人数 1奥币 +10 贡献 +2 收起 理由
基迪奥-李泽标 + 10 + 2 楼主V5!

查看全部评分

新的一天加油!
回复

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
15
注册时间
2016.4.20
在线时间
466 小时

突出贡献优秀版主论坛元老


发表于 2019.6.4 07:59:01 | 显示全部楼层
厉害,只是有重叠的exon也要处理
回复 支持 反对

使用道具 举报

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
402
注册时间
2018.4.19
在线时间
896 小时

推广达人宣传达人


发表于 2019.6.4 08:48:10 | 显示全部楼层
点赞~
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2018.1.3
在线时间
338 小时

热心会员活跃会员


 楼主| 发表于 2019.6.4 11:27:18 | 显示全部楼层
Wuii 发表于 2019.6.4 07:59
厉害,只是有重叠的exon也要处理

你的意思是多个exon重叠成一个?
新的一天加油!
回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
15
注册时间
2016.4.20
在线时间
466 小时

突出贡献优秀版主论坛元老


发表于 2019.6.4 16:40:32 | 显示全部楼层
大裤衩 发表于 2019.6.4 11:27
你的意思是多个exon重叠成一个?

嗯,我只是大概看下代码,你的脚本只考虑了 每个基因一个转录本的情况?
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2018.1.3
在线时间
338 小时

热心会员活跃会员


 楼主| 发表于 2019.6.4 17:46:01 | 显示全部楼层
Wuii 发表于 2019.6.4 16:40
嗯,我只是大概看下代码,你的脚本只考虑了 每个基因一个转录本的情况? ...

不是啊,所有的exon啊,代码第20行就是对这个处理的
新的一天加油!
回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
15
注册时间
2016.4.20
在线时间
466 小时

突出贡献优秀版主论坛元老


发表于 2019.6.5 22:08:50 | 显示全部楼层
大裤衩 发表于 2019.6.4 17:46
不是啊,所有的exon啊,代码第20行就是对这个处理的

恩,可能我没有说清楚。假设 一个gtf文件中,
有gene_id A , 对应了 两个 transcript_id 1 2 ,分别有9和10个exon,每个长度为10bp,那么按照你的脚本逻辑,计算结果是190bp,而 A 事实上长度应该是 100bp。
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2018.1.3
在线时间
338 小时

热心会员活跃会员


 楼主| 发表于 2019.6.5 23:25:48 | 显示全部楼层
Wuii 发表于 2019.6.5 22:08
恩,可能我没有说清楚。假设 一个gtf文件中,
有gene_id A , 对应了 两个 transcript_id 1 2 ,分别有9和 ...

我算明白你的意思了,一个基因多个转录本若有相同的exon,我这个计算方法确实是计算重复了。
哈哈,些太急促,没考虑到这一点儿,嗯,归根到底还是没把gtf文件记录理解透彻,
新的一天加油!
回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
15
注册时间
2016.4.20
在线时间
466 小时

突出贡献优秀版主论坛元老


发表于 2019.6.6 09:02:29 | 显示全部楼层
大裤衩 发表于 2019.6.5 23:25
我算明白你的意思了,一个基因多个转录本若有相同的exon,我这个计算方法确实是计算重复了。
哈哈,些太 ...

恩,挺好,只是提醒。加油整吧。
回复 支持 反对

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2018.1.3
在线时间
338 小时

热心会员活跃会员


 楼主| 发表于 2019.6.6 10:16:04 | 显示全部楼层
Wuii 发表于 2019.6.6 09:02
恩,挺好,只是提醒。加油整吧。

嗯,改了下,效率及其低下
[Python] 纯文本查看 复制代码
import re
import sys

patG=re.compile(r'gene_id\s+"(.+?)";')
Gene=dict()
uniq=[]
d=open(sys.argv[1],'r')
for each in d:
        if each.startswith("#"):
                continue
        line=each.strip().split()
        typ=line[2]
        posinfo=line[3]+"-"+line[4]

        if typ != "exon":
                continue
        start=line[3]
        end=line[4]
        each_len=int(end)-int(start)+1
        mgene=patG.search(each)
        geneName=mgene.groups()[0]
        if posinfo not in uniq:
                Gene.setdefault(geneName,[]).append(each_len)
                uniq.append(posinfo)


for each in Gene:
        print(f'{each}\t{sum(Gene[each])}')
新的一天加油!
回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
15
注册时间
2016.4.20
在线时间
466 小时

突出贡献优秀版主论坛元老


发表于 2019.6.6 23:07:41 | 显示全部楼层
大裤衩 发表于 2019.6.6 10:16
嗯,改了下,效率及其低下
[mw_shl_code=python,true]import re
import sys

恩,可以继续修改。需要考虑,转录本A 的 Exon 和 转录本B 的 Exon 是存在Overlap的情况
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表