Поиск лучшего ответного удара в одном файле BLAST с использованием python

У меня есть выходной файл BLAST outfmt 6 в стандартном формате, я хочу найти способ просмотреть файл, выбрать каждое попадание, найти его ответное попадание и расшифровать, какое попадание лучше всего хранить.

Например:

d = {}
for line in input_file:
    term = line.split('\t')
    qseqid = term[0]
    sseqid = term[1]
    hit = qseqid, sseqid
    recip_hit = sseqid, qseqid
    for line in input_file:
        if recip_hit in line:
            compare both lines
done

Пример ввода (разделители табуляцией):

Seq1    Seq2    80    1000   10    3   1    1000    100    1100    0.0    500
Seq2    Seq1    95    1000   10    3   100    1100    1    1000    1e-100    500

Может ли кто-нибудь дать представление о том, как эффективно решить эту проблему?

Спасибо заранее


python bioinformatics blast
person Gloom    schedule 06.02.2018    source источник
comment
Мои извинения, сейчас добавлю информацию, нет, просто ищу направление   -  person Gloom    schedule 06.02.2018
comment
В этом контексте я предполагаю, что каждая пара Seq1, Seq2 и их обратный аналог не могут встречаться во входном файле более одного раза, верно?   -  person Mr. T    schedule 06.02.2018
comment
Во всяком случае, не в этом наборе данных   -  person Gloom    schedule 06.02.2018
comment
извините - в данном случае это список списков, однако, на мой взгляд, открытый файл был бы гораздо лучшим подходом   -  person Gloom    schedule 06.02.2018
comment
Вы должны обновить свой вопрос с дополнительной информацией. Люди склонны не читать комментарии. И имхо читать прямо из файла не очень хорошая идея. До самой последней строки вы не знаете, какие строки можно отбросить, потому что дубликата нет. Таким образом, вы должны прочитать файл дважды: один раз, чтобы сделать индекс совпадающих пар, и один раз, чтобы прочитать строки для этих совпадающих пар.   -  person Mr. T    schedule 06.02.2018

Ответы (1)


arrow_upward
2
arrow_downward

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

#create a dictionary to store pairs
line_dict = {}
#iterate over your file
for line in open("test.txt", "r"):
    line = line[:-1].split("\t")
    #ignore line, if not at least one value apart from the two sequence IDs
    if len(line) < 3:
        continue
    #identify the two sequences
    seq = tuple(line[0:2])
    #is reverse sequence already in dictionary?
    if seq[::-1] in line_dict:
        #append new line
        line_dict[seq[::-1]].append(line)
    else:
        #create new entry
        line_dict[seq] = [line]

#remove entries, for which no counterpart exists
pairs = {k: v for k, v in line_dict.items() if len(v) > 1}

#and do things with these pairs
for pair, seq in pairs.items():
    print(pair, "found in:")
    for item in seq:
        print(item)

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

Аналогичный подход — хранить все данные в рабочей памяти — использует pandas. Это должно быть быстрее, так как алгоритмы сортировки оптимизированы для панд. Еще одно преимущество pandas заключается в том, что все ваши другие значения уже находятся в столбцах pandas, поэтому дальнейший анализ упрощается. Я определенно предпочитаю версию для панд, но я не знаю, установлена ​​ли она в вашей системе. Чтобы упростить общение, я присвоил a и b столбцам, содержащим последовательности Seq1 и Seq2.

import pandas as pd
#read data into a dataframe
#not necessary: drop the header of the file, use custom columns names
df = pd.read_csv("test.txt", sep='\t', names=list("abcde"), header = 0)

#create a column that joins Seq1 - Seq2 or Seq2 - Seq1 to Seq1Seq2
df["pairs"] = df.apply(lambda row: ''.join(sorted([row["a"], row["b"]])), axis = 1)
#remove rows with no matching pair and sort the database
only_pairs = df[df["pairs"].duplicated(keep = False)].sort_values(by = "pairs")

print(only_pairs)
person Mr. T    schedule 06.02.2018
comment
Это именно то, что я искал! Большое спасибо! - person Gloom; 06.02.2018