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

Hi, My protein of interest is a multiple domain protein, and I would like to prepare a "structure based multiple sequence alignment" of it (around 1000 sequences).

Since the available structures of my protein are in different conformational states, I had to make multiple sequence alignments for EACH DOMAIN SEPARATELY. Now I have 6 high quality alignments, corresponding to the six domains of my protein.

The question is that: How can I merge/attach/join these files together? Notes:

  • I am new to programming!
  • All files are in Fasta format (also i have them in Clustal format)
  • The length of ALL sequences in each file is the same (but the length of sequences in different files are different. For example in File 1, the length of all sequences is 120, and in file 2 all of them are 78)
  • Sequences are arranged base on the "ID" in these files (i.e. the sequence IDs in the files are the same & are in the same order)
  • I've used cat command, but it just added sequences at the end of file!
  • Example files are presented below (I have shown only the first 4 sequences. there are around 1000 sequences). they are not real sequences.

I appreciate your help regarding to this matter, and please let me know if you need any more information.

Cheers, Ali

#

#

File 1:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
MLGWIAKKII
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
MLGWIAKKII
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
MLGWIAKKII
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
MLGWIAKKII

File 2:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
--GWIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
--GWIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
--GWIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
--GWIAKKI-AAAA

File 3:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
--GGYYYKLIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
--GGYYYKLIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
--GGYYYKLIAKKI-AAAA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
--GGYYYKLIAKKI-AAAA

File 4:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
--------HHWSGYYYKLIAKKI-AAAA--D
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
--------HHWSGYYYKLIAKKI-AAAA--D
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
--------HHWSGYYYKLIAKKI-AAAA--D
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
--------HHWSGYYYKLIAKKI-AAAA--D

File 5:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
AEE
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
EEE
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
EEE
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
EEE

File 6:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
A
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
A
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
A
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
A

So the final file should look like this: Result file:

>gi|6226123|sp|O67718.1|SECA_AQUAE/1-88 RecName: Full=Protein translocase subunit SecA
MLGWIAKKII--GWIAKKI-AAAA--GGYYYKLIAKKI-AAAA--------HHWSGYYYKLIAKKI-AAAA--DAEEA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA1
MLGWIAKKII--GWIAKKI-AAAA--GGYYYKLIAKKI-AAAA--------HHWSGYYYKLIAKKI-AAAA--DEEEA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA3
MLGWIAKKII--GWIAKKI-AAAA--GGYYYKLIAKKI-AAAA--------HHWSGYYYKLIAKKI-AAAA--DEEEA
>gi|81837980|sp|Q821L5.1|SECA_CHLCV/1-78 RecName: Full=Protein translocase subunit SecA4
MLGWIAKKII--GWIAKKI-AAAA--GGYYYKLIAKKI-AAAA--------HHWSGYYYKLIAKKI-AAAA--DEEEA
multiple-alignment
• 9,824 views • link
updated by User000 User000 • written by akhalili26 • akhalili26

do you want to merge the sequences together or files? basically, it is the same ID?

• 0 views • link
written by User000 User000

show us a few lines of your fasta files please...

• 226 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Updated my question. please let me know if you need any more information.

• 0 views • link
written by akhalili26 • akhalili26

That looks like CLUSTAL format, or something, but is not at all fasta.

• 247 views • link
written by pld pld

Does all of your 1000 sequences have their name ending with "SecA"?

• 0 views • link
written by Ashutosh Pandey Ashutosh Pandey

Not necessarily, they have they own Unique name "ID". at the beginning of each "name" or ID, there is ">" character (which is not shown here! something wrong with the forum)

• 0 views • link
written by akhalili26 • akhalili26

Next time, it is better to past sequences as code rather than quotes. I've updated my post to pull the gi number from your sequence and concatenate the alignments on that. Give it a try.

• 28 views • link
written by pld pld

Hey thank you very much. I am new to scripting, specially Python. can you please tell me how can I insert multiple files (six files) as input file here?

• 0 views • link
written by akhalili26 • akhalili26

Check my post, I've added a way to do this.

• 0 views • link
written by pld pld

I tried your script. and ran this command:

python script3.py allfasta.fa new.fa

I got this error:

File "script3.py", line 4

infile = sys.argv[1]
     ^
IndentationError: expected an indented block
• 0 views • link
written by akhalili26 • akhalili26

Oops, seems like I goofed up the code when putting it into my post. I just need to add the spaces back again.

• 229 views • link
written by pld pld

Thanks, i got another error:

Traceback (most recent call last):

File "script4.py", line 30, in <module>
  __main__()
File "script4.py", line 16, in __main__
  key = re.search('(?<=gi|)([0-9]*)(?=|)', x[0]).group(0)
File "/sw/lib/python2.7/re.py", line 142, in search
  return _compile(pattern, flags).search(string)
File "/sw/lib/python2.7/re.py", line 242, in _compile
  raise error, v # invalid expression
sre_constants.error: look-behind requires fixed-width pattern
• 0 views • link
written by akhalili26 • akhalili26

Forgot to escape the pipes when I updated it, it works now.

• 0 views • link
written by pld pld

4 answers

Here is a python script that will concat your alignments, assuming each alignment has the same id in the fasta header

import sys
import re
def __main__():
    infile  = sys.argv[1]
    outfile = sys.argv[2]

    with open(infile, 'rb') as fi:
        seqs = fi.read().split('>')[1:]
        seqs = [x.split('\n')[:2] for x in seqs]

    merge = dict()

    #index 0 is header
    #index 1 is seq
    for x in seqs:
        key = re.search('(?<=gi\|)([0-9]*)(?=\|)', x[0]).group(0)
        try:
            merge[key] = merge[key] + x[1]

        except KeyError:
            merge[key] = x[1]


    #write results to file
    with open(outfile, 'w') as fi:
       for x in merge.keys():
           fi.write('>' + x + '\n')
           fi.write(merge[x] + '\n')

__main__()

To put your individual files into a single file:

cat f1.fasta, f2.fasta, ... > bigfile.fasta

You would then run the above script with:

python myscript.py bigfile.fasta outputfile.fasta

Just remember to adjust the file names here to whatever you're working with.

• 273 views • link
written by pld pld

Sorry for the errors in the code, I verified the regex but typed it into the post, hence the errors. Everything works now.

• 0 views • link
written by pld pld

If I understand correctly, you have six files (one for each domain) and each file having 1000 sequences.

1) Excluding the clustalw file for the first domain, run the following command for all the other five files. Make sure you have backed up your files so if something goes wrong, your original files remain unchanged.

sed -i ';s/.* //' File_name_2

sed -i ';s/.* //' File_name_3

........

sed -i ';s/.* //' File_name_6

2) Then use join command in unix and concatenate the first file with the second , third .. sixth file.

join File_name_1 File_name_2 ... File_name_6 > Mega_file

Hope it works.

• 233 views • link
written by Ashutosh Pandey Ashutosh Pandey

I'm new to scripting! I got this when running sed -i ';s/.* //' domain2_cut_clustal

sed: 1: "domain2_cut_clustal": extra characters at the end of d command

• 0 views • link
written by akhalili26 • akhalili26

There should be a gap between "'" and domain_cut_clustal

• 0 views • link
written by Ashutosh Pandey Ashutosh Pandey

Thanks for your help. It does the job of combining the files. BUT, When I tried to use the final "Mega_file", i realized that it is not a Clustal file! since in the middle of the file some of the lines don't meet the Clustal's standard. (i. e. some of the lines have less than 60 sequence symbols. This is because, for example, the last lines of file 2, which are in the middle of the "Mega_file", have 24 sequence symbols (which i believe will be interpreted as the end of file!))

any suggestion for turning this "Mega file" into a Fasta file or Clustal file?

• 0 views • link
written by akhalili26 • akhalili26

Essentially, what you want to do regardless of programming language is:

  1. Parse a set of "fasta-like" files (they are FASTA multiple alignment format which for all intents and purposes is fasta)
  2. Extract the ID, description and sequence
  3. For each unique ID + description, concatenate the associated sequences
  4. Print results

So here is a quick, dirty and ugly solution which employs BioPerl SeqIO. I'm assuming a Linux-like system with a command line. BioPerl installation on Ubuntu-based Linux is as simple as sudo apt-get install bioperl.

First, assuming that your sequence files are the only files in their directory and they end with the suffix ".fa", concatenate them into one file:

cat *.fa > allseqs.fa

Now the BioPerl; save this as merge_seqs.pl:

#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;

my %seqs; # a hash to hold the results
my $file = shift or die("Usage = merge_seqs.pl <fasta file>\n");
my $seqio = Bio::SeqIO->new(-file => $file, -format => "fasta");

while(my $seq = $seqio->next_seq) { 
  $seqs{$seq->display_id." ".$seq->description} .= $seq->seq;
}
# now print
foreach my $key (keys %seqs) { 
  print ">$key\n$seqs{$key}\n";
}

and run like so (output to the file merged.fa):

perl merge_seqs.pl allseqs.fa > merged.fa
• 0 views • link
written by Neilfws Neilfws

This ought to do it:

unix$ cat file* | perl -pe 's/\r\n?/\n/' | perl -ne 'chomp;if(/^>/) { $id = $_ } else { $content{$id}.=$_} END{print map { "$id\n$content{$id}\n" } keys %content}'  > file.all

It will work even if the order of the sequences inside the files are flipped around or if the lengths are different inside the files. It requires all files to be small enough to fit in memory.

Edit:

Modified to accept Apple \r and Microsoft \r\n line separation, too.

• 0 views • link
written by ole.tange ole.tange

I ran your script. It took the last sequence of the last file and repeated it over 3000 times in the "file.all" file!

• 0 views • link
written by akhalili26 • akhalili26

I have re-tested it on the 6 files you give, and I still get the correct result. Is there something very different hiding in your other files? Can you reproduce my results by running it on your 6 example files?

• 0 views • link
written by ole.tange ole.tange

Were these files ever saved on Windows? That might be the issue.

• 31 views • link
written by pld pld

That might just explain the results. Edited answer to support both Windows and Mac. Does it work now?

• 0 views • link
written by ole.tange ole.tange

Log in to answer this question.

More posts like this