包含基因组图的图形

Bio.Graphics 模块依赖于第三方 Python 库 ReportLab。虽然 ReportLab 专注于生成 PDF 文件，但它也可以创建封装的 PostScript (EPS) 和 (SVG) 文件。除了这些基于矢量的图像之外，在安装了一些其他依赖项（如 Python Imaging Library (PIL)）的情况下，ReportLab 也可以输出位图图像（包括 JPEG、PNG、GIF、BMP 和 PICT 格式）。

基因组图

简介

Bio.Graphics.GenomeDiagram 模块已添加到 Biopython 1.50 中，此前它是一个独立的 Python 模块，依赖于 Biopython。基因组图在 Pritchard 等人（2006 年）的《生物信息学》杂志出版物中进行了描述 [Pritchard2006]，其中包含一些示例图像。这里有一个旧手册的 PDF 版本 https://biopython.pythonlang.cn/DIST/docs/GenomeDiagram/userguide.pdf，其中还有一些其他示例。

顾名思义，基因组图旨在绘制整个基因组，尤其是原核生物基因组，以线性图（可选择性地分解为片段以更好地适应）或圆形轮状图的形式。请查看 Toth 等人（2006 年）[Toth2006] 中的图 2，以了解一个很好的示例。它也证明非常适合绘制更小基因组（如噬菌体、质粒或线粒体）的非常详细的图形，例如，请查看 Van der Auwera 等人（2009 年）[Vanderauwera2009] 中的图 1 和 2（显示有额外的手动编辑）。

如果您已将基因组加载为包含大量 SeqFeature 对象的 SeqRecord 对象（例如，从 GenBank 文件中加载，请参阅第 序列注释对象 章和第 序列输入/输出 章），则此模块最易于使用。

图、轨道、特征集和特征

基因组图使用一组嵌套的对象。在顶层，您有一个图对象，表示沿水平轴（或圆形）的序列（或序列区域）。一个图可以包含一个或多个轨道，这些轨道垂直（或在圆形图上呈放射状）堆叠显示。这些轨道通常具有相同的长度，并且表示相同的序列区域。您可以使用一个轨道来显示基因位置，另一个轨道来显示调控区域，第三个轨道来显示 GC 百分比。

最常用的轨道类型将包含特征，这些特征捆绑在特征集中。您可以选择对所有 CDS 特征使用一个特征集，对 tRNA 特征使用另一个特征集。这不是必需的 - 它们都可以放在同一个特征集中，但这使得更新选定特征的属性变得更容易（例如，使所有 tRNA 特征变为红色）。

构建完整图主要有两种方法。首先，自上而下的方法，您创建一个图对象，然后使用其方法添加轨道，并使用轨道方法添加特征集，并使用其方法添加特征。其次，您可以分别创建各个对象（以任何适合您代码的顺序），然后将它们组合在一起。

自上而下的示例

我们将从 GenBank 文件中读取的 SeqRecord 对象绘制整个基因组（请参阅第 序列输入/输出 章）。此示例使用来自鼠疫耶尔森氏菌微鼠生物变种的 pPCP1 质粒，该文件包含在 Biopython 单元测试中，位于 GenBank 文件夹下，或在线从我们的网站获取 NC_005816.gb。

from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO

record = SeqIO.read("NC_005816.gb", "genbank")

我们使用的是自上而下的方法，因此在加载序列后，我们接下来创建一个空的图，然后添加一个（空的）轨道，并在其中添加一个（空的）特征集

gd_diagram = GenomeDiagram.Diagram("Yersinia pestis biovar Microtus plasmid pPCP1")
gd_track_for_features = gd_diagram.new_track(1, name="Annotated Features")
gd_feature_set = gd_track_for_features.new_set()

现在是激动人心的部分 - 我们获取 SeqRecord 中的每个基因 SeqFeature 对象，并使用它在图上生成一个特征。我们将用蓝色对其进行着色，在深蓝色和浅蓝色之间交替。

for feature in record.features:
    if feature.type != "gene":
        # Exclude this feature
        continue
    if len(gd_feature_set) % 2 == 0:
        color = colors.blue
    else:
        color = colors.lightblue
    gd_feature_set.add_feature(feature, color=color, label=True)

现在我们开始实际制作输出文件。这分为两步，首先，我们调用 draw 方法，该方法使用 ReportLab 对象创建所有形状。然后，我们调用 write 方法，该方法将这些形状渲染到请求的文件格式。请注意，您可以以多种文件格式输出

gd_diagram.draw(
    format="linear",
    orientation="landscape",
    pagesize="A4",
    fragments=4,
    start=0,
    end=len(record),
)
gd_diagram.write("plasmid_linear.pdf", "PDF")
gd_diagram.write("plasmid_linear.eps", "EPS")
gd_diagram.write("plasmid_linear.svg", "SVG")

此外，只要您安装了依赖项，您也可以制作位图，例如

gd_diagram.write("plasmid_linear.png", "PNG")

预期的输出如 图 9 所示。

图 9 鼠疫耶尔森氏菌微鼠生物变种质粒 pPCP1 的简单线性图。

请注意，我们设置为 4 的 fragments 参数控制基因组被分成多少个片段。

如果您想绘制圆形图，请尝试以下操作

gd_diagram.draw(
    format="circular",
    circular=True,
    pagesize=(20 * cm, 20 * cm),
    start=0,
    end=len(record),
    circle_core=0.7,
)
gd_diagram.write("plasmid_circular.pdf", "PDF")

预期的输出如 图 10 所示。

图 10 鼠疫耶尔森氏菌微鼠生物变种质粒 pPCP1 的简单圆形图。

这些图不是很有吸引力，但我们才刚刚开始。

自下而上的示例

现在让我们生成完全相同的图形，但使用自下而上的方法。这意味着我们直接创建不同的对象（并且这几乎可以按任何顺序完成），然后将它们组合在一起。

from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO

record = SeqIO.read("NC_005816.gb", "genbank")

# Create the feature set and its feature objects,
gd_feature_set = GenomeDiagram.FeatureSet()
for feature in record.features:
    if feature.type != "gene":
        # Exclude this feature
        continue
    if len(gd_feature_set) % 2 == 0:
        color = colors.blue
    else:
        color = colors.lightblue
    gd_feature_set.add_feature(feature, color=color, label=True)
# (this for loop is the same as in the previous example)

# Create a track, and a diagram
gd_track_for_features = GenomeDiagram.Track(name="Annotated Features")
gd_diagram = GenomeDiagram.Diagram("Yersinia pestis biovar Microtus plasmid pPCP1")

# Now have to glue the bits together...
gd_track_for_features.add_set(gd_feature_set)
gd_diagram.add_track(gd_track_for_features, 1)

您现在可以像以前一样调用 draw 和 write 方法来生成线性图或圆形图，使用上述自上而下示例末尾的代码。这些图形应该完全相同。

没有 SeqFeature 的特征

在上面的示例中，我们使用 SeqRecord 的 SeqFeature 对象来构建我们的图（另请参阅第 特征、位置和位置对象 节）。有时您没有 SeqFeature 对象，而只是想要绘制的特征的坐标。您必须创建最小的 SeqFeature 对象，但这很容易

from Bio.SeqFeature import SeqFeature, SimpleLocation

my_seq_feature = SeqFeature(SimpleLocation(50, 100, strand=+1))

对于链，使用 +1 表示正链，-1 表示反链，None 表示两者。以下是一个简短的独立示例

from Bio.SeqFeature import SeqFeature, SimpleLocation
from Bio.Graphics import GenomeDiagram
from reportlab.lib.units import cm

gdd = GenomeDiagram.Diagram("Test Diagram")
gdt_features = gdd.new_track(1, greytrack=False)
gds_features = gdt_features.new_set()

# Add three features to show the strand options,
feature = SeqFeature(SimpleLocation(25, 125, strand=+1))
gds_features.add_feature(feature, name="Forward", label=True)
feature = SeqFeature(SimpleLocation(150, 250, strand=None))
gds_features.add_feature(feature, name="Strandless", label=True)
feature = SeqFeature(SimpleLocation(275, 375, strand=-1))
gds_features.add_feature(feature, name="Reverse", label=True)

gdd.draw(format="linear", pagesize=(15 * cm, 4 * cm), fragments=1, start=0, end=400)
gdd.write("GD_labels_default.pdf", "pdf")

输出显示在 图 11 的顶部（以默认特征颜色，淡绿色显示）。

请注意，我们在这里使用了 name 参数来指定这些特征的标题文本。这将在下一节中详细讨论。

特征标题

回想一下，我们使用以下代码（其中 feature 是一个 SeqFeature 对象）将特征添加到图中

gd_feature_set.add_feature(feature, color=color, label=True)

在上面的示例中，SeqFeature 注释用于为特征选择合理的标题。默认情况下，将使用 SeqFeature 对象的限定词字典中的以下可能条目：gene、label、name、locus_tag 和 product。更简单地说，您可以直接指定一个名称

gd_feature_set.add_feature(feature, color=color, label=True, name="My Gene")

除了每个特征标签的标题文本外，您还可以选择字体、位置（默认情况下为标记的开始位置，您也可以选择中间或结尾）和方向（仅限线性图，默认情况下旋转 \(45\) 度）。

# Large font, parallel with the track
gd_feature_set.add_feature(
    feature, label=True, color="green", label_size=25, label_angle=0
)

# Very small font, perpendicular to the track (towards it)
gd_feature_set.add_feature(
    feature,
    label=True,
    color="purple",
    label_position="end",
    label_size=4,
    label_angle=90,
)

# Small font, perpendicular to the track (away from it)
gd_feature_set.add_feature(
    feature,
    label=True,
    color="blue",
    label_position="middle",
    label_size=6,
    label_angle=-90,
)

将这三个片段中的每一个与上一节中的完整示例结合起来，应该会得到类似于 图 11 中的轨道的内容。

图 11 显示标签选项的简单基因组图。

淡绿色的顶部图显示了默认标签设置（请参阅第 没有 SeqFeature 的特征 节），而其余部分则显示了标签大小、位置和方向的变化（请参阅第 特征标题 节）。

我们在这里没有显示出来，但是您也可以设置 label_color 来控制标签的颜色（在第 一个不错的示例 节中使用）。

您会注意到默认字体非常小 - 这是有道理的，因为您通常会在页面上绘制许多（小的）特征，而不是像这里所示的几个大的特征。

特征标记

以上示例都只使用了该功能的默认符号，一个普通的方框，这是 GenomeDiagram 上一个公开发布的独立版本中唯一可用的符号。当 GenomeDiagram 添加到 Biopython 1.50 中时，包括了箭头符号。

# Default uses a BOX sigil
gd_feature_set.add_feature(feature)

# You can make this explicit:
gd_feature_set.add_feature(feature, sigil="BOX")

# Or opt for an arrow:
gd_feature_set.add_feature(feature, sigil="ARROW")

Biopython 1.61 添加了三个新的符号，

# Box with corners cut off (making it an octagon)
gd_feature_set.add_feature(feature, sigil="OCTO")

# Box with jagged edges (useful for showing breaks in contains)
gd_feature_set.add_feature(feature, sigil="JAGGY")

# Arrow which spans the axis with strand used only for direction
gd_feature_set.add_feature(feature, sigil="BIGARROW")

这些符号如 图 12 所示。大多数符号都适合在一个边界框（由默认的 BOX 符号给出）内，要么在正向或反向链的轴线上方或下方，要么跨越轴线（高度加倍）用于无链特征。BIGARROW 符号有所不同，它总是跨越轴线，方向取自特征的链。

图 12 展示不同符号的简单 GenomeDiagram。

箭头符号

我们在上一节介绍了箭头符号。有两个额外的选项可以调整箭头的形状，首先是箭头杆的厚度，它以边界框高度的比例给出。

# Full height shafts, giving pointed boxes:
gd_feature_set.add_feature(feature, sigil="ARROW", color="brown", arrowshaft_height=1.0)
# Or, thin shafts:
gd_feature_set.add_feature(feature, sigil="ARROW", color="teal", arrowshaft_height=0.2)
# Or, very thin shafts:
gd_feature_set.add_feature(
    feature, sigil="ARROW", color="darkgreen", arrowshaft_height=0.1
)

结果如 图 13 所示。

图 13 展示箭头杆选项的简单 GenomeDiagram。

其次，箭头头的长度 - 以边界框高度的比例给出（默认为 \(0.5\)，或 \(50\%\)）。

# Short arrow heads:
gd_feature_set.add_feature(feature, sigil="ARROW", color="blue", arrowhead_length=0.25)
# Or, longer arrow heads:
gd_feature_set.add_feature(feature, sigil="ARROW", color="orange", arrowhead_length=1)
# Or, very very long arrow heads (i.e. all head, no shaft, so triangles):
gd_feature_set.add_feature(feature, sigil="ARROW", color="red", arrowhead_length=10000)

结果如 图 14 所示。

图 14 展示箭头头选项的简单 GenomeDiagram。

Biopython 1.61 添加了一个新的 BIGARROW 符号，它总是跨越轴线，对于反向链指向左侧，否则指向右侧。

# A large arrow straddling the axis:
gd_feature_set.add_feature(feature, sigil="BIGARROW")

上面为 ARROW 符号显示的所有杆和箭头头选项也可以用于 BIGARROW 符号。

一个不错的示例

现在让我们回到来自鼠疫耶尔森氏菌微鼠生物变种的 pPCP1 质粒，以及在第自顶向下的示例 节中使用的自顶向下方法，但利用我们现在讨论的符号选项。这次我们将使用箭头表示基因，并将它们与无链特征（作为普通方框）重叠，以显示一些限制酶切位点的位点。

from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO
from Bio.SeqFeature import SeqFeature, SimpleLocation

record = SeqIO.read("NC_005816.gb", "genbank")

gd_diagram = GenomeDiagram.Diagram(record.id)
gd_track_for_features = gd_diagram.new_track(1, name="Annotated Features")
gd_feature_set = gd_track_for_features.new_set()

for feature in record.features:
    if feature.type != "gene":
        # Exclude this feature
        continue
    if len(gd_feature_set) % 2 == 0:
        color = colors.blue
    else:
        color = colors.lightblue
    gd_feature_set.add_feature(
        feature, sigil="ARROW", color=color, label=True, label_size=14, label_angle=0
    )

# I want to include some strandless features, so for an example
# will use EcoRI recognition sites etc.
for site, name, color in [
    ("GAATTC", "EcoRI", colors.green),
    ("CCCGGG", "SmaI", colors.orange),
    ("AAGCTT", "HindIII", colors.red),
    ("GGATCC", "BamHI", colors.purple),
]:
    index = 0
    while True:
        index = record.seq.find(site, start=index)
        if index == -1:
            break
        feature = SeqFeature(SimpleLocation(index, index + len(site)))
        gd_feature_set.add_feature(
            feature,
            color=color,
            name=name,
            label=True,
            label_size=10,
            label_color=color,
        )
        index += len(site)

gd_diagram.draw(format="linear", pagesize="A4", fragments=4, start=0, end=len(record))
gd_diagram.write("plasmid_linear_nice.pdf", "PDF")
gd_diagram.write("plasmid_linear_nice.eps", "EPS")
gd_diagram.write("plasmid_linear_nice.svg", "SVG")

gd_diagram.draw(
    format="circular",
    circular=True,
    pagesize=(20 * cm, 20 * cm),
    start=0,
    end=len(record),
    circle_core=0.5,
)
gd_diagram.write("plasmid_circular_nice.pdf", "PDF")
gd_diagram.write("plasmid_circular_nice.eps", "EPS")
gd_diagram.write("plasmid_circular_nice.svg", "SVG")

预期输出如 显示选定限制酶切位点的 pPCP1 质粒线性图。 和 显示选定限制酶切位点的 pPCP1 质粒环状图。 所示。

图 15 显示选定限制酶切位点的 pPCP1 质粒线性图。

图 16 显示选定限制酶切位点的 pPCP1 质粒环状图。

多条轨道

到目前为止，所有示例都只使用了一条轨道，但你可以有多条轨道 - 例如，在一條軌道上顯示基因，在另一條軌道上顯示重複區域。在本例中，我们将按比例并排显示三个噬菌体基因组，灵感来自 Proux 等人（2002）[Proux2002] 的图 6。我们需要以下三个噬菌体的 GenBank 文件。

• NC_002703 - 乳酸乳球菌噬菌体 Tuc2009，完整基因组（\(38347\) bp）

• AF323668 - 噬菌体 bIL285，完整基因组（\(35538\) bp）

• NC_003212 - 无害李斯特菌 Clip11262，完整基因组，我们只关注整合的原噬菌体 5（长度相似）。

如果你愿意，可以使用 Entrez 下载这些文件，有关更多详细信息，请参阅第EFetch：从 Entrez 下载完整记录 节。对于第三条记录，我们已经确定噬菌体整合到基因组中的位置，并切片记录以提取它（保留特征，请参阅第切片 SeqRecord 节），并且还必须反向互补以匹配前两个噬菌体的方向（同样保留特征，请参阅第反向互补 SeqRecord 对象 节）。

from Bio import SeqIO

A_rec = SeqIO.read("NC_002703.gbk", "gb")
B_rec = SeqIO.read("AF323668.gbk", "gb")
C_rec = SeqIO.read("NC_003212.gbk", "gb")[2587879:2625807].reverse_complement(name=True)

我们模仿的图形使用不同的颜色表示不同的基因功能。一种方法是编辑 GenBank 文件以记录每个特征的颜色偏好 - 这是桑格的 Artemis 编辑器所做的，GenomeDiagram 应该能够理解。然而，这里我们只会硬编码三个颜色列表。

请注意，GenBank 文件中的注释与 Proux 等人显示的注释不完全匹配，他们在一些未注释的基因上进行了绘制。

from reportlab.lib.colors import (
    red,
    grey,
    orange,
    green,
    brown,
    blue,
    lightblue,
    purple,
)

A_colors = (
    [red] * 5
    + [grey] * 7
    + [orange] * 2
    + [grey] * 2
    + [orange]
    + [grey] * 11
    + [green] * 4
    + [grey]
    + [green] * 2
    + [grey, green]
    + [brown] * 5
    + [blue] * 4
    + [lightblue] * 5
    + [grey, lightblue]
    + [purple] * 2
    + [grey]
)
B_colors = (
    [red] * 6
    + [grey] * 8
    + [orange] * 2
    + [grey]
    + [orange]
    + [grey] * 21
    + [green] * 5
    + [grey]
    + [brown] * 4
    + [blue] * 3
    + [lightblue] * 3
    + [grey] * 5
    + [purple] * 2
)
C_colors = (
    [grey] * 30
    + [green] * 5
    + [brown] * 4
    + [blue] * 2
    + [grey, blue]
    + [lightblue] * 2
    + [grey] * 5
)

现在开始绘制它们 - 这次我们在图中添加三条轨道，并且注意到它们被赋予了不同的开始/结束值以反映它们的长度不同（这需要 Biopython 1.59 或更高版本）。

from Bio.Graphics import GenomeDiagram

name = "Proux Fig 6"
gd_diagram = GenomeDiagram.Diagram(name)
max_len = 0
for record, gene_colors in zip([A_rec, B_rec, C_rec], [A_colors, B_colors, C_colors]):
    max_len = max(max_len, len(record))
    gd_track_for_features = gd_diagram.new_track(
        1, name=record.name, greytrack=True, start=0, end=len(record)
    )
    gd_feature_set = gd_track_for_features.new_set()

    i = 0
    for feature in record.features:
        if feature.type != "gene":
            # Exclude this feature
            continue
        gd_feature_set.add_feature(
            feature,
            sigil="ARROW",
            color=gene_colors[i],
            label=True,
            name=str(i + 1),
            label_position="start",
            label_size=6,
            label_angle=0,
        )
        i += 1

gd_diagram.draw(format="linear", pagesize="A4", fragments=1, start=0, end=max_len)
gd_diagram.write(name + ".pdf", "PDF")
gd_diagram.write(name + ".eps", "EPS")
gd_diagram.write(name + ".svg", "SVG")

预期输出如 图 17 所示。

图 17 三个噬菌体的三轨线性图。

这显示了乳酸乳球菌噬菌体 Tuc2009 (NC_002703)、噬菌体 bIL285 (AF323668) 和无害李斯特菌 Clip11262 (NC_003212) 的原噬菌体 5（请参阅第多条轨道 节）。

我很好奇为什么在原始手稿中，最下面的噬菌体中没有标记红色或橙色的基因。另一个重要的一点是，这里噬菌体以不同的长度显示 - 这是因为它们都被绘制到相同的比例（它们确实有不同的长度）。

与已发布图形的关键区别在于，他们在相似蛋白质之间使用颜色编码的链接 - 这也是我们在下一节中将要做的。

轨道之间的交叉链接

Biopython 1.59 添加了在轨道之间绘制交叉链接的功能 - 既有我们将在下面展示的简单线性图，也有分为片段的线性图和环状图。

继续上一节中受 Proux 等人 2002 年 [Proux2002] 图 6 启发的示例，我们需要一个基因对之间的交叉链接列表，以及要使用的分数或颜色。实际上，你可能会从 BLAST 文件中以计算方式提取它，但在这里我手动输入了它们。

我的命名约定继续将三个噬菌体称为 A、B 和 C。以下是我们要在 A 和 B 之间显示的链接，以元组列表的形式给出（百分比相似度得分、A 中的基因、B 中的基因）。

# Tuc2009 (NC_002703) vs bIL285 (AF323668)
A_vs_B = [
    (99, "Tuc2009_01", "int"),
    (33, "Tuc2009_03", "orf4"),
    (94, "Tuc2009_05", "orf6"),
    (100, "Tuc2009_06", "orf7"),
    (97, "Tuc2009_07", "orf8"),
    (98, "Tuc2009_08", "orf9"),
    (98, "Tuc2009_09", "orf10"),
    (100, "Tuc2009_10", "orf12"),
    (100, "Tuc2009_11", "orf13"),
    (94, "Tuc2009_12", "orf14"),
    (87, "Tuc2009_13", "orf15"),
    (94, "Tuc2009_14", "orf16"),
    (94, "Tuc2009_15", "orf17"),
    (88, "Tuc2009_17", "rusA"),
    (91, "Tuc2009_18", "orf20"),
    (93, "Tuc2009_19", "orf22"),
    (71, "Tuc2009_20", "orf23"),
    (51, "Tuc2009_22", "orf27"),
    (97, "Tuc2009_23", "orf28"),
    (88, "Tuc2009_24", "orf29"),
    (26, "Tuc2009_26", "orf38"),
    (19, "Tuc2009_46", "orf52"),
    (77, "Tuc2009_48", "orf54"),
    (91, "Tuc2009_49", "orf55"),
    (95, "Tuc2009_52", "orf60"),
]

同样，对于 B 和 C

# bIL285 (AF323668) vs Listeria innocua prophage 5 (in NC_003212)
B_vs_C = [
    (42, "orf39", "lin2581"),
    (31, "orf40", "lin2580"),
    (49, "orf41", "lin2579"),  # terL
    (54, "orf42", "lin2578"),  # portal
    (55, "orf43", "lin2577"),  # protease
    (33, "orf44", "lin2576"),  # mhp
    (51, "orf46", "lin2575"),
    (33, "orf47", "lin2574"),
    (40, "orf48", "lin2573"),
    (25, "orf49", "lin2572"),
    (50, "orf50", "lin2571"),
    (48, "orf51", "lin2570"),
    (24, "orf52", "lin2568"),
    (30, "orf53", "lin2567"),
    (28, "orf54", "lin2566"),
]

对于第一个和最后一个噬菌体，这些标识符是基因座标签，对于中间的噬菌体，没有基因座标签，所以我使用了基因名称。以下小助手函数使我们能够使用基因座标签或基因名称查找特征。

def get_feature(features, id, tags=["locus_tag", "gene"]):
    """Search list of SeqFeature objects for an identifier under the given tags."""
    for f in features:
        for key in tags:
            # tag may not be present in this feature
            for x in f.qualifiers.get(key, []):
                if x == id:
                    return f
    raise KeyError(id)

我们现在可以将这些标识符对列表转换为 SeqFeature 对，从而找到它们的位置坐标。我们现在可以将所有这些代码以及以下代码片段添加到前面的示例中（就在 gd_diagram.draw(...) 行之前 - 请参阅完成的示例脚本 Proux_et_al_2002_Figure_6.py，该脚本包含在 Biopython 源代码的 Doc/examples 文件夹中）以将交叉链接添加到图形中。

from Bio.Graphics.GenomeDiagram import CrossLink
from reportlab.lib import colors

# Note it might have been clearer to assign the track numbers explicitly...
for rec_X, tn_X, rec_Y, tn_Y, X_vs_Y in [
    (A_rec, 3, B_rec, 2, A_vs_B),
    (B_rec, 2, C_rec, 1, B_vs_C),
]:
    track_X = gd_diagram.tracks[tn_X]
    track_Y = gd_diagram.tracks[tn_Y]
    for score, id_X, id_Y in X_vs_Y:
        feature_X = get_feature(rec_X.features, id_X)
        feature_Y = get_feature(rec_Y.features, id_Y)
        color = colors.linearlyInterpolatedColor(
            colors.white, colors.firebrick, 0, 100, score
        )
        link_xy = CrossLink(
            (track_X, feature_X.location.start, feature_X.location.end),
            (track_Y, feature_Y.location.start, feature_Y.location.end),
            color,
            colors.lightgrey,
        )
        gd_diagram.cross_track_links.append(link_xy)

这段代码中有几个重要的部分。首先，GenomeDiagram 对象有一个 cross_track_links 属性，它只是一个 CrossLink 对象的列表。每个 CrossLink 对象采用两组轨道特定的坐标（这里以元组形式给出，你可以选择使用 GenomeDiagram.Feature 对象代替）。你可以选择提供一个颜色、边框颜色，以及说这个链接是否应该反向绘制（对于显示反转很有用）。

你还可以看到我们如何将 BLAST 百分比同一性得分转换为颜色，在白色（\(0\%\)）和深红色（\(100\%\)）之间进行插值。在本例中，我们没有遇到任何交叉链接重叠的问题。解决此问题的一种方法是在 ReportLab 中使用透明度，方法是使用其 alpha 通道设置的颜色。但是，这种类型的阴影颜色方案与重叠透明度结合起来会很难理解。预期输出如 图 18 所示。

图 18 三轨线性图加上基本的交叉链接。

三条轨道显示了乳酸乳球菌噬菌体 Tuc2009 (NC_002703)、噬菌体 bIL285 (AF323668) 和无害李斯特菌 Clip11262 (NC_003212) 的原噬菌体 5，以及按百分比同一性着色的基本交叉链接（请参阅第轨道之间的交叉链接 节）。

在 Biopython 中还可以做很多事情来帮助改进这张图。首先，在这种情况下，交叉链接是在以链特定方式绘制的蛋白质之间。在特征轨道上添加一个背景区域（使用“BOX”符号的特征）以扩展交叉链接会有所帮助。此外，我们可以减少特征轨道的垂直高度，以便为链接分配更多空间 - 一种方法是为空轨道分配空间。此外，在这种情况下，由于没有大的基因重叠，我们可以使用跨轴的 BIGARROW 符号，这使我们能够进一步减少轨道所需的垂直空间。这些改进在示例脚本 Proux_et_al_2002_Figure_6.py 中进行了演示，该脚本包含在 Biopython 源代码的 Doc/examples 文件夹中。预期输出如 图 19 所示。

图 19 三轨线性图加上阴影交叉链接。

三条轨道显示了乳酸乳球菌噬菌体 Tuc2009 (NC_002703)、噬菌体 bIL285 (AF323668) 和无害李斯特菌 Clip11262 (NC_003212) 的原噬菌体 5，以及按百分比同一性着色的交叉链接（请参阅第轨道之间的交叉链接 节）。

除此之外，你可能想要在矢量图像编辑器中手动完成的收尾工作包括微调基因标签的放置，以及添加其他自定义注释，例如突出显示特定区域。

虽然在本例中并不十分必要，因为没有交叉链接重叠，但在 ReportLab 中使用透明颜色是一种非常有用的技术，用于叠加多个链接。但是，在这种情况下应避免使用阴影颜色方案。

更多选项

您可以控制刻度线以显示比例 - 毕竟每个图表都应该显示其单位，以及灰色轨迹标签的数量。

此外，我们到目前为止只使用了 FeatureSet。GenomeDiagram 还具有 GraphSet，可用于显示折线图、条形图和热图（例如，在与特征平行的轨道上显示 GC% 图）。

这些选项目前尚未涵盖，因此我们现在将您转至独立版 GenomeDiagram 附带的 用户指南 (PDF)（但请先阅读下一节）以及文档字符串。

转换旧代码

如果您使用独立版 GenomeDiagram 编写了旧代码，并且想将其切换到 Biopython 中包含的新版本，那么您需要进行一些更改 - 最重要的是更改导入语句。

此外，旧版本的 GenomeDiagram 仅使用颜色和中心的英国拼写（colour 和 centre）。您需要更改为美式拼写，尽管多年来 Biopython 版本的 GenomeDiagram 支持这两种拼写。

例如，如果您以前有

from GenomeDiagram import GDFeatureSet, GDDiagram

gdd = GDDiagram("An example")
...

您可以像这样切换导入语句

from Bio.Graphics.GenomeDiagram import FeatureSet as GDFeatureSet, Diagram as GDDiagram

gdd = GDDiagram("An example")
...

希望这足够了。从长远来看，您可能希望切换到新名称，但您需要更改更多代码

from Bio.Graphics.GenomeDiagram import FeatureSet, Diagram

gdd = Diagram("An example")
...

或者

from Bio.Graphics import GenomeDiagram

gdd = GenomeDiagram.Diagram("An example")
...

如果您遇到困难，请在 Biopython 邮件列表中寻求建议。需要注意的是，我们尚未包含旧模块 GenomeDiagram.GDUtilities。这包括了许多与 GC% 相关的函数，这些函数可能会在稍后合并到 Bio.SeqUtils 中。

染色体

Bio.Graphics.BasicChromosome 模块允许绘制染色体。Jupe 等人（2012 年）[Jupe2012]（开放获取）中有一个示例，使用颜色突出显示不同的基因家族。

简单染色体

这是一个非常简单的示例 - 我们将使用拟南芥。

图 20 拟南芥的简单染色体图。

您可以跳过这一部分，但首先我从 NCBI 的 FTP 站点 ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Arabidopsis_thaliana/ 下载了五个测序的染色体作为五个独立的 FASTA 文件，然后使用 Bio.SeqIO 解析它们以找出它们的长度。您可以使用 GenBank 文件来完成此操作（接下来的示例使用它们来绘制特征），但是如果您只需要长度，使用整个染色体的 FASTA 文件更快

from Bio import SeqIO

entries = [
    ("Chr I", "CHR_I/NC_003070.fna"),
    ("Chr II", "CHR_II/NC_003071.fna"),
    ("Chr III", "CHR_III/NC_003074.fna"),
    ("Chr IV", "CHR_IV/NC_003075.fna"),
    ("Chr V", "CHR_V/NC_003076.fna"),
]
for name, filename in entries:
    record = SeqIO.read(filename, "fasta")
    print(name, len(record))

这给出了五个染色体的长度，我们现在将在以下简短演示 BasicChromosome 模块中使用它们

from reportlab.lib.units import cm
from Bio.Graphics import BasicChromosome

entries = [
    ("Chr I", 30432563),
    ("Chr II", 19705359),
    ("Chr III", 23470805),
    ("Chr IV", 18585042),
    ("Chr V", 26992728),
]

max_len = 30432563  # Could compute this from the entries dict
telomere_length = 1000000  # For illustration

chr_diagram = BasicChromosome.Organism()
chr_diagram.page_size = (29.7 * cm, 21 * cm)  # A4 landscape

for name, length in entries:
    cur_chromosome = BasicChromosome.Chromosome(name)
    # Set the scale to the MAXIMUM length plus the two telomeres in bp,
    # want the same scale used on all five chromosomes so they can be
    # compared to each other
    cur_chromosome.scale_num = max_len + 2 * telomere_length

    # Add an opening telomere
    start = BasicChromosome.TelomereSegment()
    start.scale = telomere_length
    cur_chromosome.add(start)

    # Add a body - using bp as the scale length here.
    body = BasicChromosome.ChromosomeSegment()
    body.scale = length
    cur_chromosome.add(body)

    # Add a closing telomere
    end = BasicChromosome.TelomereSegment(inverted=True)
    end.scale = telomere_length
    cur_chromosome.add(end)

    # This chromosome is done
    chr_diagram.add(cur_chromosome)

chr_diagram.draw("simple_chrom.pdf", "Arabidopsis thaliana")

这应该创建一个非常简单的 PDF 文件，如 图 20 所示。此示例故意简短。下一个示例显示了感兴趣特征的位置。

带注释的染色体

图 21 显示 tRNA 基因的拟南芥染色体图。

从前面的示例继续，我们还将显示 tRNA 基因。我们将通过解析五个拟南芥染色体的 GenBank 文件来获取它们的位置。您需要从 NCBI 的 FTP 站点 ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Arabidopsis_thaliana/ 下载这些文件，并保留子目录名称或编辑下面的路径

from reportlab.lib.units import cm
from Bio import SeqIO
from Bio.Graphics import BasicChromosome

entries = [
    ("Chr I", "CHR_I/NC_003070.gbk"),
    ("Chr II", "CHR_II/NC_003071.gbk"),
    ("Chr III", "CHR_III/NC_003074.gbk"),
    ("Chr IV", "CHR_IV/NC_003075.gbk"),
    ("Chr V", "CHR_V/NC_003076.gbk"),
]

max_len = 30432563  # Could compute this from the entries dict
telomere_length = 1000000  # For illustration

chr_diagram = BasicChromosome.Organism()
chr_diagram.page_size = (29.7 * cm, 21 * cm)  # A4 landscape

for index, (name, filename) in enumerate(entries):
    record = SeqIO.read(filename, "genbank")
    length = len(record)
    features = [f for f in record.features if f.type == "tRNA"]
    # Record an Artemis style integer color in the feature's qualifiers,
    # 1 = Black, 2 = Red, 3 = Green, 4 = blue, 5 =cyan, 6 = purple
    for f in features:
        f.qualifiers["color"] = [index + 2]

    cur_chromosome = BasicChromosome.Chromosome(name)
    # Set the scale to the MAXIMUM length plus the two telomeres in bp,
    # want the same scale used on all five chromosomes so they can be
    # compared to each other
    cur_chromosome.scale_num = max_len + 2 * telomere_length

    # Add an opening telomere
    start = BasicChromosome.TelomereSegment()
    start.scale = telomere_length
    cur_chromosome.add(start)

    # Add a body - again using bp as the scale length here.
    body = BasicChromosome.AnnotatedChromosomeSegment(length, features)
    body.scale = length
    cur_chromosome.add(body)

    # Add a closing telomere
    end = BasicChromosome.TelomereSegment(inverted=True)
    end.scale = telomere_length
    cur_chromosome.add(end)

    # This chromosome is done
    chr_diagram.add(cur_chromosome)

chr_diagram.draw("tRNA_chrom.pdf", "Arabidopsis thaliana")

它可能会警告您标签过于靠近 - 查看 Chr I 的正向链（右侧），但它应该创建一个彩色的 PDF 文件，如 图 21 所示。