python - 使用 python 在单个 BLAST 文件中找到最佳相互命中

Question

我有一个标准格式的 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

谁能提供任何有关如何有效解决此问题的见解？

提前致谢

Answer 1

您可以解决您的问题以找到这些对并比较如下行:

#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 的另一个优点是所有其他值都已经在 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)