Python脚本之从gtf文件中提取splice sites

代码如下:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# extract_splice_sites.py
# @author Aaron Liu
# @description 从gtf文件中提取剪切位点(参考hsiat2_extract_splice_sites.py)
# @created 2020-09-04T19:26:12.427Z+08:00
# @last-modified 2020-09-05T17:33:50.935Z+08:00

from argparse import ArgumentParser, FileType

def extract_splice_sites(gtf_file):
    trans = {}
    for line in gtf_file:
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        if '#' in line:
            line = line.strip('#')[0].strip()

        chr, source, feature, start, end, score, strand, \
            frame, values = line.split('\t')
        start, end = int(start), int(end)

        if feature != 'exon' or start >= end:
            continue

        value_attr = {}
        for attr in values.split(';'):
            gene, _, transcript = attr.strip().partition(' ')
            value_attr[gene] = transcript.strip('"')

        if "gene_id" not in value_attr \
            or "transcript_id" not in value_attr:
            continue

        transcript_id = value_attr['transcript_id']
        gene_id = value_attr['gene_id']
        if transcript_id not in trans:
            trans[transcript_id] = [gene_id, chr, strand, [[start, end]]]
        else:
            trans[transcript_id][3].append([start, end])

    for transid, [gene_id, chr, strand, exon] in trans.items():
        exon.sort()
        tem_exon = [exon[0]]
        for i in range(1, len(exon)):
            if exon[i][0] - tem_exon[-1][1] <= 5:
                tem_exon[-1][1] = exon[i][1]
            else:
                tem_exon.append(exon[i])
        trans[transid] = [gene_id, chr, strand, tem_exon]

    junction = set()
    for gene, chrom, strand, exon in trans.values():
        exon.sort()
        for i in range(1, len(exon)):
            junction.add(
                (gene, chrom, exon[i - 1][1] + 1, exon[i][0] - 1, strand))

    junction = sorted(junction)
    return junction


if __name__ == "__main__":
    parse = ArgumentParser(description='Extract splice sites from a GTF file')
    parse.add_argument('-i',
                       '--gtf_file',
                       required=True,
                       nargs='?',
                       type=FileType('r'),
                       help="a gtf file is required")
    parse.add_argument('-o',
                       '--out_file',
                       required=True,
                       nargs='?',
                       type=FileType('w'),
                       help='output file is required')

    args = parse.parse_args()

    splicesites = extract_splice_sites(args.gtf_file)
    for row in splicesites:
        args.out_file.write(row[0] + '\t' + row[1] + '\t' + str(row[2]) +
                            '\t' + str(row[3]) + '\t' + row[4] + '\n')

实例:

./extract_splice_sites.py -i merge.gtf -o ss.txt

输出:

...
MSTRG.10013 chr04   2534745 2535071 +
MSTRG.10013 chr04   2535164 2535818 +
MSTRG.10013 chr04   2535988 2536609 +
MSTRG.10013 chr04   2536732 2537314 +
MSTRG.10013 chr04   2536732 2537381 +
MSTRG.10013 chr04   2537444 2537540 +
MSTRG.10013 chr04   2537444 2537551 +
MSTRG.10013 chr04   2537651 2537730 +
MSTRG.10013 chr04   2537789 2537870 +
...
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,761评论 5 460
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,953评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,998评论 0 320
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,248评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,130评论 4 356
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,145评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,550评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,236评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,510评论 1 291
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,601评论 2 310
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,376评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,247评论 3 313
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,613评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,911评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,191评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,532评论 2 342
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,739评论 2 335