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

【问题标题】：How to intersect two FASTA files using Python如何使用 Python 相交两个 FASTA 文件
【发布时间】：2021-07-05 11:17:30
【问题描述】：

请给出您的提示 - 我在这段代码上花了很多时间，但未能解决它。

我有一个代码以我想要的格式提供输出，但是它显示的结果并非 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（6个字符窗口）
进行 6mer 的反向补码（这不是问题）
在file2 中查找任何序列（这些是目标序列）以查看其中是否有任何反向互补 6mer，然后将这些序列相互标记为 1，否则标记为 0
以表格形式报告结果。

所以基于上面的示例序列 - 来自 seq1 (CGCGAT) 的反向补码 6mer 之一在 file2 的 seq9 内 - 来自 seq20 (TTCATG) 的反向补码 6mer 之一也在 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)

我不确定是代码的第一部分还是第二部分导致了问题。如果您能帮助我修复我的代码，那就太好了。谢谢

【问题讨论】：

如果您正在开发 Python 解决方案，为什么要使用 bash 和 awk 标签？
我也可以使用 bash 和 awk - 如果有快捷方式。给您带来的不便敬请谅解

标签： python bash awk output intersection

【解决方案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

希望它有效！

【讨论】：