How to convert Foldseek .m8 alignment files to FASTA format?

Question

I want to use sequence alignments found by Foldseek in PROSS. I know there are other ways to get alignments (also tried with results of HHPred), but I had in mind that Foldseek may find remote homologs that would still be structurally relevant.

Currently Foldseek only gives alignment results in m8 format, however PROSS only accepts FASTA alignment files.

I tried looking into BioPython but it seems it doesn't open m8 files, only m9. I also tried importing the alignment in Clustal X but it can't read it.

I could probably use Foldseek CURL API to get the alignment as json and make a script to write the results in FASTA format, but this seems like a heavy handed solution.

What would be other tools I could try to convert m8 to FASTA?

Answer 1

tl;dr

Use an awk script (below) or "one-liner"

awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta

Caveat emptor: I've tested this on my own Foldseek outputs and it works for my pedestrian purposes of multisequence alignment with COBALT. I haven't tested it with PROSS but expect it should work for you too.

Long answer:

Using this definition of m8 (which is a modification of the BLAST m8 format):

https://github.com/steineggerlab/foldseek#output-search

And this definition of FASTA:

https://en.wikipedia.org/wiki/FASTA_format

We see we want to take lines like this (I've added [...] here to show what I've cut off for brevity):

job.pdb_A   Some matched structure name 51.100  311 138 8   2   305 25  328 1.000   2.790E-35   1272    306 334 VDVSVSIFINKIYGVNTLEQTYKVDGYIVAQWTGKPR--KTPGDKPLIVENTQIERWINNGLWVPALEFINVVGSPDTGNKRLMLFPDGRVIYNARFLGSFSNDMDFRLFPFDRQQFVLELEPFSYNNQQLRF--SDIQVYTENIDN[...] -27.488,-8.287,-35.993,-26.838,-4.500,-36.421,-25.278,-3.958,-32.938,-22.398,-6.436,-33.555,-21.375,-4.861,-36.892,-21.373,-1.326,-35.385,-19.237,-2.585,--18.044,-1.289,8.073,-14.345,-2.072,7.367,-15.099,-5.282,5.328,-16.872,-3.358,2.522,-16.876,-5.018,-0.77 [...] MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]   1918945 Vibrio spartinae

and map them to entries like:

>some matched structure name
MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]

We can do this with an awk script (The awk manual is your friend here https://www.gnu.org/software/gawk/manual/html_node/Long.html). Put the following in a file somewhere handy and call it using awk -f m8_to_fasta.awk my_m8_file.m8 > my_fasta_file.fasta you may want to combine this with head, tail, and/or grep to reduce the number of sequences that end up in your output.

#! /bin/awk    
{
        name = "> "
        # start on the second field, foldseek will put the job name in $1
        field = 2 
        # Check each field until we see numbers
        while ( $field + 0 != $field ) { 
           name = name " " $field
           field = field+1 
        }
        # name now contains the name string which is the first line of a fasta entry
        
        field = NF
        # Now find the matched sequence
        # Sometimes an entry ENDS in a number, so we strip that
        if ( $field + 0 == $field ) { field = field -1 }
        while ( $field + 0 != $field ) {
            field = field - 1
        }
        # field is now pointing at some undocumented number. Decrement once more to get the sequence
        field = field - 1
        # field is now pointing at the sequence
        seq = $field
        
        #Now I/we want to sanitize the sequence to remove any gaps
        gsub(/-/, "", seq)
        
        print name
        print seq
    }

Or we can coerce this into a monstrous one-liner:

 awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta