如何读取基因对列表并为每一行编写一个fasta文件答案

【问题标题】：How to read a list of gene pairs and write a fasta file for each line如何读取基因对列表并为每一行编写一个fasta文件
【发布时间】：2023-02-09 19:23:19
【问题描述】：

我是生物信息学的新手，非常希望能得到一些帮助！

我有一个大的 multi-fasta 文件（genes.faa），像这样：

>gene1_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene2_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene3_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene4_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
(...)

以及基因对列表 (gene.pairs.txt)，每行两个基因由制表符分隔：

gene13_A \t gene33_B
gene2_A \t gene48_B
gene56_A \t gene2_B

我需要一种方法来读取基因对列表并为基因对列表的每一行创建一个 fasta 文件。所以，在这种情况下，我会有 3 个 fasta 文件（输出 fasta 文件的名称并不重要），如下所示：

法斯塔1

>gene13_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene33_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC

法斯塔2

>gene2_A 
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene48_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC

法斯塔3

>gene56_A 
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene2_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC

我试图用 python 编写脚本，但我找不到一种方法来循环读取列表并为每一行编写一个 fasta 文件。非常感谢您的帮助！

【问题讨论】：

请编辑问题以向我们展示您最近尝试的代码以及您遇到困难的地方。另请参阅：How to Ask 和 help center。此外，您可能想使用Biopython，特别是Bio.SeqIO。 Biopython 可以很容易地安装，例如使用conda。

标签： list bioinformatics biopython fasta

【解决方案1】：

我的尝试，肯定有更快更好的方法来完成它，这让我想知道是否有一种方法可以跳过大字典的创建：sequences = { i.id : i for i in SeqIO.parse('big_fasta_2.fa', 'fasta')}。

我正在使用 Biopython 库来解析 fasta 文件并将它们写入https://biopython.org/、https://github.com/biopython/biopython；反正：

输入'big_fasta_2.fa'：

>gene1_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene2_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene2_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene3_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene4_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAC
>gene13_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAA
>gene33_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAY
>gene48_B
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAW
>gene56_A
MCTGTRNKIIRTCDNCRKRKIKCDRKRPAP

输入"gene_pairs_3.txt"：

gene13_A    gene33_B
gene1344_A  gene33_B
gene2_A gene48_B
gene23333_A gene48_B
gene56_A    gene2_B

代码：

from Bio import SeqIO,  __version__

print('Biopython version : ', __version__)


sequences = { i.id : i for i in SeqIO.parse('big_fasta_2.fa', 'fasta')}

print(sequences)


file = open("gene_pairs_3.txt","r")


cnt = 1
for line in file:
    
    
    a, b  =  line.split()
    print('++++++++++++')
    print('pairs N° : ', cnt)
    print(a)
    print(b)
    
    if a in sequences:
        print('ok A')
        print(sequences[a])
        
        if b in sequences:
            print('ok B')
            print(sequences[b])
            
            SeqIO.write([sequences[a],sequences[b]] , 'Fasta'+str(cnt)+'.fa' , 'fasta')
            
            print('
written file : ' ,'Fasta'+str(cnt)+'.fa' )
            
            cnt += 1
        else:
            
            print('No B')
    
    else:
        print('No A')
        continue
    
    print('-----------
')

查看输出文件，看看它们是否符合您的预期。

【讨论】：