搜索
查看: 2274|回复: 3

生信编程实战代码注释【第三题-hg38每条染色体。。】

[复制链接]

10

主题

11

帖子

400

积分

中级会员

Rank: 3Rank: 3

积分
400
发表于 2017-7-9 09:23:38 | 显示全部楼层 |阅读模式
[Python] 纯文本查看 复制代码
'''
首先声明【代码来自生信编程直播第三题-hg38每条染色体基因转录本的分布-python-东】
我是新手,只是学习了,做了代码注释,[b][color=#ff0000]如果注释有错,欢迎小伙伴指出,以防祸害别人[/color][/b]
本节课目的:对这个文件[url]ftp://ftp.ensembl.org/pub/releas[/url] ... RCh38.87.chr.gtf.gz  进行如下统计:
染色体上基因数量分布
基因长度分布情况
基因的转录本数量分布
转录本的外显子数量统计
外显子坐标统计
#注释者:默存,20170708
'''
# -*- coding: UTF-8 -*-  
import sys
import re 
args = sys.argv
class Genome_info:# Genome_info类(基类)为后续共有属性,单独列出来,后续都继承
    def __init__(self): #python中定义一个类方法时规定:第一个参数必须为self,__init__方法:属于python语言类的构造函数,一个类只能有一个__init__方法
        self.chr = ""#属于每个对象级别的变量,在调用的时候一定要带上self表明属于当前对象。
        self.start = 0
        self.end = 0
 
 #继承语法   class 派生本类名(基类名) eg:class Gene(Genome_info)
class Gene(Genome_info):#Gene类继承Genome_info类的属性
    def __init__(self):#
        Genome_info.__init__(self)#在调用基类的方法时,需要加上基类的类名前缀,且需要带上self参数,变量继承时,
        self.orientation = ""
        self.id = ""

class Transcript(Genome_info):
    def __init__(self):
        Genome_info.__init__(self)
        self.id = ""
        self.parent = ""

class Exon(Genome_info):
    def __init__(self):
        Genome_info.__init__(self)
        self.parent = ""
 
def main(args):#主函数统计gtf中我们所需要信息
    """
一个输入参数:
    第一个参数为物种gtf文件
    :return:
    """
    list_chr = []
    list_gene = {}#含有gene_id
    list_transcript = {}#含有transcript_id
    list_exon = []
    
    with open(args[1]) as fp_gtf:#读入gtf文件,args[1]代表gtf文件
        for line in fp_gtf:#遍历gtf
            if line.startswith("#"):#过滤掉#开头的注释
                continue            
            lines = line.strip("\n").split("\t")#strip() 移除行末换行符\n,使用split()"\t"分割
            chr = lines[0]#染色体号
            type = lines[2] #类型:gene  transcript  exon
            start = int(lines[3])#转化下数据类型
            end = int(lines[4])
            orientation = lines[6]
            attr = lines[8]
            if not re.search(r'protein_coding', attr):#正则寻找编码蛋白的序列类型
                continue 
            if not chr in list_chr :
                list_chr.append(chr)#如果list_chr中没有,append()添加入
 
            if type == "gene":
                gene = Gene()#初始化一个Gene类
                id = re.search(r'gene_id "([^;]+)";?', attr).group(1)#group(1)代表([^;]+)
                gene.chr = chr 
                gene.start = start
                gene.end = end
                gene.id = id
                gene.orientation = orientation
                list_gene[id] = gene#将数值添加到list_gene字典中
                
            elif type == "transcript":
                transcript = Transcript()
                id = re.search(r'transcript_id "([^;]+)";?', attr).group(1)#transcript_id
                parent = re.search(r'gene_id "([^;]+)";?', attr).group(1)#gene_id作为parent
                if not parent in list_gene:#排除没有所属基因的转录本
                    continue
                transcript.chr = chr
                transcript.start = start
                transcript.end = end
                transcript.id = id
                transcript.parent = parent
                list_transcript[id] = transcript#将数值添加到list_transcript字典中
 
            elif type == "exon":
                exon = Exon()
                parent = re.search(r'transcript_id "([^;]+)";?', attr).group(1)#parent是transcript
                if not parent in list_transcript:
                    continue
                exon.chr = chr
                exon.start = start
                exon.end = end
                exon.parent = parent 
                list_exon.append(exon)#将数值添加到list_exon列表中
 #取出我们要的信息
    chr_gene(list_gene)#染色体上基因数量分布
    gene_len(list_gene)#基因长度分布情况
    gene_transcript(list_transcript)#基因的转录本数量分布
    transcript_exon(list_exon)#转录本的外显子数量统计
    exon_pos(list_exon)#外显子坐标统计
 
 #取出我们需要的信息
def chr_gene(list_gene):
    """
    染色体上基因数量分布
    :param list_gene:
    :return:
    """
 
    print("染色体上基因数量分布")
    #统计
    count_gene = {}#计数器{info.chr,info}
    for info in list_gene.values():#遍历list_gene,只取values
        chr = info.chr#list_gene中的value info对应的染色体号
        if chr in count_gene:
            count_gene[info.chr] += 1 
        else:
            count_gene[info.chr] = 1
            print(count_gene)
            
     #输出
    with open("chr_gene.txt", 'w') as fp_out:
        for chr, num in count_gene.items():
            print("\t".join([chr, str(num)]) + "\n")#染色体号和基因之间用"\t"连接
            fp_out.write("\t".join([chr, str(num)]) + "\n")

def gene_len(list_gene):#遍历字典list_gene
    """
    基因长度分布情况
 
    :param list_gene:
    :return:
    """
 
    print ("基因长度分布情况")
    with open("gene_len.txt", 'w') as fp_out:
        for gene_id, info in list_gene.items():
            len = info.end - info.start + 1
            fp_out.write("\t".join([gene_id, str(len)]) + "\n")
            print("\t".join([gene_id, str(len)]) + "\n")
 
def gene_transcript(list_transcript):
    """
    基因的转录本数量分布
 
    :param list_transcript:
    :return:
    """
 
    print("基因的转录本数量分布")
    count_transcript = {}
    for info in list_transcript.values():
        gene_id = info.parent
        if gene_id in count_transcript:
            count_transcript[gene_id] += 1
        else:
            count_transcript[gene_id] = 1
    with open("gene_transcript.txt", 'w') as fp_out:
        for gene_id, num in count_transcript.items():
            print("\t".join([gene_id, str(num)]) + "\n")
            fp_out.write("\t".join([gene_id, str(num)]) + "\n")
 
 
def transcript_exon(list_exon):
    """
    转录本的外显子数量统计
 
    :param list_exon:
    :return:
    """
 
    print("转录本的外显子数量统计")
    count_exon = {}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id in count_exon:
            count_exon[transcript_id] += 1
        else:
            count_exon[transcript_id] = 1
    with open("transcript_exon.txt", 'w') as fp_out:
        for transcript_id, num in count_exon.items():
            print("\t".join([transcript_id, str(num)]) + "\n")
            fp_out.write("\t".join([transcript_id, str(num)]) + "\n")
 
def exon_pos(list_exon):
    """
    外显子坐标统计
 
    :param list_exon:
    :return:
    """
 
    print("外显子坐标统计")
    count_exon = {}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id in count_exon:
            count_exon[transcript_id] += ",%s-%s" % (str(exon.start), str(exon.end))#python格式化输出,不是第一次,加,
        else:
            count_exon[transcript_id] = "%s-%s" % (str(exon.start), str(exon.end))#第一次,没,
    with open("exon_pos.txt", 'w') as fp_out:
        for transcript_id, pos in count_exon.items():
            print("\t".join([transcript_id, pos]) + "\n")
            fp_out.write("\t".join([transcript_id, pos]) + "\n") 
 
if __name__ == "__main__":#每个模块都有自己的名称__name__,如果模块直接执行,__name__="__main__";此语句用于模块测试中使用。
    main(args)
   
   #脚本执行:python gtf_count.py Homo_sapiens.GRCh38.87.chr.gtf
'''         
这一节补充学习的:
推荐使用git bash for windows  [url]https://git-for-windows.github.io/[/url]
python中的继承使用:
[url]http://www.cnblogs.com/Joans/archive/2012/11/09/2757368.html[/url]
[url]http://www.cnblogs.com/Joans/archive/2012/11/06/2749993.html[/url]
 '''





上一篇:构建nr或者nt数据库的子库
下一篇:Rsamtools
回复

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-7-20 05:15:15 | 显示全部楼层
请问.group()那个是什么语法啊
回复 支持 反对

使用道具 举报

10

主题

11

帖子

400

积分

中级会员

Rank: 3Rank: 3

积分
400
 楼主| 发表于 2017-7-20 07:34:05 | 显示全部楼层
朝曦曦 发表于 2017-7-20 05:15
请问.group()那个是什么语法啊

group()这个是正则表达式里面的,group[1]代表前面()里面正则匹配的东西,如果多次匹配的话,还有group(2).。。group(n)
回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-8-1 14:26:22 | 显示全部楼层
默存 发表于 2017-7-20 07:34
group()这个是正则表达式里面的,group[1]代表前面()里面正则匹配的东西,如果多次匹配的话,还有gro ...

谢谢楼主!
我还有一个问题一直没搞懂,就是main(args)一直理解不透彻。我也上网查了,似乎是能看懂,自己不会用。楼主方便稍微说两句解释一下吗?
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:31 , Processed in 0.036259 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.