Изменить местоположение функции генбанка

Question

Редактировать: я знаю feature.type даст ген /CDS и feature.qualifiers затем даст "db_xref"/"locus_tag"/"вывод" и т. д. Есть ли feature. объект, который позволит мне получить доступ к местоположению (например: [5240:7267](+)) прямо?

Этот URL дает немного больше информации, хотя я не могу понять, как использовать его для своих целей... http://biopython.org/DIST/docs/api/Bio.SeqFeature.SeqFeature-class.html

Исходное сообщение:

Я пытаюсь изменить расположение объектов в файле GenBank. По сути, я хочу изменить следующий бит файла GenBank:

 gene            5240..7267
                 /db_xref="GeneID:887081"
                 /locus_tag="Rv0005"
                 /gene="gyrB"
 CDS             5240..7267
                 /locus_tag="Rv0005"
                 /inference="protein motif:PROSITE:PS00177"
                 ...........................

в

 gene            5357..7267
                 /db_xref="GeneID:887081"
                 /locus_tag="Rv0005"
                 /gene="gyrB"
 CDS             5357..7267
                 /locus_tag="Rv0005"
                 /inference="protein motif:PROSITE:PS00177"
                 .............................

Обратите внимание на изменения с 5240 до 5357

Пока что, отыскивая Интернет и Stackru, я имею:

from Bio import SeqIO
gb_file = "mtbtomod.gb"
gb_record = SeqIO.parse(open(gb_file, "r+"), "genbank")
rvnumber = 'Rv0005'
newstart = 5357

final_features = []

for record in gb_record:
  for feature in record.features:
    if feature.type == "gene":
        if feature.qualifiers["locus_tag"][0] == rvnumber:
            if feature.location.strand == 1:
                feature.qualifiers["amend_position"] = "%s:%s" % (newstart, feature.location.end+1)
            else:
                # do the reverse for the complementary strand
    final_features.append(feature)
  record.features = final_features
  with open("testest.gb","w") as testest:
    SeqIO.write(record, testest, "genbank")

Это в основном создает новый классификатор с именем "mend_position"... однако, я бы хотел изменить местоположение напрямую (с созданием или без создания нового файла...)

Rv0005 - это просто пример locus_tag, который мне нужно обновить. У меня есть еще около 600 мест для обновления, что объясняет необходимость скрипта. Помогите!

Answer 1

Хорошо, у меня есть кое-что, что теперь полностью работает. Я выложу код на случай, если кому-нибудь понадобится что-то подобное

__author__ = 'Kavin'

from Bio import SeqIO
from Bio import SeqFeature
import xlrd
import sys
import re

workbook = xlrd.open_workbook(sys.argv[2])
sheet = workbook.sheet_by_index(0)
data = [[sheet.cell_value(r, c) for c in range(sheet.ncols)] for r in range(sheet.nrows)]

# Create dicts to store TSS data
TSS = {}
row = {}
# For each entry (row), store the startcodon and strand information
for i in range(2, sheet.nrows - 1):
    if data[i][5] < -0.7:   # Ensures BASS score is within significant range
        Gene = data[i][0]
        row['Direction'] = str(data[i][3])
        row['StartCodon'] = int(data[i][4])
        TSS[str(Gene)] = row
        row = {}
    else:
        i += 1

# Create an output filename based on input filename
outfile_init = re.search('(.*)\.(\w*)', sys.argv[1])
outfile = str(outfile_init.group(1)) + '_modified.' + str(outfile_init.group(2))

final_features = []
for record in SeqIO.parse(open(sys.argv[1], "r"), "genbank"):
    for feature in record.features:
        if feature.type == "gene" or feature.type == "CDS":
            if TSS.has_key(feature.qualifiers["locus_tag"][0]):
                newstart = TSS[feature.qualifiers["locus_tag"][0]]['StartCodon']
                if feature.location.strand == 1:
                    feature.location = SeqFeature.FeatureLocation(SeqFeature.ExactPosition(newstart - 1),
                                                                  SeqFeature.ExactPosition(
                                                                      feature.location.end.position),
                                                                  feature.location.strand)
                else:
                    feature.location = SeqFeature.FeatureLocation(
                        SeqFeature.ExactPosition(feature.location.start.position),
                        SeqFeature.ExactPosition(newstart), feature.location.strand)
        final_features.append(feature)  # Append final features
    record.features = final_features
    with open(outfile, "w") as new_gb:
        SeqIO.write(record, new_gb, "genbank")

Это предполагает использование, такое как python program.py <genbankfile> <excel spreadsheet>

Это также предполагает электронную таблицу следующего формата:

Синоним гена Tuberculist_annotated_start Ориентация Re-annotated_start BASS_score

Rv0005 gyrB 5240 + 5357 -1,782

Rv0012 Rv0012 14089 + 14134 -1,553

Rv0018c pstP 23181 - 23172 -2,077

Rv0032 bioF2 34295 + 34307 -0,842

Rv0037c Rv0037c 41202 - 41163 -0,554

Answer 2

Я думаю, что это может быть сделано с нативным синтаксом биопиона, регулярное выражение не требуется, минимальный рабочий пример здесь:

from Bio import SeqIO
from Bio import SeqFeature
import copy

gbk = SeqIO.read('./test_gbk', 'gb')

index = 1

feature_to_change = copy.deepcopy(gbk.features[index]) #depends what feature you want to change, 
#can also be done with loop if you want to change them all, or by some function...

new_start = 0
new_end = 100

new_feature_location = SeqFeature.FeatureLocation(new_start, new_end, feature.location.strand) #create a new feature location object

feature_to_change.location = new_feature_location #change old feature location

del gbk.features[index] #remove changed feature

gkb.features.append(feature_to_change) #add identical feature with new location

gbk.features = sorted(gbk.features, key = lambda feature: feature.location.start) # if you want them sorted by the start of the location, which is the usual case

SeqIO.write(gbk, './test_gbk_with_new_feature', 'gb')

Answer 3

Итак, вы можете попробовать что-то вроде ниже. Так как количество изменений будет равно количеству CDS/ генов, найденных в файле. Вы можете прочитать местоположение / позиции из CSV / текстового файла и сделать список, как я сделал вручную change_values,

import re
f = open("test.txt")
change_values=["1111", "2222"]
flag = True
next = 0
for i in f.readlines():
    if i.startswith(' CDS') or i.startswith(' gene'):
        out = re.sub(r"\d+", str(change_values[next]), i)
                #Instead of print write
        print out
        flag = not flag
        if flag==True:
            next += 1
    else:
                #Instead of print write
        print i

Пример файла Amy test.txt выглядит следующим образом:

 gene            5240..7267
                 /db_xref="GeneID:887081"
                 /locus_tag="Rv0005"
                 /gene="gyrB"
 CDS             5240..7267
                 /locus_tag="Rv0005"
                 /inference="protein motif:PROSITE:PS00177"
                 ...........................

 gene            5240..7267
                 /db_xref="GeneID:887081"
                 /locus_tag="Rv0005"
                 /gene="gyrB"
 CDS             5240..7267
                 /locus_tag="Rv0005"
                 /inference="protein motif:PROSITE:PS00177"
                 ...........................

Надеюсь, это решит вашу проблему. Ура!