从gtf文件中计算每个基因exon总长python脚本，给有需要的人

大裤衩

代码如下：

[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])}')

Wuii

厉害，只是有重叠的exon也要处理

基迪奥-李泽标

点赞~

大裤衩

Wuii 发表于 2019.6.4 07:59
厉害，只是有重叠的exon也要处理

你的意思是多个exon重叠成一个?

Wuii

大裤衩 发表于 2019.6.4 11:27
你的意思是多个exon重叠成一个?

嗯，我只是大概看下代码，你的脚本只考虑了 每个基因一个转录本的情况？

大裤衩

Wuii 发表于 2019.6.4 16:40
嗯，我只是大概看下代码，你的脚本只考虑了 每个基因一个转录本的情况？ ...

不是啊，所有的exon啊，代码第20行就是对这个处理的

Wuii

大裤衩 发表于 2019.6.4 17:46
不是啊，所有的exon啊，代码第20行就是对这个处理的

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

大裤衩

Wuii 发表于 2019.6.5 22:08
恩，可能我没有说清楚。假设 一个gtf文件中，
有gene_id A , 对应了 两个 transcript_id 1 2 ，分别有9和 ...

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

Wuii

大裤衩 发表于 2019.6.5 23:25
我算明白你的意思了，一个基因多个转录本若有相同的exon，我这个计算方法确实是计算重复了。
哈哈，些太 ...

恩，挺好，只是提醒。加油整吧。

大裤衩

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])}')

Wuii

大裤衩 发表于 2019.6.6 10:16
嗯，改了下，效率及其低下
[mw_shl_code=python,true]import re
import sys

恩，可以继续修改。需要考虑，转录本A 的 Exon 和 转录本B 的 Exon 是存在Overlap的情况