本文主要工作：
对SBT基因家族的基因结构进行可视化呈现
4.2 基因结构可视化
可视化基因结构，无非就是展示基因区域的组成：3’,5’-UTR，exon或展示CDS，intro n等等，有关它们的详细定义可以参考文章。在这里我们想要可视化展示3’,5’-UTR，exon和基因的关系。需要的原始文件是基因组注释gff3文件，我们想要的信息只有SBT家族基因长度，基因id，exon的长度，在基因上的起始和终止位置，它属于哪个基因等信息。因此，要对gff3文件进行一定的处理，处理思路如下：
import re
gene_dict = {"gene_start":""}
output = ""
output_type = ["exon", "five_prime_UTR", "three_prime_UTR"]
with open("./stat_input.txt", "r") as f:
while True:
line = f.readline()
if not line:
break
else:
line_list = re.split("\t|\n", line)
if line_list[1] == "gene":
gene_dict["gene_start"] = line_list[2]
elif line_list[1] in output_type:
region_start = str(int(line_list[2]) - int(gene_dict["gene_start"]))
region_end = str(int(line_list[3]) - int(gene_dict["gene_start"]))
output += line_list[0] + "\t" + line_list[1] + "\t" + region_start + "\t" + region_end + "\n"
file_output = open("gene_structure_temp.txt", "w")
file_output.write(output)
file_output.close()
除此之外，我们还需要记录基因长度的文件作为绘图的输入文件，这个比较好解决，在之前的脚本中就已经直接生成了，现在我们按照绘制蛋白质模体的思路绘制相同的图:
library(data.table)
library(tidyverse)
# 载入文本
cds_length <- fread("./DATA/aco.sbt.gene.length.txt", header = F)
cds_length <- mutate(cds_length, V3=0)
gene_structure <- fread("./DATA/gene_structure_temp.txt", header = F)
# 画图
ggplot() +
geom_segment(data = gene_structure,
aes(x = as.numeric(V3),
y = V1,
xend = as.numeric(V4),
yend = V1,
color = V2),
linewidth = 2.5,
position = position_nudge(y = 0.2)) +
scale_color_brewer(palette = "Set2") +
geom_segment(data = cds_length,
aes(x = as.numeric(V3),
y = V1,
xend = as.numeric(V2),
yend = V1),
color = "grey",
linewidth = 1) +
scale_x_continuous(expand = c(0,0)) +
labs(y = "Family",
x = "Length",
color = "Structure") +
theme_classic()
| 基因家族分析系列（持续更新）
0.基因家族分析（0）：概念明晰
1.基因家族分析（1）：数据下载与处理
2.基因家族分析（2）：基因家族鉴定与蛋白质性质简单分析
3.基因家族分析（3）：序列比对与进化树构建
4.基因家族分析（4）：基因家族蛋白质模体鉴定与可视化
