Trim multifasta.pl: Difference between revisions

From RAJ INFO
Jump to navigation Jump to search

Raj (talk | contribs)

2 user groups
11,120 edits
VisualWikitext
Created page with 'Category: Bioinformatics <source lang=perl> #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in seque…'
 
No edit summary
 
Line 3: Line 3:
<source lang=perl>
<source lang=perl>
#! /usr/bin/perl -w
#! /usr/bin/perl -w
# *************************************************************
# *************************************************************
# This Programme does the following in sequence
# This Programme does the following in sequence
# -> Performs regular Blast
# -> Parses a multi fasta
# -> Infile as folder and write the progress
# -> Removes the sequence if they have <50bp and of 'N'
# -> Capured as individual files.
# Author : Rajkumar (itc@rajkumar.in)
# Author : Rajkumar (itc@rajkumar.in)
# Release: JUL 2009
# Release: JUL 2009
Line 16: Line 16:
use Bio::Perl;
use Bio::Perl;
use Bio::SeqIO;
use Bio::SeqIO;
use Parallel::ForkManager;
# At the top: Time Start
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);


# Variables
# Variables
my $path           = "/home/raj/bio";
my $path       = "/home/raj/bio";
my $path1           = "$path/user/raja";
my $path1       = "$path/user/raja";
my $path2           = "$path1/blast/BES-NR";
my $path2       = "$path1/blast/BES-NR";
my $multifasta     = "$path2/BES.fa";
my $multifasta = "$path2/BES.fa";
my $infiledir      = "$path2/multifasta";
my $resdir      = "$path2/multifasta";
my $outdir          = "$path2/blastout";
my $db              = "$path/data/nr/nr";
my ($i, $total)    = $_;
my (@files, @files_)= @_;
 
 
# Multi-threading
my $midstop    = "16"; # If midstop is 0 it deals with all the files
my $cpu        = "1";
my $pm          = new Parallel::ForkManager(16);
# Getting the Grip of folders
# Getting the Grip of folders
if (-e $outdir  ) {`rm -fr $outdir`;}
if (-e $resdir)   {`rm -fr $resdir`;}
if (! -e $outdir) {`mkdir $outdir`;}
if (! -e $resdir) {`mkdir $resdir`;}


# Functions


opendir (DIR, "$infiledir") or die $!;
# making several Fasta files of GDNA
@files = readdir(DIR);
my $in  = Bio::SeqIO->new(-file => "$multifasta", -format => 'Fasta');
close(DIR);
my $identifier;
my $sequence;
my $c_files = @files - 2;
while ( my $seq = $in->next_seq() ) {
$identifier = $seq->id;
$i = 0;
  $sequence = $seq->seq;
foreach my $file(@files) {
 
next if ($file eq "."  or $file eq "..");
  #Trim the identifier
$i++;
  $identifier =~ s/_Size:[0-9]+bp//g;  
if ($i > $midstop) {last };
my ($fasta) = $_;
$fasta = ">$identifier\n$sequence\n";
$pm->start and next;
my $percent = ($i/$c_files)*100;
$percent = sprintf("%.2f", $percent);
 
print "\r\tInitiated [$i of $c_files] | $percent% ";
chdir $infiledir;
`blastall -p blastx -d $db -i $file -m 9 -o $outdir/\/$file.out -a $cpu`;
# -p program
# -d database
# -i query file
# -e evalue
# -m 8 tab delimited out
# -o outfile
# -a number of CPUs
# -q  Penalty for a nucleotide mismatch (blastn only)  default = -3
# -v  Number of database sequences to show one-line
# -K  Number of best hits from a region to keep
print "Done.";
$pm->finish;


}
my $count_seq = length($sequence);
my $count_N = ($sequence =~ tr/N//);
  my $math = $count_seq * 0.2;
 
$sequence =~ s/-/N/;
if ($count_seq < 50){ next; }
if ($count_N > $math){ next; }
open (FASTA, ">".$resdir."/".$identifier.".fa");
print FASTA $fasta;
close FASTA;
$pm->wait_all_children;
# print "$identifier\t$count_seq\n";
 
print "\nAll Done!\n\n";
 
}
# At the End: Time Calculation
my ($sec_,$min_,$hour_,$mday_,$mon_,$year_,$wday_,$yday_,$isdst_) = localtime(time);
my $days = $mday_-$mday;
my $hours = $hour_-$hour;
my $mins = $min_-$min;
my $secs = $sec_-$sec;
print "\n------------------------\nDuration: ";
if ($days>0) { print "$days days";}
if ($hours>0) { print ": $hours hours";}
if ($mins>0) { print ": $mins mins";}
if ($secs>0) { print ": $secs secs";}
print "\n------------------------\n";
 
# END
# END




</source>
</source>

Latest revision as of 11:28, 16 July 2010


#! /usr/bin/perl -w
# *************************************************************
# This Programme does the following in sequence
# -> Parses a multi fasta
# -> Removes the sequence if they have <50bp and of 'N'
# -> Capured as individual files.
# Author : Rajkumar (itc@rajkumar.in)
# Release: JUL 2009
# *************************************************************

# Libraries
	use strict;
	use Bio::Perl;
	use Bio::SeqIO;

# Variables
	my $path        = "/home/raj/bio";
	my $path1       = "$path/user/raja";
	my $path2       = "$path1/blast/BES-NR";
	my $multifasta  = "$path2/BES.fa";
	my $resdir      = "$path2/multifasta";
	
# Getting the Grip of folders
	if (-e $resdir)   {`rm -fr $resdir`;}
	if (! -e $resdir) {`mkdir $resdir`;}


# making several Fasta files of GDNA
	my $in  = Bio::SeqIO->new(-file => "$multifasta", -format => 'Fasta');
	my $identifier;
	my $sequence;
	
	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
  		
  		#Trim the identifier
  		$identifier =~ s/_Size:[0-9]+bp//g; 
		
		my ($fasta) = $_;
		$fasta = ">$identifier\n$sequence\n";

		my $count_seq = length($sequence);
		my $count_N = ($sequence =~ tr/N//);
  		my $math = $count_seq * 0.2;
  		
		$sequence =~ s/-/N/;
		if ($count_seq < 50){ next; }
		if ($count_N > $math){ next; }
				 
		open (FASTA, ">".$resdir."/".$identifier.".fa");
		print FASTA $fasta;
		close FASTA;		
	
		# print "$identifier\t$count_seq\n";
	
	}
# END
Retrieved from "https://www.rajkumar.in/index.php?title=Trim_multifasta.pl&oldid=1861"