This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Make multi-fasta using perl

Header and sequence are coming in loop:

$header = "seq1";

$sequence = "GTGGTCAGAAAAGTGGCCTCAGCCCATTTCCAGCATCCTATGGCTATCTAAGTAGGAAGACGATGAGCTATGGAA";

Sequence is single line without newline.

I want output multi-fasta like this below. Sequence should have newline after each 60 bases.

>seq1
GTGGTCAGAAAAGTGGCCTCAGCCCATTTCCAG
CATCCTATGGCTATCTAAGTAGGAAGACGATGA
GCTATGGAA

>seq2
AAGTGGCCTCAGCCCATTTCCAGCATCCTATGGC
TATCTAAGTAGGAAGACGATGAGCTATGGAA

Thanks

perl multi-fasta fasta
• 6,165 views • link
updated by Ram Ram • written by sun.nation sun.nation

Can you please give me some hint to create loop in perl? I need that at the end of my perl code.

Thanks Ram and Dylan

• 0 views • link
updated by Ram Ram • written by sun.nation sun.nation

Look at the Text::Wrap module if you need a perl specific and easy method.

• 0 views • link
updated by Ram Ram • written by dylan.storey dylan.storey

Thanks. I will try that

• 0 views • link
written by sun.nation sun.nation

This was asked before here and the solution (pre-v5.10) is:

$sequence =~ s/(.{60})/$1\n/gs;

That formats the sequence, then you just need to join the id and sequence with a newline when you print.

• 0 views • link
updated by Ram Ram • written by SES SES

Thanks a lot SES.

It was a fantastic solution with a line code. It worked good when I added this in my script.

• 0 views • link
written by sun.nation sun.nation

I think the following lines could be an alternative solution.

print ">$header\n";
$seq_len = length($sequence);

for ($i = 0; $i < $seq_len; $i + 60) {
    $segment = substr($sequence, $i, 60);
    print "$segment\n";
}
• 0 views • link
updated by Ram Ram • written by biolab biolab

Yes, the loop is a good basic alternative, but the RegEx is better, this being Perl.

• 0 views • link
written by Ram Ram

The problem with the regex is that it is limited to a max sequence length of ~32k, which gives you trouble with contigs, etc. Here is a loop solution that should be faster than regex.

my $single_line_seq = "ATTGC"x1000;
my $line_width = 80;
print $_,"\n" for unpack("(A$line_width)*", $single_line_seq);
• 0 views • link
updated by Ram Ram • written by thackl thackl

3 answers

Here is a comparison of the two methods:

Here are the results:

$ perl biostars140255.pl 2>/dev/null
Benchmark: timing 10 iterations of regex, unpack...
     regex: 17 wallclock secs (15.48 usr +  1.44 sys = 16.92 CPU) @  0.59/s (n=10)
    unpack: 66 wallclock secs (29.63 usr + 34.99 sys = 64.62 CPU) @  0.15/s (n=10)
       s/iter unpack  regex
unpack   6.46     --   -74%
regex    1.69   282%     --

The benchmark shows that regex is about 6X faster than unpack. Also, there is no length limit (the input in the script is a 100 Mb string). Are you referring to a 32-bit system with those limits?

• 0 views • link
updated by Ram Ram • written by SES SES

I had more than 200k in a single fasta, regex worked good. Thanks

• 0 views • link
written by sun.nation sun.nation

How I see SES now:

image: dwight meme

• 0 views • link
written by Ram Ram

I was looking for a mic drop gif but just gave up myself.

• 0 views • link
written by dylan.storey dylan.storey

Ha! I might have been a bit terse but I honestly wasn't sure until I did the benchmark so I wanted to share. Also, I use 5.20 in the script because there is an awesome optimization (called "copy on write") that is lightening fast for copying strings. Here is an example script:

my $x = "A" x 1_000_000;
my $y = $x for 1...1_000_000;

That makes a long string and copies it a bunch of times in memory. Super useful for genome analysis. With Perl version 5.10.1 that takes 250.8 seconds on my machine. With Perl 5.20.2, that takes 0.15 seconds!

• 0 views • link
updated by Ram Ram • written by SES SES

Okay, my bad - I shouldn't have mentioned speed in the first place. unpack can be faster, especially if I'm not too lazy to put print outside the loop. I am also sure, one could optimize things further...

'unpack' => sub {
    say STDERR join("\n", unpack("(A60)*", $seq));
}

         Rate  regex unpack
regex  1.74/s     --   -47%
unpack 3.27/s    88%     --

Regarding the limit of 32766, I got confused there. It only applies to a maximum size of the capture group and hasn't any effect here, unless you want your lines to be longer than 32 kbp ;)

• 0 views • link
updated by Ram Ram • written by thackl thackl

Either use a loop and introduce a new line characters after every 60 characters in the sequence, or use Picard's NormalizeFasta if you'd like a tried and tested method.

• 0 views • link
written by Ram Ram

You can also try cat *.fasta | fold > all.fasta if you know that the headers are less than 60 bp.

• 0 views • link
updated by Ram Ram • written by dylan.storey dylan.storey

Can you please give me some hint to create loop in perl? I need that at the end of my perl code.

• 0 views • link
updated by Ram Ram • written by sun.nation sun.nation 
fasta_formatter -i input.fasta -o output.fasta -w 60

fasta_formatter: http://hannonlab.cshl.edu/fastx_toolkit/commandline.html#fasta_formatter_usage

• 0 views • link
updated by Ram Ram • written by Prakki Rama Prakki Rama

Log in to answer this question.

More posts like this