ggplot集锦

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

2021-08-02  本文已影响0人  所以suoyi
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
上一篇下一篇

猜你喜欢

热点阅读