extract/parse 大型 multifasta 使用 table (csv, tsv) 进行比对

Question

我经常需要使用从另一个 program/code.

生成的 table 将大型 multifasta 解析为单独的 multifastas 以进行下游对齐

我有一个大的 multifasta (seq.fa):

>sp1_gene1
ATTAC
>sp1_gene2
CCATTA
...
>sp2_gene1
ATTAC
>sp1_gene2
TCGAGT

我有一个 tsv 文件，第一列是基因座名称，随后的列是 headers 的列表。每行中的字段数可能不相等，因为一个物种可能没有它。但是我可以轻松地为每个物种添加 headers 并为缺失数据添加 NA 或类似的东西。 table (genes.tsv):

geneA    sp1_gene3    sp2_gene1
geneB    sp1_gene5    sp2_gene7
...

我想使用基因 table 来创建具有 headers 和序列的单个 multifastas（最好使用第一列中的名称）和序列来获得如下内容：

cat geneA.fa
>sp1_gene3
ATTAC
>sp2_gene1
ATTAC
...
cat geneB.fa
>sp1_gene5
TCGAGT
>sp2_gene7
ATTAC
...

我熟悉 bash（awk、grep、sed），并且仍在学习 R 和 python 用于生物信息学。我最初将 table 拆分为 bash 中的单个文件，将 fasta 转换为 csv，然后 grepping 并加入，但它真的很乱而且并不总是有效。关于可以执行此操作的脚本或包的任何建议？谢谢！

Answer 1

我认为这可以解决您的问题：

sequences = {}

with open("seq.fa") as my_fasta:
    for header in my_fasta:
        seq = next(my_fasta)
        sequences[header[1:].rstrip()] = seq.rstrip()

with open("genes.tsv") as my_tsv:
    for line in my_tsv:
        splitted_line = line.split()
        gene_writer = open("/your/output/Dir/" + splitted_line[0] + ".fa", "w")
        for gene in splitted_line[1:]:
            if gene in sequences:
                gene_writer.write(">" + gene + "\n")
                gene_writer.write(sequences[gene] + "\n")
            else:
                print(gene, "in tsv file but not in fasta")
        gene_writer.close()

分解：

sequences = {}

with open("seq.fa") as my_fasta:
    for header in my_fasta:
        seq = next(my_fasta)
        sequences[header[1:].rstrip()] = seq.rstrip()

这将创建一个字典 sequences，其中键为基因名称，并为序列赋值。像这样：

{'sp1_gene1': 'ATTAC', 'sp1_gene2': 'TCGAGT', 'sp2_gene1': 'ATTAC'}

代码的第二部分遍历 TSV 文件，并为每一行创建一个新的 .fa 文件，并将 fasta 格式的序列添加到该文件。

希望这对您有所帮助。 :)

extract/parse 大型 multifasta 使用 table (csv, tsv) 进行比对

extract/parse large multifasta into alignments using table (csv, tsv)

bioinformatics

biopython

fasta

bioconductor