如何使用 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

我想做的是：

从file1
进行6mers的反向互补（这不是问题）
在file2中寻找任何序列（这些是目标序列），看看是否有任何反向互补6mers在里面，然后用1标记这些序列，其他用0
将结果报告为 table。

所以基于上面的示例序列 - 来自 seq1 (CGCGAT) 的反向互补 6mers 之一在 file2 的 seq9 内部 - 也是来自 seq20 (TTCATG) 的反向互补 6mers 之一在 file2.

的 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

您的代码中似乎存在错误。改变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

希望有用！

如何使用 Python 交叉两个 FASTA 文件

How to intersect two FASTA files using Python

python

bash

awk

intersection

output