美文网首页ggplot集锦
matplotlib:根据gtf文件画出基因全部的转录本结构

matplotlib:根据gtf文件画出基因全部的转录本结构

作者: 所以suoyi | 来源:发表于2021-08-02 14:58 被阅读0次
    import numpy as np
    from matplotlib import pyplot as plt
    
    面向对象的界面:显式创建图形和轴,并在它们上调用方法
    x = np.linspace(0, 2, 100)
    fig, ax = plt.subplots()
    ax.set_xlabel('x label')  # x轴
    ax.set_ylabel('y label')  # y轴
    ax.set_title("Simple Plot")  # 图名
    ax.plot(x, x, label='linear')  # 画x
    ax.plot(x, x**2, label='quadratic')  # 画x方
    ax.plot(x, x**3, label='cubic')  # 画x立方
    ax.legend()  # 画图例
    
    pyplot 界面:依靠pyplot自动创建和管理图形和轴,并使用pyplot函数进行绘图
    x = np.linspace(0, 2, 100)
    plt.xlabel('x label')
    plt.ylabel('y label')
    plt.title("Simple Plot")
    plt.plot(x, x, label='linear')  # Plot some data on the (implicit) axes.
    plt.plot(x, x**2, label='quadratic')  # etc.
    plt.plot(x, x**3, label='cubic')
    plt.legend()
    
    image.png
    from matplotlib import pyplot as plt
    
    fig = plt.figure(1)  # 新建一个名叫 Figure1的画图窗口
    fig.patch.set_alpha(1)  # 设置不透明度,默认为1,完全不透明
    fig.patch.set_facecolor('w')  # 自定义背景色 "w" 白色
    

    根据gtf文件画出某个基因全部的转录本结构

    http://www.biotrainee.com/thread-624-1-1.html 此片段改编自生信技能树论坛生信编程直播第五题:根据GTF画基因的多个转录本结构 发表于 2017-10-4 00:41:45 的 源氏 的回答

    # https://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz
    def find_target_data(gtf, gene):
        from collections import defaultdict
        gene_transInfo = defaultdict(list)
        with open(gtf, 'r') as f:
            """1    ensembl_havana  gene    65419   71585   .       +       .       gene_id "ENSG00000186092"; 
                    gene_version "6"; gene_name "OR4F5"; gene_source "ensembl_havana"; gene_biotype "protein_coding";"""
            for line in f:
                if f'gene_name "{gene.upper()}"' in line:
                    line_spt = line.strip().split('\t')
                    try:
                        chr, db, record, start, end, score, strand, phase, info = line_spt
                        gene_transInfo['start'].append(start)
                        gene_transInfo['end'].append(end)
                        gene_transInfo['record'].append(record)
                    except:
                        pass
        if not gene_transInfo:
            print('\n\n There is some wrong with your gene name!\n')
            raise NameError('your gene is not exit')
    
        return gene_transInfo
    
    
    def draw_gene_structure(gene, gene_transInfo,  png_path):
        from matplotlib import pyplot as plt
        from matplotlib.ticker import MultipleLocator
        color_dic = {'exon': '#00896C', 'CDS': '#78552B', 'start_codon': '#CB1B45', 'stop_codon': 'black',
                     'five_prime_utr': '#F19483', 'three_prime_utr': '#2EA9DF'}
        linewith_dic = {'exon': 8.0, 'CDS': 6.0, 'start_codon': 8.0, 'stop_codon': 8.0, 'five_prime_utr': 4.0, 'three_prime_utr': 4.0}
    
        gene_start = min(map(int, gene_transInfo['start'])) - 500   # 两边边扩大些
        gene_end = max(map(int, gene_transInfo['end'])) + 500
    
        fig = plt.figure(figsize=(12, 6))  # 建个画板,画板大小
        ax = fig.add_subplot()  # 画板上建个画布,这里可以建多个画布
        ax.set_xlim(int(gene_start) - 500, int(gene_end) + 500)  # 画x轴
        ax.ticklabel_format(useOffset=False, style='plain')  # x,y轴禁用科学记数法
    
        t = 0
        record = gene_transInfo['record']
        for i in range(len(record)):
            if record[i] == 'transcript':
                t += 1
                ax.plot([int(gene_transInfo['start'][i]), int(gene_transInfo['end'][i])], [t, t], color="black")
            elif record[i] == 'gene':
                pass
            else:
                ax.plot([int(gene_transInfo['start'][i]), int(gene_transInfo['end'][i])], [t, t], color=color_dic[record[i]], linewidth=linewith_dic[record[i]])
    
        ax.set_title(f"the transcripts of {gene} gene")
        ymajorLocator = MultipleLocator(1)  # y轴主间距设为1,事实上按数量怎么顺眼怎么来
        ax.yaxis.set_major_locator(ymajorLocator)  #
        ax.set_ylim(0, t + 3)  # +3 单纯是让y轴长点
    
        # 添加图例
        import matplotlib.patches as mpatches
        legend_lt = []
        for region in color_dic:
            legend = mpatches.Patch(color=color_dic[region], label=region, linewidth=linewith_dic[region])
            legend_lt.append(legend)
        ax.legend(handles=legend_lt)
        #  这个图例画的有点丑 可以再调整调整
    
        fig.savefig(png_path, dpi=150)
        plt.show()
    
    
    gtf_path = r'D:\database\hg38\Homo_sapiens.GRCh38.100.chr.gtf'
    gene = 'ATP7B'
    gene_transInfo = find_target_data(gtf_path, gene)
    png_path = f'D:\result\{gene}.png'
    draw_gene_structure(gene, gene_transInfo, png_path)
    
    
    image.png

    调整坐标轴刻度间隔 (本片段为CSDN博主「南辰以北」的原创文章)

    from matplotlib.ticker import MultipleLocator, FormatStrFormatter
    
    xmajorLocator   = MultipleLocator(1)   x轴主刻度
    ax.xaxis.set_major_locator(xmajorLocator)
    
    ymajorLocator   = MultipleLocator(1)   y轴主刻度
    ax.yaxis.set_major_locator(ymajorLocator)
    
    xminorLocator   = MultipleLocator(0.25)   x轴副刻度
    ax.xaxis.set_minor_locator(xminorLocator)
    
    yminorLocator   = MultipleLocator(0.25)   y轴副刻度
    ax.yaxis.set_minor_locator(yminorLocator)
    ————————————————
    版权声明:本文为CSDN博主「南辰以北」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    原文链接:https://blog.csdn.net/weixin_34498545/article/details/112631706
    

    相关文章

      网友评论

        本文标题:matplotlib:根据gtf文件画出基因全部的转录本结构

        本文链接:https://www.haomeiwen.com/subject/smxjvltx.html