python - Biopython：在中间序列中添加部分具有对齐的特征

Question

我想在前一个序列的中间添加一段序列（在 gb 文件中），并且所有特征仍然由旧序列索引。

例如：

上一个序列： ATAGCCATTGAATGTGTGTG TGTCCTA GAGGGCCTAAAA

特征：misc_feature 补充（20..27）/gene="Py_ori+A"

我添加TTTTTT位置 10。

新序列：ATAGCCATTG TTTTTT AAGTGTGTG TGTCCTA GAGGGCCTAAAA

特征：misc_feature 补充（26 .. 33）/gene="Py_ori+A"

特征的索引发生了变化，因为特征必须仍然与段 TGTCCTA 有关。我想将新序列保存在新的 gb 文件中。

是否有任何 biopython 函数或方法可以在旧序列中间添加序列段并将添加段的长度添加到添加段之后的特征索引中？

标签： pythonbioinformaticsbiopython

Answer 1

TL;博士

调用+您的切片段（例如a + b）。只要您没有切入某个功能，就应该没问题。

---

长版：

BioPython 支持特征加入。只需调用a + b相应的SeqRecord类即可完成（功能是SeqRecord对象的一部分而不是Seq类。）。

关于具有特征的切片序列，有一个怪癖需要注意。如果您碰巧在特征中进行切片，则该特征将不会出现在结果中SeqRecord。

我试图说明以下代码中的行为。

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation


# THIS IS OK

a = SeqRecord(
    Seq('ACGTA'),
    id='a',
    features=[
        SeqFeature(FeatureLocation(2,4,1), id='f1')
    ]
)

b = SeqRecord(
    Seq('ACGTA'),
    id='b',
    features=[
        SeqFeature(FeatureLocation(2,4,1), id='f2')
    ]
)

c = a + b

print('seq a')
print(a.seq)
print(a.features)

print('\nseq b')
print(b.seq)
print(b.features)

print("\n two distinct features joined in seq c")
print(c.seq)
print(c.features)
print("notice how the second feature has now indices (7,9), instead of 2,4\n")


# BEWARE
# slicing into the feature will remove the feature !

print("\nsliced feature removed")
d = a[:3]
print(d.seq)
print(d.features)
# Seq('ACG')
# []

# However slicing around the feature will preserve it
print("\nslicing out of the feature will preserve it")
e = c[1:6]
print(e.seq)
print(e.features)

输出

seq a
ACGTA
[SeqFeature(FeatureLocation(ExactPosition(2), ExactPosition(4), strand=1), id='f1')]

seq b
ACGTA
[SeqFeature(FeatureLocation(ExactPosition(2), ExactPosition(4), strand=1), id='f2')]

 two distinct features joined in seq c
ACGTAACGTA
[SeqFeature(FeatureLocation(ExactPosition(2), ExactPosition(4), strand=1), id='f1'), SeqFeature(FeatureLocation(ExactPosition(7), ExactPosition(9), strand=1), id='f2')]
notice how the second feature has now indices (7,9), instead of 2,4


sliced feature removed
ACG
[]

slicing out of the feature will preserve it
CGTAA
[SeqFeature(FeatureLocation(ExactPosition(1), ExactPosition(3), strand=1), id='f1')]