This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Merging of sequence chunks

Hello!

I would like to merge all amino acid sequences with identical names in a given fasta file (trimmed alignment). I know that this can be done for manually selected sequences with Aliview, however, I am looking for a tool or script (bash, perl, python) that allows me to do that for hundreds of files without manual selection. Non-matching amino acids should be replaced by "X". If the script/tool allows to set up a sequence similarity and overlapping length threshold, even better :) Any help would be much appreciated. Thank you!

My file:

>Blue
MADKLTRIAIVNHDKCKPKKCRQECKKSCPVVRMGKLCIEVTPQSKIAWISETLCIGCGI
>Red
CIKKCPFGALSIVNLPSNLEKETTHRYCAI------------------------------
>Red
-----------------------THRYCANAFKLHRLPIPRPGEVL--------------
>Red
--------------------------------------------------TNGIGKSTAL
>Yellow
KILAGKQKPNLGKYDDPPDWQEILTYFRGSELQNYFTKILEDDLKAIIKPQYVDQIPKAA
>Green
KGTVGSILDRKDETKTQAI-----------------------------------------
>Green
-------------TKKQAIVCQQLDLTHLKERNVEDLSGGELQRFACAVVCIQDQICKKI

Desired output:

>Blue
MADKLTRIAIVNHDKCKPKKCRQECKKSCPVVRMGKLCIEVTPQSKIAWISETLCIGCGI
>Red
CIKKCPFGALSIVNLPSNLEKETTHRYCAXAFKLHRLPIPRPGEVL----TNGIGKSTAL                                                          
>Yellow
KILAGKQKPNLGKYDDPPDWQEILTYFRGSELQNYFTKILEDDLKAIIKPQYVDQIPKAA
>Green
KGTVGSILDRKDETKXQAIVCQQLDLTHLKERNVEDLSGGELQRFACAVVCIQDQICKKI
alignment sequence
• 3,245 views • link
written by jan • jan

jan : Please test the perl solution offered by @JC. Then accept it as well (you can accept multiple answers) if all looks well.

• 0 views • link
written by GenoMax GenoMax

2 answers

using bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

and the following script:

final Map<String,List<FastaSequence>> id2seqs= new HashMap<>();
stream().forEach(S->{
    String n = S.getName();
    if(!id2seqs.containsKey(n)) id2seqs.put(n,new ArrayList<>());
    id2seqs.get(n).add(S);
    });

for(final List<FastaSequence> seqs: id2seqs.values())
{
println(">"+seqs.get(0).getName());
int len = seqs.stream().mapToInt(S->S.length()).max().getAsInt();
for(int i=0;i< len;i++)
    {
    char c='\0';
    for(FastaSequence s:seqs)
        {
        if(i>=s.length() || s.charAt(i)=='-'  ||  c==s.charAt(i)) continue;
        if(c=='\0')
            {
            c=s.charAt(i);
            }
        else
            {   
            c='X';
            }
        }
    print(c=='\0'?'-':c);
    }
println();
}

usage:

$ java -jar  dist/bioalcidaejdk.jar -f script.txt in.fa

output:

>Red
CIKKCPFGALSIVNLPSNLEKETTHRYCAXAFKLHRLPIPRPGEVL----TNGIGKSTAL
>Blue
MADKLTRIAIVNHDKCKPKKCRQECKKSCPVVRMGKLCIEVTPQSKIAWISETLCIGCGI
>Yellow
KILAGKQKPNLGKYDDPPDWQEILTYFRGSELQNYFTKILEDDLKAIIKPQYVDQIPKAA
>Green
KGTVGSILDRKDETKXQAIVCQQLDLTHLKERNVEDLSGGELQRFACAVVCIQDQICKKI
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks Pierre! Unfortunately, however, I don't seem to be able to install BioAlcidaejdk (there is no bioalcidaejdk.jar in dist/). Can you maybe provide some guidelines how to install this? Is it also running on a Mac?

• 0 views • link
written by jan • jan

what happened after you typed :

make bioalcidaejdk
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Sorry, my fault. After running sudo apt install default-jdk, it worked. The tool is just great! Thank you very much :-) Could the script also be modified in a way that it accepts headers that are identical only before the "@" sign (e.g., >red@abcd and >red@efgh)?

• 0 views • link
written by jan • jan

change String n = S.getName();

to

String n = S.getName();
int a=n.indexOf('@');
if(a!=-1) n=n.substring(0,a);
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Perl option:

#!/usr/bin/perl

use strict;
use warnings;

my %seqs = ();

# load sequences
$/ = "\n>"; # ingest all fasta sequences
    while (<>) {
    s/>//g;
    my ($name, @seq) = split (/\n/, $_);
    #$name =~ s/@.+//; # uncomment to edit names like red@abc, red@edf
    my $seq = join ("", @seq);
    push @{ $seqs{$name} }, $seq;
}

# process all sequences
foreach my $seq (keys %seqs) {
    print ">$seq\n";
    my $nseqs = scalar @{ $seqs{$seq} };
    if ($nseqs == 1) { # for just one sequence
        print "$seqs{$seq}[0]\n";
    }
    else {
        my $len = length $seqs{$seq}[0];
        for (my $i=0; $i<=$len-1; $i++) {
            my %b = ();
            for (my $j=0; $j<=$nseqs-1; $j++) {
                $b{ substr ($seqs{$seq}[$j], $i, 1) }++;
            }
            my $obs = join ("", keys %b);
            $obs =~ s/-//;
            if    (length $obs == 1) { print "$obs"; }
            elsif (length $obs == 0) { print    "-"; }
            else                     { print    "X"; }
        }
        print "\n";
    }
}

Usage:

$ perl mergeSeqs.pl < test_seqs.fa
>Red
CIKKCPFGALSIVNLPSNLEKETTHRYCAXAFKLHRLPIPRPGEVL----TNGIGKSTAL
>Green
KGTVGSILDRKDETKXQAIVCQQLDLTHLKERNVEDLSGGELQRFACAVVCIQDQICKKI
>Blue
MADKLTRIAIVNHDKCKPKKCRQECKKSCPVVRMGKLCIEVTPQSKIAWISETLCIGCGI
>Yellow
KILAGKQKPNLGKYDDPPDWQEILTYFRGSELQNYFTKILEDDLKAIIKPQYVDQIPKAA
• 0 views • link
written by JC JC

Thanks JC, the script is doing the job!

• 0 views • link
written by jan • jan

Log in to answer this question.

More posts like this