本文主要工作：
对SBT基因家族的顺式作用元件进行可视化

4.3 顺式作用元件可视化

顺式作用元件是位于基因组上的非编码序列，它一般位于基因序列的附近，主要作用是调控基因的表达。在文章中，作者的构思巧妙，将顺势作用元件依据其对应的植物生长发育功能分为三类，进行可视化展示。仿照文章思路，我们首先选取SBT基因家族各基因上游2000bp的序列。这里我写了一个脚本来完成这件事。需要用到的文件是凤梨基因组文件和SBT基因家族的gtf文件。

#!/usr/bin/env python

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
import re

gene_family_dict = {}

# 1.gtf文件处理
with open("./aco.sbt.gtf", "r") as f:
    while True:
        file_line = f.readline()
        
        if not file_line:
            break

        else:
            file_line_list = re.split("\n|\t|\"", file_line)
            
            if file_line_list[2] == "transcript":
                gene_name = file_line_list[-3]
                chromosome = file_line_list[0]
                start = file_line_list[3]
                gene_family_dict[gene_name] = {"chr":chromosome, "start":start, "name":gene_name}

# 2.fasta文件读入
sequences_dict = SeqIO.index("../11.data/aco.fa", "fasta")

# 3.输出上游序列
for gene,info in gene_family_dict.items():
    start_point = int(info["start"]) - 2001
    end_point = int(info["start"]) - 1
    upstream = sequences_dict[info["chr"]].seq[start_point:end_point]
    upstream = SeqRecord(upstream, id=info["name"], description = "")

    with open("aco.sbt.fa", "a") as output_fa:
        SeqIO.write(upstream, output_fa, "fasta")

当我们得到提取的序列后，在PlantCare中上传该序列，并进行鉴定。

在这里只需要填写你的电子邮件地址以及你需要鉴定的序列，我们在这里上传上文提到的提取的上游2000bp的序列。随后在电子邮件中我们会得到一个表格，这个表格非常杂乱无章，在记事本中，打开来是这样子：

library(data.table)
library(tidyverse)
library(patchwork)

# 1.读入文件
cis <- fread("./DATA/plantCARE_output_PlantCARE_20618.txt",
             sep = "\t", header = F)

# 2.文件处理
cis_element <- c("ABRE", "AuxRR-core", "ERE", "CGTCA-motif", "GARE-motif", 
                 "TCA-element", "TGA-element", "ARE", "LTR", "MBS", "MYB",
                 "TC-rich", "W-box", "WUN-motif", "AE-box", "ACE", "ATCT-motif",
                 "I-box", "Box-4", "G-box", "GATA-motif")

new_cis <- filter(cis, V2 %in% cis_element) %>%
  select(name = V1, type = V2, length = V5) 
new_cis <- aggregate(length~type + name, new_cis, sum)
unique(new_cis$name)

# 画图
p1 <- ggplot() +
  geom_col(data = filter(new_cis, name %in% unique(new_cis$name)[1:17]),
           aes(x = name, fill = type, y = length),
           position = "stack",
           color = "black") +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, 
                                   vjust = 0.5)) 

p2 <- ggplot() +
  geom_col(data = filter(new_cis, name %in% unique(new_cis$name)[18:36]),
           aes(x = name, fill = type, y = length),
           position = "stack",
           color = "black") +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, 
                                   vjust = 0.5)) 

p3 <- ggplot() +
  geom_col(data = filter(new_cis, name %in% unique(new_cis$name)[37:54]),
           aes(x = name, fill = type, y = length),
           position = "stack",
           color = "black") +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, 
                                   vjust = 0.5)) 

p4 <- p1 / p2 / p3
p4 <- p4 + plot_layout(guides = 'collect')

最后得到我们想要的图像，虽然颜色很丑，但是这个是可以调整的，总体上思路是正确的。

卖家秀（论文原图）：