Perl script to convert raw biological strings in a file (DNA or amino acids strings) to FASTA format

Hello, this simple script is an update for another one I wrote eons ago.
Download the script <- seq2fas.pl

NOTE: It depends on the UNIX commands "pr","sed","tr" and "fold" to work out ;)

The functionality is very simple, it receives a file like this:

ACCTTACGCC
AGTAACGTAG
TTAGTATATA
ACCTACGATA
AAACAGGCCC
ACCGCTAGAT
AGCCCCATCC
CCGGTATACC
AGCGGACCCC
AACAACCCCC

And prints the output in FASTA format:

>Random_seq-1
ACCTTACGCC
>Random_seq-2
AGTAACGTAG
>Random_seq-3
TTAGTATATA
>Random_seq-4
ACCTACGATA
>Random_seq-5
AAACAGGCCC
>Random_seq-6
ACCGCTAGAT
>Random_seq-7
AGCCCCATCC
>Random_seq-8
CCGGTATACC
>Random_seq-9
AGCGGACCCC
>Random_seq-10
AACAACCCCC

This script is very useful if you are working with artificial sets produced in R, Perl, etc.

To run the script, just download it, change its permission to be executable and run it:

$ chmod +x seq2fas.pl
$ perl seq2fas.pl

Here is the USAGE:

perl seq2fas.pl randomDNAsequences.dna random_seqs 60 randomSeqs.fas

CODE:

#!/usr/bin/perl
 
################################################################################
# seq2fas.pl
# This script takes an input file like this:
#
#    ACCTTACGCC
#    AGTAACGTAG
#    TTAGTATATA
#    ACCTACGATA
#    AAACAGGCCC
#    ACCGCTAGAT
#    AGCCCCATCC
#    CCGGTATACC
#    AGCGGACCCC
#    AACAACCCCC
#
# and prints an output file in fasta format like this:
#
#    >1
#    ACCTTACGCC
#    >2
#    AGTAACGTAG
#    >3
#    TTAGTATATA
#    >4
#    ACCTACGATA
#    >5
#    AAACAGGCCC
#    >6
#    ACCGCTAGAT
#    >7
#    AGCCCCATCC
#    >8
#    CCGGTATACC
#    >9
#    AGCGGACCCC
#    >10
#    AACAACCCCC
#
################################################################################
# Author: Benjamin Tovar
################################################################################

use warnings;
use strict;

my $USAGE = "
USAGE:

seq2fas.pl    

EXAMPLE: 

seq2fas.pl randomDNAsequences.dna random_seqs 60 randomSeqs.fas

";

my $user_in = shift or die $USAGE;
my $fasta_header = shift or die $USAGE;
my $width = shift or die $USAGE;
my $output_file = shift or die $USAGE;

system("pr -n:3 -t -T $user_in | sed 's/^[ ]*/>$fasta_header-/' | tr \":\" \"\n\" | fold -w $width > $output_file");


exit;

Author: Benjamin Tovar

fasta2clustal.pl <- Perl script for convert aligned and gapless FASTA files to aligned and gapless CLUSTAL formated files

Fasta to Clustal:

Hello there, recently I am working with literally thousands of PWM (Positional Weight Matrices) generated from MEME and MotifClick output.

So, to use the motifs in FASTA format to generate PWM I need to align them first to generate CLUSTAL formated files but errors occur when gaps delivered from the alignment appear..

Yep, the UGENE PWM builder does not like gaps and shows up errors building the PWM.

The situation is that you already got an aligned and gapless FASTA file and you only need to convert it to CLUSTAL without "the need" to align the sequences once more and deal with gaps.

So, here is my Perl script called: fasta2clustal.pl

STEP 1 <- To execute it, just change the permission of the file to run as a program:

$ chmod +x fasta2clustal.pl

STEP 2 <- To run the program (you can copy and paste it inside your bin path or run the script locally):

# From the bin folder (The folder must contain the target FASTA files):

$ fasta2clustal.pl <input_fasta_file.fasta>

EXAMPLE: fasta2clustal.pl fasta_sequences.fasta

# From the local folder (Which contain the script and the target FASTA files)

$ ./fasta2clustal.pl <input_fasta_file.fasta>

EXAMPLE: ./fasta2clustal.pl fasta_sequences.fasta

SHORT TUTORIAL:

RUNNING THE SCRIPT

INPUT FILE:

OUTPUT FILE:

Code:

#!/usr/bin/perl


#       fasta2clustal.pl
#       
#       Copyright 2011 Benjamin Tovar
#       
#       This program is free software; you can redistribute it and/or modify
#       it under the terms of the GNU General Public License as published by
#       the Free Software Foundation; either version 2 of the License, or
#       (at your option) any later version.
#       
#       This program is distributed in the hope that it will be useful,
#       but WITHOUT ANY WARRANTY; without even the implied warranty of
#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#       GNU General Public License for more details.
#       
#       You should have received a copy of the GNU General Public License
#       along with this program; if not, write to the Free Software
#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#       MA 02110-1301, USA.
#
################################################################################
#
# DATE: 27/July/2011
# AUTHOR: Benjamin Tovar
#
# This program converts aligned and gapless FASTA file to
# a gapless CLUSTAL file, ideal for convert motifs to PWM.
#
################################################################################

use warnings;
use strict;

###### Print USAGE if no argument is given
my $usage = "\nUSAGE: fasta2clustal.pl 
EXAMPLE: fasta2clustal.pl dna_sequence.fa\n\n";

###### Get the arguments:  
my $input_fasta_file = shift or die $usage;

###### Let the computer decide:

 my $dna_sequence ='';
 my $sequence_name='';
 my $output_file_name='';

 open(INPUT_FILE,$input_fasta_file);
 my @file = ;
 close INPUT_FILE;

  # ClUSTAL HEADER
  print "CLUSTAL W 2.0 multiple sequence alignment\n\n";

 foreach my $line(@file){
  
  # Discard empty lines
  if($line =~ /^\s*$/){
  next;
  }
  
  else{
    # Extract sequence names
    if($line =~ /^>/) {
     $sequence_name = $line;
     $sequence_name =~ s/\s//g;
     $sequence_name =~ s/\>//g;
     $sequence_name =~ s/\Start_position/_Start_position/g;
     $sequence_name =~ s/\:/_/g;
     chomp $sequence_name;
    }

    # Extract DNA sequence
    else{
    $dna_sequence = $line;
    chomp $dna_sequence;
  # Print results
                printf("%-50s %10s \n", "$sequence_name","$dna_sequence");
    } 
  }
 }

# Powred by #!CrunchBang Linux
# Benjamin Tovar

exit;

Benjamin

Perl script to delete repeated entries of a plain text file using the Linux/UNIX command "awk"

This script called "delete_repeats.pl" will help you if you are looking for a more user friendly alternative to execute the following lines inside a Terminal:

$ awk ' !x[$0]++' input_file > output_file

I know! the syntax is very logic, simple and probably is not necessary to deserve a Perl script to automatize this process of deleting entries with this simple line of code.

But, the good thing about writing a script that ask me for the input and output and then automatically replace those values inside the system command awd is that I no longer need to type the whole line.

I copied my script inside my "bin" folder, execute it and then I just worry about type the correct name of the input file and type the output file name and that is all ;)

Here is an example:

Input file:

Now, lets execute the Perl script with the following line inside a terminal (remember that we must go inside the folder that contains our input file first and be sure that our script is inside the same folder):


$ perl delete_repeats.pl


NOTE: I execute the script this way:


$ delete_repeats.pl


Because I copied the script to my /home/benjamin/bin folder so is no longer necessary to copy it script to every folder where I want to execute it. This way lets the Terminal to recognize it and execute it inside every folder that I want to use it :D

And the output :D

Benjamin

Bug fixed: Perl script to generate "n" random DNA strings of "n" length with a desired nucleotide probability distribution defined by the user and export them in a file in FASTA format

After a short review of my own script "random_dna_strings.pl" I realized that is more useful to set the parameter of C content distribution probability (which is the last parameter of probability distribution that is asked to the user) automatically to:

line 116:

my $C_content = (1-($A_content+$T_content+$G_content))

Rather than ask the user for typing manually the remaining value that finally satisfy the 1.00 of total probability value.

Total p value = 1 = (A content + T content + G content + C content)


C content = (1 - (A content + T content + G content))

Benjamin

Perl script to generate "n" random DNA strings of "n" length with a desired nucleotide probability distribution defined by the user and export them in a file in FASTA format

How many of you that are working with artificial DNA sets and have to use the tools that are available over the Web have to edit the output because the FASTA header of such programs does not correspond to your needs, the output is not in FASTA format or the lazy reason to copy and paste the output sequences inside a manually generated file because it depends on important ATP molecules that can be used in another tasks.

Well, this Perl script called "random_dna_strings.pl" can be downloaded by clicking the link.

This script takes arguments given by the user such the nucleotide frequencies of each nucleotide (in a scale from 0.0 to 1.0), generates a "n" number of sequences of "n" length with a FASTA header also given by the user and finally prints an output file in FASTA format.

So, here is a short Tutorial about how to use it:

Well, imagine that we want to generate 100 random DNA sequences of 100 nucleotides length with a nucleotide distribution equal to every nucleotide, this is "A=0.25, T=0.25, G=0.25 and C=0.25" with the FASTA header of "just a bunch of random sequences of equal random prob nt distribution" in a file called "ran_set_1.fa".

STEP 1 <- Open a Terminal inside the output folder (where you want to put the output file)

STEP 2 <- Execute the Perl script (This time, I copied the script to my bin folder, if you have questions about how to run your scripts this way, please visit this entry.

Command:


$ random_dna_strings.pl

STEP 3 <- Please type the number of iterations (How many random sequences do you want)

In this Tutorial I want to generate 100 artificial sequences.

STEP 4 <- Please type the length of the random DNA strings (how many nucleotides length)

In this Tutorial I want each sequence to be 100 nucleotides long

STEP 5 <- Please type the probability distribution of A, T, G and C content:

In this Tutorial I want each nucleotide got an equal probability distribution (this means 1/4 of probability of A,T,G or C at each string position)

STEP 6 <- Please, type the name of the fasta header for each sequence (is not necessary to put the >)

In this Tutorial, I want that the FASTA header be: "just a bunch of random sequences of equal random prob nt distribution"

STEP 7 <- Please, type the name of the output file:

In this Tutorial, I want to name it "ran_set_1.fa"

STEP 8 <- Enjoy the results

Another example:

Now I want to generate 100 random DNA sequences of 100 nucleotides length with a nucleotide distribution of "A=0.7, T=0.2, G=0.1 and C=0.0" with the FASTA header of "random dna sequence with too many As" in a file called "random_set_A_rich.fa".

Now, take a look of the output folder:

Do you see it :D, now we got two output files in FASTA format of random made sequences that follows a probability distribution of nucleotides given by us.

NOTE: The yellow highlighting just mean the differences of the content of A nucleotides among the two files showing us that the probability distribution obviously has an impact on the nucleotide composition of the random generated DNA strings.

Benjamin

Perl script to export every line of independent DNA sequences inside a single file to a FASTA formated output file

This Perl script will take every line of the input file as a single DNA sequence and then it will create a new file with every sequence in FASTA format with the FASTA headers defined by the user.

The script is called "fasta_seq.pl" and can be downloaded by clicking the link.

So here it is a short Tutorial about how to use it:

Imagine that we got a file called "seq" that contains a single DNA sequence per row and we do not want to separate each one with spaces and name them manually.

NOTE: This script only will run under Linux/UNIX environments because it depends of the commands: "pr","sed","tr" and "fold" (Sorry to all the Windows users).

STEP 1 <- Open a Terminal and get inside the input files folder:

STEP 2 <- Execute the Perl script with

perl fasta_seq.fa


EXPLANATION OF THE PARAMETERS:

1. perl <- here you are telling the Terminal that you want to run the Perl environment
2. fasta_seq.fa <- name of the Perl script

STEP 3 <- Type the name of the input file that contains the DNA sequences (in this Tutorial, our file is called "seq")

STEP 4 <- Type the name of the fasta header for each sequence (is not necessary to put the ">" symbol).

In this Tutorial, I want that the FASTA header will be "just a dna sequence"

STEP 5 <- Type the width of the sequences (how many nucleotides per column)

In this Tutorial, I want that the width of nucleotides per column will be "45" nucleotides

STEP 6 <- Type the complete name of output file (In this Tutorial, I want that the output file name of my results will be "output_sequences.fa"

STEP 7 <- Finally, enjoy the results ;)

Now we can go into our input folder and take a look of the output file:

Do you see it?? successfully now we got an output file that contains every sequence of the input file (this script considers every line as an independent DNA string) under a fasta header named by ourselves and ready in FASTA format.


Benjamin

Simple Perl script for computing the reverse complement of a single DNA string

This entry will help you if you are looking for a Perl script that reads an input file (it can be in FASTA format or just the sequence without the FASTA header) that contains just one sequence (at least for now, later on I will upload a script that can handle n number of sequences inside a single file in FASTA format), computes the reverse complement of the string and then create an output file in FASTA format (even if the original file was not in FASTA format).

The script is named "rev_com.pl" and can be downloaded by clicking in the script name.

NOTE: In a Linux/UNIX system we can set the properties of the script to be executable and then put it inside the "bin" folder (like /usr/local/bin or /home/USER/bin), So in this way we can run the script in every directory we want without the need to copy it to every folder where we want to run the script.

Here is an example of my personal "bin" folder:

But in this short Tutorial I will run the Perl script that is inside the folder of my DNA sequences and not the one that is inside my "bin" folder.



Lets run an example to see how it works ;)

Lets imagine that the got this two files inside a folder:

1. One is in FASTA format <- D_mel.fa 
2. The other one just contains the DNA sequence <- dna_sequence

And we want to compute the reverse complement of each file and export them in a new output file in FASTA format with the FASTA header.

So, lets run the script:

NOTE: remember that we got to access the folder with the Terminal to execute the script. We can do this by using the UNIX command "cd" and finally get access to the folder that contains our files.

1. STEP 1 <- Access the folder with the terminal
2. STEP 2 <- Execute the Perl script:


$ perl rev_com.pl


Example of executing the Perl script

EXPLANATION OF THE PARAMETERS

perl <- You are telling the Terminal that you want to run a Perl script.

rev_com.pl <- name of the perl script.

STEP 3 <- Write the name of the input file (if the file has extension, do not forget to type it too):

STEP 4 <- Hit enter and run! (it was a joke):

Enjoy the results directly from the Terminal

In the bottom of the Terminal, you can read how is the name of the output file :D

Enjoy the results of another example of an input file that is not in FASTA format (it just contains the DNA sequence pasted in)

Finally, we can go to the folder and see our new files that contains the reverse complement of an input user string:

Can you see it now? The script had named the output files with the ".fa" extension and inside the file, we can see the FASTA headers in making them officially FASTA readable and recognizable files.

##############################################

If you are curious and want to know how to run the script from the "bin" folder, here is an screenshot about it:

As you can see, there is no need to tell the Terminal that we want to run a Perl script and another good thing about running programs from the bin folder is that we do not need to copy them to every input folder (where our input files lives in).

#!/usr/bin/perl
#       rev_com.pl
#       
#       Copyright 2011 Benjamin Tovar 
#       
#       This program is free software; you can redistribute it and/or modify
#       it under the terms of the GNU General Public License as published by
#       the Free Software Foundation; either version 2 of the License, or
#       (at your option) any later version.
#       
#       This program is distributed in the hope that it will be useful,
#       but WITHOUT ANY WARRANTY; without even the implied warranty of
#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#       GNU General Public License for more details.
#       
#       You should have received a copy of the GNU General Public License
#       along with this program; if not, write to the Free Software
#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#       MA 02110-1301, USA.
#
################################################################################
#
# NAME OF THE PROGRAM: rev_com.pl
# DATE: 01/Jun/2011
# AUTHOR: Benjamin Tovar
# COMMENTS: Little script that can be used for the computation of 
# the reverse complement of a DNA string and then export the resulting string
# to a file in Fasta format with fasta headers (even if the original file was not)
#
################################################################################

use warnings;
use strict;

#### script introduction

print "\n\n-------------------------------------------------------------\n\n";

print "PROGRAM DEFINITION: rev_com.pl <- PERL SCRIPT WRITTEN BY BENJAMIN TOVAR
THAT COMPUTES THE REVERSE COMPLEMENT OF DNA STRINGS AND EXPORT THE RESULTS TO AN 
OUTPUT FILE IN FASTA FORMAT (EVEN IF THE ORIGINAL FILE WAS NOT)\n\n";

print "This program will ask the user to type the name of the file 
(if the file name has an extension, write it please) that contains the DNA sequence\n\n";

#### Ask the user for the name of the file

print "Please, type the complete name of the file:\n\n";

my $user_in =;

### Remove empty spaces

chomp $user_in;

### Open the file

open(INPUT_FILE,$user_in) or die "\n\n WARNING: The file does not exist, please check the spelling, the extension and the existence of the file\n\n";

### Copy the content of the file to an array variable called "@input_file"

my @input_file = ;

### Remove the Fasta header and extract the DNA sequence

my $input_file = extract_sequence_from_fasta_data(@input_file);

### Close the opened file

close INPUT_FILE;

### Compute the reverse DNA string

my $rev_com = reverse ($input_file);

### Compute the complement DNA string

$rev_com =~ tr/ATGCatgc/TACGtacg/;

### Compute the length of the input DNA string

my $length = length ($input_file);

############## RESULTS #################################################

print "\n----------- RESULTS --------------\n\n";

print "Input file string of length $length:\n\n";

print $input_file,"\n";

print "\nOutput reverse complement of the DNA string:\n\n";

print $rev_com, "\n\n";

############## EXPORTING THE RESULTS TO A FILE IN FASTA FORMAT #########

######################## Naming the output file: #######################

### sdfgh

my $fasta_header_name = extract_fasta_header_name(@input_file);

## Remove the input file extension

$user_in =~ s/\..*//; 

# Concatenate the name of the file (without the .fa/.fasta extension) with "-rev_com.fa"

my $rev_com_name = "-revcom.fa";

my $output_name = $user_in . $rev_com_name;

# Name of the output file

my $out = "$output_name";

# Set the file handle "OUTPUT".

open (OUTPUT, ">$out"); 

# Print the results (content) of the variable "$rev_com" (this variable contains the reverse complement string) 
# into a file named "$output_name" and put the Fasta header before the output DNA string with "$fasta_header_name","-reverse_complementary\n"

print OUTPUT "$fasta_header_name","-reverse_complement\n","$rev_com";


print "-------- EXPORT THE RESULTS TO A FILE IN FASTA FORMAT ----------\n";
print "\nThe output string has been exported to the file \"$output_name\"\n\n";

exit;

################################################################################
############################### SUBROUTINES ####################################
################################################################################


################################################################################
# extract_fasta_header_name
# A subroutine to extract the FASTA header of the original input file
# and use it to name the FASTA header of the output file
################################################################################

sub extract_fasta_header_name{

    my(@fasta_file_data) = @_;

    use strict;
    use warnings;

    # Declare and initialize variables
    my $fasta_header_name = '';

    foreach my $line (@fasta_file_data) {
    
        if($line =~ /^>/) {
            
            $fasta_header_name = $line;
    
    # If the file is not in Fasta format, use the name of the file to name the fasta header           
        } else {
        
            $fasta_header_name = ">$user_in";
        
        }
        
    # Remove non-sequence data (in this case, whitespace) from $fasta_header_name string
        $fasta_header_name =~ s/\s//g;
        
    # Export the results of the subroutine to the main program    
        return $fasta_header_name;
        
    }
}

################################################################################
# extract_sequence_from_fasta_data
# A subroutine to extract FASTA sequence data from an array
# taken from James Tisdall's Beginning Perl for Bioinformatics
################################################################################

sub extract_sequence_from_fasta_data {

    my(@fasta_file_data) = @_;

    use strict;
    use warnings;

    # Declare and initialize variables
    my $sequence = '';

    foreach my $line (@fasta_file_data) {

        # discard blank line
        if ($line =~ /^\s*$/) {
            next;

        # discard comment line
        } elsif($line =~ /^\s*#/) {
            next;

        # discard fasta header line
        } elsif($line =~ /^>/) {
            next;

        # keep line, add to sequence string
        } else {
            $sequence .= $line;
        }
    }

    # remove non-sequence data (in this case, whitespace) from $sequence string
    $sequence =~ s/\s//g;

    return $sequence;
}

Benjamin.