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