This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Two multifasta files with the same sequences, but different headers. How to change headers to match?

I have a large multifasta file (about 125,000 sequences) and a smaller multifasta file (about 100 sequences). All sequences in the smaller multifasta file are found in the larger file, but the headers are different. I have many (thousands) of such smaller multifasta files. How can I search the larger file for the sequences found in the smaller and then exchange the header? I would ideally be able to print out a smaller multifasta file that would be identical to the one I started with, just with the headers found in the larger file. All sequences in both files have been linearized- that is, they are a single line. Thanks!

perl multifasta fasta

3 answers

The following might work for you.

grep -f other.fa reference.fa -B1 --no-group-separator

or

grep -f other.fa reference.fa -B1 | grep -v -- "^--$"

if --no-group-separator is not available in your version of grep.

Note that this will match substrings, which could be unwanted or undesired depending on your use case.

Thanks! Is there any way to change the grep command to only print out the headers themselves, instead of also including the associated sequences? I am not very familiar with grep. Thanks again.

You could use grep "^>" , which would get you just the headers.

Looks like the original grep solution worked? I am curious since you had large files.

I couldn't do it locally because grep runs out of memory, but I can run it successfully over the server that I have access to. I need to learn more about grep, awk, sed, etc. They seem quick, powerful, and really simple. Maybe I am deceiving myself about the simple part though! Do you have any resources to suggest for learning more about these type of commands? Thanks again.

Solution of grep -B1 by @roblogan6 is very cool, however as he/she said, it matched substrings instead of whole sequences. Besides, sequences both in the big and small files must be in single-line format.

Here's is a robust preciser solution with SeqKit:

Big file:

$ cat big.f
>seq1
ACTACGACGTC
TAGCGTA
>seq2
CGACGATCTAC
GTAGCTAGAT
>seq3
ACGTCTGACGT
>seq4 containing seq3
ACGTACGTCTG
ACGTCC

Small file:

$ cat small.fa
>seq_abc
ACTACGACGTC
TAGCGTA
>seq_123
ACGTCTGACGT

Precisely matching by sequences:

$ seqkit grep -s -i -f <(seqkit seq -s -w 0 small.fa) -w 70 big.fa    
>seq1
ACTACGACGTCTAGCGTA
>seq3
ACGTCTGACGT

Here's the "long-option" version:

seqkit grep --by-seq --ignore-case --pattern-file <(seqkit seq --seq --line-width 0 small.fa) --line-width 70 big.fa

You can use a data structure in Python called a dictionary, which lets you look things up by their name.

Here's a script that takes two linearized FASTA files reference.fa and other.fa, which makes a dictionary of sequence keys and header values from pairs of lines in reference.fa, and reads through other.fa to look up the reference dictionary and print out any matching sequence entries:

#!/usr/bin/env python

import sys

header = None
sequence = None

ref_seq_dict = dict()
ref_line_counter = 0    
with open('reference.fa') as ref_filehandle:
    for line in ref_filehandle:
        if ref_line_counter % 2 == 0:
            header = line.rstrip()
        else:
            sequence = line.rstrip()
            ref_seq_dict[sequence] = header
        ref_line_counter += 1

other_line_counter = 0
with open('other.fa') as other_filehandle:
    for line in other_filehandle:
        line = line.rstrip()
        if other_line_counter % 2 == 1:
            try:
                header = ref_seq_dict[line]
                sequence = line
                sys.stdout.write("%s\n%s\n" % (header, sequence))
            except KeyError:
                sys.stderr.write("Warning: Sequence in other file [%s] not found in reference!\n" % (sequence))
        other_line_counter += 1

You could use a similar approach in Perl, using a hash table.

Hash tables use a fair bit of memory, but I doubt that would be an issue for most computers these days. In a worst-case scenario where there's not enough memory, you could split up the reference file into smaller FASTA files to make smaller dictionaries, doing passes of the secondary files over each sub-dictionary.

This example script shows how you'd transform one file. If you have lots of files to transform, once you have built the dictionary, you can reuse it on as many files as you'd like. You might keep an array of filenames and loop through this array in the second half of the script, processing each file in turn. There are libraries in Python for parallelizing the work, though I/O is probably going to be the bottleneck, anyway.

Another approach that uses very little memory is to sort the files by their sequences. You can then stream through a pair of files one record at a time and move the file pointer forwards as matches are found. But you have to pre-process the files into a sortable form and that is a bit more work. You also usually lose the original order of headers after sorting, whereas hash table lookups preserve the ordering in the second file.

Do you have a specific reason why you don't use Biopython SeqIO to parse the fasta files? That would avoid assumptions about the format.

Log in to answer this question.