如何使用Python将两个FASTA文件相交

Question

请给出您的提示-我在这段代码上花了很多时间，但一直无法解决它

我有一个以我想要的格式输出的代码，但是它显示的结果不是100%正确

我有两个大的FASTA文件，如下所示（#之后的信息不在实际文件中-为了更好地理解，我将它们放在这个问题中）：

f1.fa

>seq1
ATCGCGT  #--> 6mers are ATCGCG and TCGCGT // reverse complements are CGCGAT and ACGCGA
>seq20
CCATGAA  #--> 6mers are CCATGA and CATGAA // reverse complements are TCATGG and TTCATG
>seq3
AACATGA  #--> 6mers are AACATG and ACATGA // reverse complements are CATGTT and TCATGT

f2.fa

>seq9
NNCGCGATNNCGCGATNN
>seq85
NNTTCATG

我想做的是：

1. 从file1读取reach序列的6个字符（6个字符窗口）
2. 对6mers进行反向补充（这不是问题）
3. 在file2中查找任何序列（这些是目标序列），查看是否有任何反向补码6mer在其中，然后用1和0标记这些序列
4. 以表格形式报告结果

因此，根据上面的示例序列-来自seq1（CGCGAT）的一个反向补码6mer位于文件2的seq9内-来自seq20（TTCATG）的一个反向补码6mer位于文件2的seq85内

因此，所需的输出应如下所示：

         seq1  seq20  seq3
name                    
seq9      1      0     0
seq85     0      1     0

但目前，我得到：

         seq1  seq20  seq3
name                    
seq9      0      0     0
seq85     0      1     0

我的代码是：

from Bio import SeqIO
import pandas as pd


def same_seq(a_record, brecord):
    window_size = 7
    for j in range(0, len(a_record.seq) - window_size + 1):  
        for k in (a_record.seq.reverse_complement()[j: j + 6].split()):  
            return brecord.seq.find(k) != -1


if __name__ == '__main__':
    records = list(SeqIO.parse("f1.fa", "fasta"))
    target_records = list(SeqIO.parse("f2.fa", "fasta"))
    rows_list = []
    for target_record in target_records:
        new_row = {'name': target_record.name}
        for record in records:
            if same_seq(record, target_record):
                new_row[record.name] = 1
            else:
                new_row[record.name] = 0
        rows_list.append(new_row)
    df = pd.DataFrame(rows_list)
    df = df.set_index(["name"])
    print(df)

我不确定是代码的第一部分还是第二部分导致了问题。如果你能帮我修改代码就太好了。多谢各位

Answer 1

看起来你的代码中有一个bug。按如下方式更改same_seq函数：

def same_seq(a_record, brecord):
    window_size=6
    for j in range(len(a_record.seq)- window_size):
        for k in (a_record.seq.reverse_complement()[j: j + 6].split()):
            if brecord.seq.find(k) != -1:
                    return 1
    return 0

希望它能起作用