使用表（csv，tsv）提取/解析大型multifasta到对齐答案

【问题标题】：extract/parse large multifasta into alignments using table (csv, tsv)使用表（csv，tsv）提取/解析大型multifasta到对齐
【发布时间】：2018-09-23 00:31:42
【问题描述】：

我经常需要使用从其他程序/代码生成的表格将大型 multifasta 解析为单独的 multifasta 以进行下游对齐。

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

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

我有一个 tsv 文件，第一列中包含基因座名称，后续列中包含标题列表。每行中的字段数可能不相等，因为一个物种可能没有它。但是我可以很容易地为每个物种添加标题，并为丢失的数据输入NA 或类似的东西。表（genes.tsv）：

geneA    sp1_gene3    sp2_gene1
geneB    sp1_gene5    sp2_gene7
...

我想使用基因表来创建带有标题和序列的单个 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 的生物信息学。我最初是在 bash 中将表格拆分为单独的文件，将 fasta 转换为 csv，然后 grepping 并加入，但它真的很乱，并不总是有效。关于可以做到这一点的脚本或包的任何建议？谢谢！

【问题讨论】：

抱歉，我没有时间回答这个问题。可以使用biopython 来完成。使用biopython 读入seq.fa 文件，然后循环遍历genes.tsv 文件并使用if 语句检查fasta 标头是否与基因列表匹配，然后附加到列表中。

标签： bioinformatics biopython fasta bioconductor

【解决方案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 格式的序列添加到该文件中。

希望这会有所帮助。 :)

【讨论】：