Bioinf
view release on metacpan or search on metacpan
#______________________________________________________________
my $clu=${$_[0]} || $_[0];
$Lean_output=${$_[1]} || $_[1];
if($verbose){
print "\n# convert_clu_to_msp : \"$clu\" is given
and I am processing it with clu_to_sso_to_msp\n" if defined $clu;
}
my %clus=%{&open_clu_files(\$clu)};
my @clusters= keys %clus;
my $num_of_cluster=@clusters=@{&sort_by_cluster_size(\@clusters)};
print "# (i) $0: convert_clu_to_msp: No. of cluster=$num_of_cluster after open_clu_files \n" if $verbose;
&show_array(\@clusters) if $verbose;
&show_hash(\%clus) if $verbose;
@possible_extensions=('msp', 'msp.gz', 'msso', 'msso.gz','fsso', 'pbla', 'pbla.gz',
'ssso', 'fso', 'out', 'prot.sso', 'prot.ts');
@U_L_case=('\U', '\L');
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# Making each SINGLE linkage clu to MSP file format to be ready for divclus
#______________________________________________________________________________
for($i=0; $i< @clusters; $i++){
my (@seq_names, @final_files, $clus_name, $big_out_msp, @msp_hashes);
$clus_name=$clusters[$i];
unless($single_file_name=~/\S/){
$big_out_msp="$clus_name\_cluster\.msp"; #<<<----- final output name
}else{
$big_out_msp=$single_file_name;
}
push(@written_msp_files, $big_out_msp); ## This is the output of this sub
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If $clus_name.msp is already there, skip
#_____________________________________________
if( (-s $big_out_msp) > 100 and !$over_write ){
print "\n# (i) convert_clu_to_msp : $big_out_msp MSP file already exists, skipping\n";
print "# Use \$over_write option \'o\' to start all over again or \n";
print "# delete clustering files like XX-XX_cluster.clu to go on\n";
next ;
}
$num_of_seq_member=@seq_names=split(/ +/, $clus{$clusters[$i]}); # @seq_names has (HIU001, HI002, HI333, MJ111, etc)
print "# $0: convert_clu_to_msp: No. of seq member=$num_of_seq_member after split \n" if $verbose;
FOR0: for($j=0; $j < @seq_names; $j++){
my($sub_dir_head, $file_name_low, $file_name_up, $file_name_prot_low,
$file_name_prot_up, $file_name_low_gz, $file_name_up_gz,
$file_name_prot_low_gz, $file_name_prot_up_gz);
$each_seq_name=$seq_names[$j];
my @poss_sub_dir_heads=('.'); ## <<<<------- This is critically important, when 'D' opt is not used!
if($each_seq_name=~/(\S+)_\d+\-\d+$/){
$each_seq_name_range=$each_seq_name;
$each_seq_name=$1;
@name_types=($each_seq_name, $each_seq_name_range);
}else{
@name_types=($each_seq_name);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Here I take chars from the sequ names, as dirs have fragments of chars
#_______________________________________________________________________________
for($s=1; $s <= $subdir_char_size ; $s++){ ## here, number 2 indicates, I check single or 2 char sub dir names
$sub_dir_head= substr($seq_names[$j], 0, $s);
push(@poss_sub_dir_heads, "\L$sub_dir_head") if (-d "\L$sub_dir_head" );
push(@poss_sub_dir_heads, "\U$sub_dir_head") if (-d "\U$sub_dir_head" );
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Checking all the possible subdirectories to crop all the sso files
#_______________________________________________________________________________
FOR1: for($p=0; $p <= @poss_sub_dir_heads; $p++){ ## Default has '.' will make things like '././fam_8_8.pbla.gz'
$subd=$poss_sub_dir_heads[$p]; ## Also, the '<=' not '<' cures the same problem.
FOR2 : for($e=0; $e < @possible_extensions; $e++){
$ext=$possible_extensions[$e];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This makes all the possible lower upper case names
#______________________________________________________
for( $u=0; $u < @U_L_case; $u++){
for($v=0; $v <@name_types; $v++){
$each_seq_name=$name_types[$v];
if($U_L_case[$u]=~/U/){ $each_seq_name="\U$each_seq_name";
}else{ $each_seq_name="\L$each_seq_name"; }
if(-s "$each_seq_name\.$ext"){
push(@final_files, "$each_seq_name\.$ext" ) ;
$found_search_prog_exention_used=$ext;
$found_real_subdir_name=$subd; ## This is to report the name of the actual subd found
$found_search_prog_exention_used=$ext;
next FOR0
}elsif(-s "$each_seq_name\.$ext\.gz"){
push(@final_files, "$each_seq_name\.$ext\.gz" ) ;
$found_search_prog_exention_used=$ext;
$found_real_subdir_name=$subd; ## This is to report the name of the actual subd found
$found_search_prog_exention_used=$ext;
next FOR0
}else{
$file_wanted="\.\/$subd\/$each_seq_name\.$ext";
if(-s $file_wanted){
push( @final_files, $file_wanted);
$found_real_subdir_name=$subd; ## This is to report the name of the actual subd found
$found_search_prog_exention_used=$ext;
next FOR0
}elsif(-s "$file_wanted\.gz"){
push( @final_files, "$file_wanted\.gz");
$found_search_prog_exention_used=$ext;
$found_real_subdir_name=$subd; ## This is to report the name of the actual subd found
next FOR0;
}
}
}
}
} # FOR2
} # FOR1
print @final_files, "\n";
} # FOR0
#print "\n# @final_files \n=============> $big_out_msp \n\n";
if(@final_files < 1){
print "\n# convert_clu_to_msp :LINE no.: ", __LINE__, " ERROR: \@final_files is empty. Serious error\n";
print "\n If you have sub dir which have more than 2 chars as names, you may increase the default 2 to 3 in the above\n";
# convert_clu_to_sso_to_msp
# Options : USE, convert_clu_to_sso_to_msp, this is obsolute now
# Category :
# Version : 1.7
#--------------------------------------------------------------------------------
sub clu_to_sso_to_msp{
my($i, $j, $k, $s, $u, $p, $m, $n, $y, @possible_extensions, @list,
@final_files, @U_L_case, $file, @file, @written_msp_files);
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Opening cluster file (xx.clu)
# %clus looks like this: 2-507 YGR041W YLR353W
# 3-308 YDR222W YDR346C YLR225C
# 2-184 YCL066W YCR040W
#______________________________________________________________
my $clu=${$_[0]} || $_[0];
if($verbose){
print "\n# clu_to_sso_to_msp : \"$clu\" is given
and I am processing it with clu_to_sso_to_msp\n" if defined $clu;
}
my %clus=%{&open_clu_files(\$clu)};
my @keys= keys %clus;
my $num_of_cluster=@keys=@{&sort_by_cluster_size(\@keys)};
print "# $0: clu_to_sso_to_msp: No. of cluster=$num_of_cluster after open_clu_files \n" if $verbose;
&show_array(\@keys) if $verbose;
&show_hash(\%clus) if $verbose;
@possible_extensions=('msp', 'sso', 'msso', 'msso.gz','fsso', 'ssso', 'fso', 'out', 'prot.sso', 'prot.ts');
@U_L_case=('\U', '\L');
for($i=0; $i< @keys; $i++){
my (@list, @final_files, $clus_name, $big_out_msp, @msp_hashes);
$clus_name=$keys[$i];
unless($single_file_name=~/\S/){
$big_out_msp="$clus_name\_cluster\.msp"; #<<<----- final output name
}else{
$big_out_msp=$single_file_name;
}
push(@written_msp_files, $big_out_msp); ## This is the output of this sub
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If $clus_name.msp is already there, skip
#_____________________________________________
if( (-s $big_out_msp) > 100 and !$over_write ){
print "\n# clu_to_sso_to_msp : $big_out_msp MSP file already exists, skipping\n";
print "# Use \$over_write option \'o\' to start all over again or \n";
print "# delete clustering files like XX-XX_cluster.clu to go on\n";
next ;
}
$num_of_seq_member=@list=split(/ +/, $clus{$keys[$i]}); # @list has (HIU001, HI002, HI333, MJ111, etc)
print "# $0: clu_to_sso_to_msp: No. of seq member=$num_of_seq_member after split \n" if $verbose;
FOR0: for($j=0; $j < @list; $j++){
my($sub_dir_head, $file_name_low, $file_name_up, $file_name_prot_low, @sub_dir_heads,
$file_name_prot_up, $file_name_low_gz, $file_name_up_gz,
$file_name_prot_low_gz, $file_name_prot_up_gz);
$each_seq_name=$list[$j];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Here I take chars from the sequ names, as dirs have fragments of chars
#_______________________________________________________________________________
for($s=1; $s <=2 ; $s++){ ## here, number 2 indicates, I check single or 2 char sub dir names
$sub_dir_head= substr($list[$j], 0, $s);
push(@sub_dir_heads, "\L$sub_dir_head") if (-d "\L$sub_dir_head" );
push(@sub_dir_heads, "\U$sub_dir_head") if (-d "\U$sub_dir_head" );
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Checking all the possible subdirectories to crop all the sso files
#_______________________________________________________________________________
FOR1: for($p=0; $p < @sub_dir_heads; $p++){
$subd=$sub_dir_heads[$p];
FOR2 : for($e=0; $e < @possible_extensions; $e++){
$ext=$possible_extensions[$e];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This makes all the possible lower upper case names
#______________________________________________________
for( $u=0; $u < @U_L_case; $u++){
if($U_L_case[$u]=~/U/){ $each_seq_name="\U$each_seq_name";
}else{ $each_seq_name="\L$each_seq_name"; }
if(-s "$each_seq_name\.$ext"){ push(@final_files, "$each_seq_name\.$ext" ) ; next FOR0 }
elsif(-s "$each_seq_name\.$ext\.gz"){ push(@final_files, "$each_seq_name\.$ext\.gz" ) ; next FOR0 }
else{
$file_wanted="\.\/$subd\/$each_seq_name\.$ext";
if(-s $file_wanted){
push( @final_files, $file_wanted); next FOR0 }
elsif(-s "$file_wanted\.gz"){
push( @final_files, "$file_wanted\.gz");
next FOR0
}
}
}
} # FOR2
} # FOR1
} # FOR0
print "\n# @final_files \n=============> $big_out_msp \n\n" if $verbose;
if(@final_files < 1){
print "\n# clu_to_sso_to_msp :LINE no.: ", __LINE__, " ERROR: \@final_files is empty. Serious error\n";
print "\n If you have sub dir which have more than 2 chars as names, you may increase the default 2 to 3 in the above\n";
next;
}
# $write_each_msp_to_disk='w';
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Check if small msp files have already made in previous steps
#________________________________________________________________
if($final_files[0]=~/\.msp *$/){ ## concatenate msp into big_msp
open(BIG_MSP_FILE, ">$big_out_msp");
for($y=0; $y< @final_files; $y++){
open(SINGLE_MSP, "$final_files[$y]");
while(<SINGLE_MSP>){
print BIG_MSP_FILE $_;
}
}
close(BIG_MSP_FILE);
close(SINGLE_MSP);
push(@written_msp_files, $big_out_msp);
# clu_to_sso_to_msp
# Options :
# Category :
# Version : 1.8
#--------------------------------------------------------------------------------
sub convert_clu_to_sso_to_msp{
my($i, $j, $k, $s, $u, $p, $m, $n, $y, @possible_extensions, @list,
@final_files, @U_L_case, $file, @file, @written_msp_files);
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Opening cluster file (xx.clu)
# %clus looks like this: 2-507 YGR041W YLR353W
# 3-308 YDR222W YDR346C YLR225C
# 2-184 YCL066W YCR040W
#______________________________________________________________
my $clu=${$_[0]} || $_[0];
if($verbose){
print "\n# clu_to_sso_to_msp : \"$clu\" is given
and I am processing it with clu_to_sso_to_msp\n" if defined $clu;
}
my %clus=%{&open_clu_files(\$clu)};
my @keys= keys %clus;
my $num_of_cluster=@keys=@{&sort_by_cluster_size(\@keys)};
print "# $0: clu_to_sso_to_msp: No. of cluster=$num_of_cluster after open_clu_files \n" if $verbose;
&show_array(\@keys) if $verbose;
&show_hash(\%clus) if $verbose;
@possible_extensions=('msp', 'sso', 'msso', 'msso.gz','fsso', 'ssso', 'fso', 'out', 'prot.sso', 'prot.ts');
@U_L_case=('\U', '\L');
for($i=0; $i< @keys; $i++){
my (@list, @final_files, $clus_name, $big_out_msp, @msp_hashes);
$clus_name=$keys[$i];
unless($single_file_name=~/\S/){
$big_out_msp="$clus_name\_cluster\.msp"; #<<<----- final output name
}else{
$big_out_msp=$single_file_name;
}
push(@written_msp_files, $big_out_msp); ## This is the output of this sub
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If $clus_name.msp is already there, skip
#_____________________________________________
if( (-s $big_out_msp) > 100 and !$over_write ){
print "\n# clu_to_sso_to_msp : $big_out_msp MSP file already exists, skipping\n";
print "# Use \$over_write option \'o\' to start all over again or \n";
print "# delete clustering files like XX-XX_cluster.clu to go on\n";
next ;
}
$num_of_seq_member=@list=split(/ +/, $clus{$keys[$i]}); # @list has (HIU001, HI002, HI333, MJ111, etc)
print "# $0: clu_to_sso_to_msp: No. of seq member=$num_of_seq_member after split \n" if $verbose;
FOR0: for($j=0; $j < @list; $j++){
my($sub_dir_head, $file_name_low, $file_name_up, $file_name_prot_low, @sub_dir_heads,
$file_name_prot_up, $file_name_low_gz, $file_name_up_gz,
$file_name_prot_low_gz, $file_name_prot_up_gz);
$each_seq_name=$list[$j];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Here I take chars from the sequ names, as dirs have fragments of chars
#_______________________________________________________________________________
for($s=1; $s <=2 ; $s++){ ## here, number 2 indicates, I check single or 2 char sub dir names
$sub_dir_head= substr($list[$j], 0, $s);
push(@sub_dir_heads, "\L$sub_dir_head") if (-d "\L$sub_dir_head" );
push(@sub_dir_heads, "\U$sub_dir_head") if (-d "\U$sub_dir_head" );
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Checking all the possible subdirectories to crop all the sso files
#_______________________________________________________________________________
FOR1: for($p=0; $p < @sub_dir_heads; $p++){
$subd=$sub_dir_heads[$p];
FOR2 : for($e=0; $e < @possible_extensions; $e++){
$ext=$possible_extensions[$e];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This makes all the possible lower upper case names
#______________________________________________________
for( $u=0; $u < @U_L_case; $u++){
if($U_L_case[$u]=~/U/){ $each_seq_name="\U$each_seq_name";
}else{ $each_seq_name="\L$each_seq_name"; }
if(-s "$each_seq_name\.$ext"){ push(@final_files, "$each_seq_name\.$ext" ) ; next FOR0 }
elsif(-s "$each_seq_name\.$ext\.gz"){ push(@final_files, "$each_seq_name\.$ext\.gz" ) ; next FOR0 }
else{
$file_wanted="\.\/$subd\/$each_seq_name\.$ext";
if(-s $file_wanted){
push( @final_files, $file_wanted); next FOR0 }
elsif(-s "$file_wanted\.gz"){
push( @final_files, "$file_wanted\.gz");
next FOR0
}
}
}
} # FOR2
} # FOR1
} # FOR0
print "\n# @final_files \n=============> $big_out_msp \n\n" if $verbose;
if(@final_files < 1){
print "\n# clu_to_sso_to_msp :LINE no.: ", __LINE__, " ERROR: \@final_files is empty. Serious error\n";
print "\n If you have sub dir which have more than 2 chars as names, you may increase the default 2 to 3 in the above\n";
next;
}
# $write_each_msp_to_disk='w';
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Check if small msp files have already made in previous steps
#________________________________________________________________
if($final_files[0]=~/\.msp *$/){ ## concatenate msp into big_msp
open(BIG_MSP_FILE, ">$big_out_msp");
for($y=0; $y< @final_files; $y++){
open(SINGLE_MSP, "$final_files[$y]");
while(<SINGLE_MSP>){
print BIG_MSP_FILE $_;
}
}
close(BIG_MSP_FILE);
close(SINGLE_MSP);
push(@written_msp_files, $big_out_msp);
# Title : get_largest_file
# Usage : $largest_file_name=${&get_largest_file(@ARGV)};
# Function : checks the size of files and returns the largest
# one's name. If a file is not present in pwd or
# specified absolute path, it ignores it.
# Example :
# Keywords : choose_largest_file, largest_file, find_largest_file
# get_the_largest_file, choose_the_largest_file, get_biggest_file
# fetch_largest_file, take_largest_file, get_bigger_file, get_larger_file
# Options : _ for debugging.
# # for debugging.
# e for extract the largest from the input array
# leaving it one element less, in this case
# there will be two returning refs.
# Category :
# Version : 1.4
#--------------------------------------------------------------
sub get_largest_file{
my @in;
if(ref $_[0] eq 'ARRAY'){
@in = @{$_[0]};
}else{
@in = @_;
}
my ($largest_file, $largest, $i, $extract_opt);
for($i=0; $i< @in; $i++){
if(($in[$i]=~/^\-?e$/i)&&(!(-f $in[$i])) ){
$extract_opt=1;
splice(@in, $i, 1);
$i--;
}
}
for($i=0; $i< @in; $i++){
my $size=(-s $in[$i]);
if($size > $largest){
$largest=$size;
if($extract_opt ==1){
print "\$extract_opt is $extract_opt \n";
push(@in, $largest_file) if defined($largest_file);
$largest_file = splice(@in, $i, 1);
print "\n $largest_file \n";
$i--;
}else{
$largest_file=$in[$i];
}
}
}
if($extract_opt==1){
return(\$largest_file, \@in);
}else{ return(\$largest_file); }
}
#______________________________________________________________
# Title : get_sequence_complexity
# Usage : print "\n", ${&get_sequence_complexity(\$seq)};
# Function : caculates the single sequence's sequence complexity
# If the seq given is larger than 20, it divides it into
# frags of 20 aa and gets the average of it.
# Example : ${&get_sequence_complexity(\$seq)};
# while $seq='TTTTTACDEFGHIKLMNPQRSTVWYAAAAACCCADFADFA'
# Warning :
# Keywords : sequence_complexity, calc_sequence_complexity,
# calc_seq_complexity, get_seq_complexity, seg
# Options : _ for debugging.
# # for debugging.
# 'w=' for window size as the first arg
# Returns : Ref. for a scalar digit.
# Argument : ref. of string.
# Category :
# Version : 1.3
#--------------------------------------------------------------
sub get_sequence_complexity{
my ($complexity, @seq,$i, $j, @frag);
my $win=20;
if(ref($_[0]) eq 'ARRAY'){
@seq=@{$_[0]};
}else{
$seq=${$_[0]} || $_[0];
@seq=split(//, $seq);
}
if(defined($_[1])){ $win=${$_[1]} || $_[1]; }
if(@seq <= $win){
my (%seq, @keys);
for($i=0; $i< @seq; $i++){
$seq{$seq[$i]}++;
}
@keys= keys %seq;
$complexity=@keys/@seq;
}else{
my @frag=@{÷_array(\@seq, "s=$win")};
my @complexity=();
for($i=0; $i < @frag; $i++){
my (%seq, @keys);
my @arr=@{$frag[$i]};
for($j=0; $j< @arr; $j++){
$seq{$arr[$j]}++;
}
@keys=keys %seq;
push(@complexity, @keys/$win);
}
$complexity=${&array_average(\@complexity)};
}
return(\$complexity);
}
#______________________________________________________________________________
# Title : calc_factorial
# Usage :
# Function :
# Example :
# Keywords : calculate_factorial, get_factorial
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub calc_factorial{
my ($number, $factorial, $i);
$factorial=1;
if(ref($_[0]) eq 'SCALAR'){
$number=${$_[0]};
}else{
$number=$_[0];
}
for($i=1; $i <=$number; $i++){
$factorial=$factorial*$i;
}
return(\$factorial);
}
#______________________________________________________________
# Title : make_swiss_index
# Usage :
# Function :
# Example :
# Warning :
# Keywords :
# Options : _ for debugging.
# # for debugging.
# Returns :
# Argument :
# Category :
# Version : 1.0
#--------------------------------------------------------------
sub make_swiss_index{
my ($swiss, %index);
if(-e "$ENV{'SWDIR'}seq.dat" ){
$swiss="$ENV{'SWDIR'}seq.dat";
}elsif( -f "$ENV{'SWISS'}seq.dat" ){
$swiss="$ENV{'SWISS'}seq.dat";
}elsif( -e 'seq.dat'){
$swiss="seq.dat";
}elsif( -f "$ENV{'swiss'}seq.dat"){
# This is to fill the empty space to make complete block of seq
#________________________________________________________________
if($fill_seq_to_the_end){
$seq .="$gap_char"x($longest_seq_leng-$seq_leng); # putting '---' at the blank line end
# Above option will coredump 'seaview' seq editor
}
$seq=~s/[\-_ ]/$gap_char/g; # setting the final gap_char you like
$seq = &put_gaps_every_x_position_in_string_special($seq, 10, ' ');
sub put_gaps_every_x_position_in_string_special{
my($string); if(ref($_[0])){ $string = ${$_[0]};
}else{ $string = $_[0]; }
my($interval) = $_[1]; my($gap_char) = $_[2];
$string =~ s/(.{$interval,$interval})/$1$gap_char/g;
return($string);
}
select (MSF_FILE_OUT); ## to print out to a FILE
write MSF_FILE_OUT;
}
print "\n"; # next block starts.
}
close(MSF_FILE_OUT);
select STDOUT; # <- this is necessary to normalize output for other sub
return(\@files_created);
}
#______________________________________________________________
# Title : get_seqblock
# Usage :
# Function :
# Example : @blocks_in_hash=@{&get_seqblock(\%msf, 30)};
# Warning :
# Keywords : find_sequence_block, get_sequence_block,
# make_seq_block, make_seqblock, find_seqblock
# Options : _ for debugging.
# # for debugging.
# m= for margin length of the seqblock
# t= for threshold
# l= for min seqlet length
#
# Returns :
# Argument :
# Category :
# Version : 1.4
#--------------------------------------------------------------
sub get_seqblock{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my ($connect_gap, @seq_frag, %digitized, $verbose, %hash, $best_block_opt);
my $margin=3;
my $threshold=0.8;
my $min_seqlet_size=25;
$connect_gap=5;
my @vars=keys %vars;
for($i=0; $i< @vars; $i++){
if($vars[$i] eq 'm'){
$margin=$margin_ori=$vars{$vars[$i]};
}elsif($vars[$i] eq 't'){
$threshold=$vars{$vars[$i]};
}elsif($vars[$i] eq 'l'){
$min_seqlet_size=$vars{$vars[$i]};
}elsif($vars[$i] eq 'c'){
$connect_gap=$vars{$vars[$i]} if( defined($vars{$vars[$i]}) );
}
}
if($char_opt=~/b/){ $best_block_opt='b' }
if($char_opt=~/r/){ $range_in_name='r' }
if($char_opt=~/c/){ $connect_opt ='c' }
if($char_opt=~/v/){ $verbose='v' }
for($o=0; $o<@hash; $o++){
%hash=%{$hash[$o]};
%digitized=%{&convert_char_to_0_or_1_hash($hash[$o])};
}
%added=%{&add_columns(\%digitized)}; # 11111 + 1010101 => 2121211
&show_hash(\%added) if ($debug==1);
%blocks=%{&get_high_score_blocks(\%added,
"m=$margin", "t=$threshold", "l=$min_seqlet_size", $verbose,
"c=$connect_gap", $connect_opt, $best_block_opt, $range_in_name)};
my @keys=keys %blocks;
for($e=0; $e< @keys; $e++){
my $range="$keys[$e]\-$blocks{$keys[$e]}";
my $seq_let_leng=$blocks{$keys[$e]} - $keys[$e] + 1;
if($seq_let_leng < $min_seqlet_size){
next;
}else{
push(@RANGE, $range);
}
}
@seq_frag=&get_seq_fragments(\%hash, @RANGE,
"l=$min_seqlet_size", "$range_in_name");
return(\@seq_frag);
}
#______________________________________________________________
# Title : add_columns
# Usage :
# Function :
# Example :
# Warning : if the attached name is too long(over 12 char),
# it changes to 'Added_upX' while X is a numb.
# Keywords : add_seq_columns, add_sequence_columns,
# Options : _ for debugging.
# # for debugging.
# Returns :
# Argument :
# Category :
# Version : 1.2
#--------------------------------------------------------------
sub add_columns{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my $gap_char=' ';
for($i=0; $i<@hash; $i++){
my %hash=%{$hash[$i]};
my @names=keys %hash;
my %final_hash_out;
my @final_added;
my $out_seq_name='Add';
for($j=0; $j<@names; $j++){
$out_seq_name.= "_$names[$j]";
my $string =$hash{$names[$j]};
my @ar_string;
if($string=~/\d{1,5}[ \,]\d{1,5}[ \,]\d{1,5}/){
@ar_string =split(/$gap_char|\,/, $string );
}elsif($string=~/^\d{5,}$/){
@ar_string =split(//, $string );
}
for($s=0; $s < @ar_string; $s++){
$final_added[$s]=$ar_string[$s]+$final_added[$s];
}
}
if(length($out_seq_name) > 12){ $out_seq_name="Added_up${i}"; }
$final_hash_out{$out_seq_name}=join("$gap_char", @final_added);
push(@OUT_HASH, \%final_hash_out);
}
wantarray ? return(@OUT_HASH) : return($OUT_HASH[0]);
}
#____________________________________________________________________
# Title : get_high_score_blocks
# Usage : get_high_score_blocks(<ref. of hash for number string>)
print "\n\|";
print "\n\|";
for($m=0; $m< $largest; $m++){
print "\n\|";
for($i=0; $i<@inversed;$i++){
if($inversed[$i] > 0){
print " ";
$inversed[$i]--;
}else{
print "\*";
}
}
}
print "\n\|", "\-"x@numbers;
print "\n";
}
#______________________________________________________________
# Title : condense_number_string
# Usage :
# Function : condenses the numbers by making an average with
# given factor. If the factor is 2 on number seq
# 1334284425 , result will be 23543
# 133428442 , 23541 <-- preserved end
# Factor 3 =>
# 133428442 , (1+3+3)/3 = 2
# (4+2+8)/3 = 4,,,
# Example : @output=@{&condense_number_string(\@input, $factor)};
# with @input=qw(1 2 4 10 10 22 2 3 44 2 3); and $factor=3
# Warning :
# Keywords : compact_number_string, compact_digits, condense
# condense_string
# Options : _ for debugging.
# # for debugging.
# Returns :
# Argument :
# Category :
# Version : 1.1
#--------------------------------------------------------------
sub condense_number_string{
my @ARRAY=@{$_[0]};
my $factor = ${$_[1]} || $_[1];
my ( $i, $j, @out );
unless(defined($factor)){ $factor=1 }
for($i=0; $i< @ARRAY; $i+=$factor){
my $temp_sum;
for($j=$i; $j < ($factor+$i); $j++){
$temp_sum+=$ARRAY[$j]
}
push(@out, int($temp_sum/$factor) );
}
return(\@out);
}
#___________________________________________________________
# Title : get_seq_fragments
# Usage : @seq_frag=&get_seq_fragments(\%msf, @RANGE);
# Function : gets sequence(string) segments with defined
# ranges.
# Example :
# %test=('seq1', '1234AAAAAAAAAAAaaaaa', 'seq2', '1234BBBBBBB');
# @range = ('1-4', '5-8');
#
# %out = %{&get_seq_fragments(\%test, \@range)};
# %out => (seq1_5-8 AAAAA
# seq2_5-8 BBBBB
# seq1_1-4 1234
# seq2_1-4 1234 )
#
# Warning :
# Keywords : get_sequence_fragments,
# Options : _ for debugging.
# # for debugging.
# l= for min seqlet length
# r for adding ranges in the seq names
#
# Returns :
# Argument :
# Category :
# Version : 1.8
#-------------------------------------------------------
sub get_seq_fragments{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my $min_seqlet_size=10;
my @vars=keys %vars;
my $no_range_in_name=1;
for($i=0; $i< @vars; $i++){
if($vars[$i] eq 'l'){
$min_seqlet_size=$vars{$vars[$i]};
}
}
if($char_opt=~/v/){ print "\n \$char_opt is $char_opt @char_opt\n"; }
if($char_opt=~/n/){ $no_range_in_name = 1 }
if($char_opt=~/r/){ $no_range_in_name = 0 }
print "\nget_seq_fragments \$no_range_in_name is $no_range_in_name \n";
for($i=0; $i< @hash; $i++){
my (%out_frag, $frag_name, $range_start, $range_end, @out_hash);
my %seqs = %{$hash[$i]};
my @names = keys %seqs;
if(@names==1){
for($j=0; $j < @names; $j++){
my $seq_name = $names[$j];
my $seq = $seqs{$seq_name};
for($k=0; $k< @range; $k++){
my $range = $range[$k];
if($no_range_in_name==1){
$frag_name = "$seq_name";
}else{
$frag_name = "$seq_name\_$range";
}
#if(length($frag_name)>14 ){
# $frag_name ='x'."${j}_${range}";
#}
($range_start, $range_end)=$range=~/(\d+\.?\d*)\-(\d+\.?\d*)/;
my $frag_len = $range_end-$range_start+1;
if($frag_len < $min_seqlet_size){
next;
}
my $fragment = substr($seq, $range_start-1, $frag_len);
$out_frag{$frag_name}=$fragment;
}
}
push(@out_hash, \%out_frag);
}elsif(@names > 1){
for($k=0; $k< @range; $k++){
my %out_frag=();
my $range=$range[$k];
($range_start, $range_end)=$range=~/(\d+\.?\d*)\-(\d+\.?\d*)/;
my $frag_len = $range_end-$range_start+1;
if($frag_len < $min_seqlet_size){
next;
}
for($j=0; $j < @names; $j++){
my $seq_name=$names[$j];
my $seq = $seqs{$seq_name};
if($no_range_in_name==1){
$frag_name = "$seq_name";
}else{
$frag_name = "$seq_name\_$range";
}
#if(length($frag_name)>15 ){
# $frag_name ='x'."${j}_${range}";
#}
if($range_start==0){ $range_start++; } ## This is a bugfix
my $fragment = substr($seq, $range_start-1, $frag_len);
$out_frag{$frag_name}=$fragment;
}
push(@out_hash, \%out_frag);
}
}
}
if(@out_hash > 1){ return(@out_hash)
}elsif(@out_hash==1){ return($out_hash[0]) }
}
#________________________________________________________________________
# Title : make_standalone_subroutines
# Usage : &make_standalone_subroutines(@ARGV);
# Example : &make_standalone_subroutines(@ARGV);
# Function : Creates each subroutine derived xxx.pl file from Bio.pl or any
# given library file. If there is a file for a sub already, it
# skips.
# Class : Utility
# Keywords :
# Options :
# Author : jong@salt2.med.harvard.edu
# Category :
# Version : 1.1
#--------------------------------------------------------------------
sub make_standalone_subroutines{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my($each_sub, %out_subs, %left_out, @lib, $ver, $real_sub_entry_found);
$|=1;
for($i=0; $i < @file; $i++){
open(LIB_FILE, "<$file[$i]")|| die "\n $file[$i] <- $! \n";
my @lib =<LIB_FILE>;
FOR: for($j=0; $j < @lib; $j++){
my (%out_subs, $each_sub);
my $title_found;
#"""" Taking the headbox """""""""""""
if( ($lib[$j]=~/^#+[_\-\*]{10,120} *$/)
&&($lib[$j+1]=~/^(# *title *: *([\w\-]+))[^\.pl]/i) ){
$each_sub=$2;
$title_found =1;
if( (-s "$each_sub\.pl") > 200 ){
print (-s "$each_sub\.pl"), " ";
print " $each_sub", " exists \n";
next FOR;
}elsif((-s "$each_sub\.pl") <= 200){
my $temp;
open (TEMP, "<$each_sub\.pl");
while(<TEMP>){
if(/^#[_\-\*]{10,120} *$/){ $temp++ }
elsif(/^# *title *: *[\w\-]+[^\.pl]/i ){
$temp++;
}elsif(/^# *\w+/){
# @out=@{&push_if_not_already(\@mother_array, $adding_scalar)};
# Function : returns ref. of an array for a list of non-repetitive entry.
# Example :
# Warning :
# Keywords : add_if_not_already, add_element_if_not_already, if_not_already
# add_element_if_not_already, push_element_if_not_already,
# if_no_already_push, put_element_if_not_already, add_new_element
# add_new_items_only, push_new_items_only, push_new_elements_only
# put_if_not_already,
# Options :
# Returns : a ref. of an array.
# Argument : two references. The first should be an array ref. The 2nd can be either
# scalar or array reference.
# Category :
# Version : 1.3
#--------------------------------------------------------------------
sub push_if_not_already{
my($already_in, $already, $i, @push_items_given);
my(@out_array)=@{$_[0]};
if(ref($_[0]) ne 'ARRAY'){ print "\n push_if_not_array need ref\n"; exit; }
push(@push_items_given, ${$_[1]}) if(ref($_[1]) eq 'SCALAR');
@push_items_given=@{$_[1]} if(ref($_[1]) eq 'ARRAY');
for $already (@out_array){ ## This for is to remove repetitive
for ($i=0; $i< @push_items_given; $i++){
if($already eq $push_items_given[$i]){ splice(@push_items_given,$i, 1); $i--; }
}
}
push(@out_array, @push_items_given);
return(\@out_array);
}
#_______________________________________________________________
# Title : compare_sec_template_with_db
# Usage :
# Function :
# Example :
# Warning :
# Keywords : sec structure mapping, map sec str, map_sec_structure
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.0
#-----------------------------------------------------------
sub compare_sec_template_with_db{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
$ref_target_hash = shift @hash;
%target = %{$ref_target_hash};
@names = keys %target;
$name_target = $names[0];
if($name_target =~/\w+(\d+)/){ $ori_target_seq_len = length($1); }
@target_frag = split(/ +/, $target{$name_target} );
for($i =0; $i< @hash; $i ++){
%db=%{$hash[$i]};
@db_name = keys %db;
for($j=0; $j < @db_name; $j ++){
$name = $db_name[$j];
@db_frag = split( / +/, $db{$name} );
for($k=0; $k < @target_frag; $k ++){
if( ($target_frag[$k]=~/H(\d+)/i)&&($db_frag[$k]=~/H(\d+)/i) ){
$simple_match_output{$name}++;
$leng_diff = abs($1 - $2)/15;
$simple_match_output{$name} = $simple_match_output{$name}- $leng_diff;
}elsif( ($target_frag[$k]=~/E(\d+)/i)&&($db_frag[$k]=~/H(\d+)/i) ){
$simple_match_output{$name}--;
#$leng_diff = abs($1 - $2)/10;
#$simple_match_output{$name} = $simple_match_output{$name}- $leng_diff;
}elsif( ($target_frag[$k]=~/H(\d+)/i)&&($db_frag[$k]=~/E(\d+)/i) ){
$simple_match_output{$name}--;
#$leng_diff = abs($1 - $2)/10;
#$simple_match_output{$name} = $simple_match_output{$name}- $leng_diff;
}elsif( ($target_frag[$k]=~/E(\d+)/i)&&($db_frag[$k]=~/E(\d+)/i) ){
$simple_match_output{$name}++;
$leng_diff = abs($1 - $2)/15;
$simple_match_output{$name} = $simple_match_output{$name}- $leng_diff;
}
}
}
}
return(\%simple_match_output);
}
#___________________________________________________________________
# Title : get_peptide_occurance
# Usage :
# Function : gets the number of occurances of peptide(with given size) for
# any number of sequences given.
# Example : %stat=%{&get_peptide_occurance(\%pro_sequence, $size)};
# while %pro_sequence has one or more sequences like
# seq1 AAAAAAAAAAAA, seq2 BBBBBBBBBBBBBB, ...
# $size is number. For dipeptide=2, tripeptide=3, tetrapep=4...
# Warning :
# Keywords :
# Options :
# Returns :
# Argument : eg=> (\%ref_hash, 4)
# Category :
# Version : 1.2
#---------------------------------------------------------------
sub get_peptide_occurance{
my($k, $i, $s, $peptide, $pep_entry_num, @name, %stat);
%all=%{$_[0]};
$size=$_[1];
@name=keys %all;
for($k=0; $k<@name; $k++){
if($all{$name[$k]}=~/[BZX]/i){ next;
}else{
$seq_leng += length($all{$name[$k]});
my @seq=split(//, $all{$name[$k]});
my $all_occur_pep;
for($i=0; $i< (@seq-($size-1)); $i ++){
my $peptide;
for($s=0; $s < $size; $s++){
$peptide .= $seq[$i+$s];
}
$stat{$peptide}++;
$all_occur_pep ++;
print "\n$peptide $stat{$peptide}" if $debug==1;
$pep_entry_num=keys %stat;
if( ($debug==1)&&($pep_entry_num%100 == 0 ) ){
print "\n Present peptide entries are: $pep_entry_num out of $all_occur_pep residues \n";
}
}
}
}
return(\%stat);
}
#___________________________________________________________________
push(@bonus, $7);
}
}
if($single_array_return == 1){ @allbins=( @bin1, @bin2, @bin3, @bin4, @bin5, @bin6, @bonus); }
else{ return(\@bin1, \@bin2, \@bin3, \@bin4, \@bin5, \@bin6, \@bonus); }
}
#_________________________________________________________________
# Title : get_probable_half
# Usage :
# Function : This produces a hash ref. which is supposed to be most probable
# according to the given array. It divides array into halves
# gets the more probable half until it gets one single number.
# Example :
# Warning :
# Keywords : get_frequent_halves,
# Options :
# Returns :
# Argument : \@array
# Category :
# Version : 1.0
#-----------------------------------------------------------------
sub get_probable_half{
my (%hash, $i, @keys, @values);
%hash=%{$_[0]};
@keys= sort {$a<=$b} keys %hash;
@values= values %hash;
print "\n Hash is ",%hash,"\n";
if(@keys == 1){
return(\%hash); last;
}elsif(@keys >1){
if((@keys % 2) != 0){ # make the number even
$keys[@keys]=$keys[$#keys];
}
@first_half=@keys[0..(@keys/2-1)];
@second_half=@keys[(@keys/2)..$#keys];
print "@first_half", "\n";
print "@second_half", "\n";
my($sum1, $sum2, %hash1, %hash2);
for($i=0; $i<(@keys/2); $i++){
$sum1 +=$hash{$first_half[$i]};
$sum2 +=$hash{$second_half[$i]};
$hash1{$first_half[$i]} =$hash{$first_half[$i]};
$hash2{$second_half[$i]}=$hash{$second_half[$i]};
}
if($sum1 > $sum2){ &get_probable_half(\%hash1);}
else{ &get_probable_half(\%hash2);}
}
}
#_______________________________________________________________
# Title : divide_array
# Usage : &show_array(÷_array(\@input, 6));
# Function : divides any array to the denominator given.
# If you give array of 100 elem, with 5, you will
# get 5 arrays with 20 elem each.
# Example :
# Warning :
# Keywords : split_array_into_pieces, split_array, chop_array,
# fragment_array,
# Options : s= for dividing the array with sub array size
# eg) to get 20 elem length sub arrays from
# a big array
# @ar_ref=@{÷_array(\@array, 's=20')};
# Returns :
# Argument :
# Category :
# Version : 1.4
#-----------------------------------------------------------
sub divide_array{
my ($size,$remaining, $size_div, $s);
my @array = @{$_[0]};
my @final_array_ref=();
if(ref($_[1])){
if(${$_[1]}=~/^\d+$/){
$denominator = ${$_[1]};
}elsif(${$_[1]}=~/s=(\d+)$/i){
$size=$1;
$size_div=1;
}
}elsif($_[1]=~/^\d+$/){
$denominator = $_[1];
}elsif($_[1]=~/s=(\d+)$/i){
$size=$1;
$size_div=1;
}
if((@_ ==1)&&($denominator == 0)){
print "\n Denominator is 0, error, setting to 1\n";
$denominator = 1;
}
if($size_div==1){
while(@array){
push(@final_array_ref, [splice(@array, 0, $size)]);
}
}else{
my $frag_ar_size = int(@array/$denominator);
if($debug eq 1){ print "\n Frag arr size is : $frag_ar_size \n" }
$remaining = @array % $denominator;
if($debug eq 1){ print "\n Remnant elem size is : $remaining \n" }
for($i=0; $i < $denominator; $i++){
if($remaining > 0){
push(@final_array_ref, [splice(@array, 0, ($frag_ar_size+1),)] );
$remaining --;
}elsif(($remaining == 0)&&(@array>0)){
push(@final_array_ref, [splice(@array, 0, ($frag_ar_size),)] );
}
}
}
return(\@final_array_ref);
}
#__________________________________________________________________________
# Title : split_fasta_files
# Usage : @names_of_single_files=@{&split_fasta_files(\@files)};
# Function :
# Example :
# Keywords : divide_fasta_files, split_fasta_db_files, divide_fasta_db_files
# make_single_fasta_files, write_single_fasta, write_single_fasta_files
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.0
#----------------------------------------------------------------------------
sub split_fasta_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my (@each_single_seq_names, $out_single_fasta_file_name );
for($i=0; $i< @file; $i++){
open(FASTA, "$file[$i]");
while(<FASTA>){
if(/\> *(\S+)/){
$out_single_fasta_file_name="$1\.fa";
open(OUT_SINGLE, ">$out_single_fasta_file_name");
print OUT_SINGLE $_;
push(@each_single_seq_names, $1);
}elsif(/\w+/){
print OUT_SINGLE $_;
}
}
close OUT_SINGLE;
}
return(\@each_single_seq_names);
}
#______________________________________________________________________________
# Title : split_files
# Usage :
# Function :
# Example :
# Keywords :
# Options :
# $division_factor= by d=
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub split_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my($division_factor, $base_name);
$division_factor=2;
if($vars{'d'}=~/(\d+)/){ $division_factor=$1 };
for($i=0; $i< @file; $i++){
open(FILE, $file[$i]);
$base_name=${&get_base_names($file[$i])};
$extension=${&get_extension_names($file[$i])};
$total_line_num=@lines=<FILE>;
$splited_frag_size=int(@lines/$division_factor);
for($j=0; $j< $division_factor; $j++){
#$frag_file="$base_name\_s${j}\.split${j}";
$frag_file="$base_name\_s${j}\.$extension";
push(@splited_files, $frag_file);
open(FRAGMENT_FILE, ">$frag_file");
for($k=0; $k<= $splited_frag_size; $k++){
print FRAGMENT_FILE $lines[$k];
}
splice(@lines, 0, $splited_frag_size);
}
print FRAGMENT_FILE @lines; # Writing the very last remaining lines
close (FRAGMENT_FILE);
}
return(\@splited_files);
}
#_______________________________________________________________
# Title : split_sequence
# Usage : %out=%{&split_sequence(\%input, 2 )};
# Function : divides any string to the denominator given.
# Example : &show_array( ÷_string(\%input, 3) );
# while $input is 'seq', '12345789ABCDEFHIJKLMN'
# The output will be 'seq_1_half', '1234578'
# 'seq_2_half', '9ABCDEF'
# 'seq_3_half', 'HIJKLMN'
# Warning :
# Keywords : divide_string, split_string, chop_string, divide_sequence
# split_sequence(look at separate split_sequence sub),
# Options :
# $reverse_second_half=S by S -S
# $reverse_first_half =F by F -F
# $reverse_rest =R by R -R ## reversing all except the first
# $reverse_all =A by A -A # reverse all the fragments
# Returns :
# Argument :
# Category :
# Version : 1.3
#-----------------------------------------------------------
sub split_sequence{
my(@string, $frag_str_size, $remaining, @hash, $i, $j,
@seq_names, $denominator, %input_seq, @final_hash_ref,
$frag, $reverse_second_half, $reverse_first_half,
$reverse_rest, $reverse_all);
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Checking arguments
#________________________________________________
for($i=0; $i< @_; $i++){
if(ref($_[$i]) eq 'HASH'){
push(@hash, $_[$i]);
}elsif(ref($_[$i]) eq 'SCALAR'){
if(${$_[$i]}=~/^(\d+\.?\d*)$/){ # this can handle fraction number!!
$denominator = int($1); # like 9.5
}
}elsif($_[$i]=~/^(\d+\.?\d*)$/){ # this can handle fraction number!!
$denominator = int($1); # like 9.5
}elsif($_[$i]=~/F/){ # this can handle fraction number!!
$reverse_first_half='F';
}elsif($_[$i]=~/S/){ # this can handle fraction number!!
$reverse_second_half='S';
}elsif($_[$i]=~/R/){ # this can handle fraction number!!
$reverse_rest='R'; $reverse_second_half='s';
}elsif($_[$i]=~/A/){ # this can handle fraction number!!
$reverse_all='A'; $reverse_rest='r'; $reverse_first_half='f';
$reverse_second_half='S';
}elsif($_[$i]=~/\S/){
print "\n# $0: split_sequence, You put some strange stuff to me!!\n\n";
exit;
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~``
# When s option is given while denominator is over 2, set r option
#____________________________________________________________________
if($denominator > 2 and $reverse_second_half){
$reverse_rest='r';
}
if($denominator == 0){
print "\n Denominator is 0, error, setting to 1\n";
$denominator = 1;
}
print "\n# split_sequence: \$denominator is $denominator, with \@hash\n";
for($i=0; $i< @hash; $i++){
my %input_seq=%{$hash[$i]};
@seq_names=keys %input_seq;
my (%out_seq);
for($k=0; $k < @seq_names; $k++){
$seq_name=$seq_names[$k];
@string = split(//, $input_seq{ $seq_name });
$frag_str_size = int(@string/$denominator);
if($debug eq 1){ print "\n Frag str size is : $frag_str_size \n" }
$remaining = @string % $denominator;
for($j=0; $j < $denominator; $j++){
$frag_number=$j+1;
if($remaining > 0){
$frag=join('', splice(@string, 0, ($frag_str_size+1) ) );
$remaining --;
$seq_name_split="$seq_name\_$frag_number";
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~```
# Before final, set the seq types
#_________________________________________________________________
if($frag_number==1 and $reverse_first_half){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($frag_number==2 and $reverse_second_half){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($frag_number >2 and $reverse_rest){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($reverse_all){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}
$out_seq{$seq_name_split}=$frag;
}elsif( $remaining == 0 and @string>0){
$seq_name_split="$seq_name\_$frag_number";
$frag=join('', splice(@string, 0, $frag_str_size,) );
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~```
# Before final, set the seq types
#_________________________________________________________________
if($frag_number==1 and $reverse_first_half){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($frag_number==2 and $reverse_second_half){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($frag_number >2 and $reverse_rest){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}elsif($reverse_all){
$frag=reverse($frag);
$seq_name_split="$seq_name_split\_rv";
}
$out_seq{$seq_name_split}=$frag;
}
}
}
push(@final_hash_ref, \%out_seq);
}
wantarray? \@final_hash_ref : $final_hash_ref[0];
}
#_______________________________________________________________
# Title : divide_string
# Usage : &show_array(÷_string(\$input, 6));
# Function : divides any string to the denominator given.
# Example : &show_array( ÷_string(\$input, 3) );
# while $input is '12345789ABCDEFHIJKLMN'
# The output will be '1234578 9ABCDEF HIJKLMN'
# Warning :
# Keywords : divide_string, split_string, chop_string, divide_sequence
# split_sequence(look at separate split_sequence sub),
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.4
#-----------------------------------------------------------
sub divide_string{
my (@array, $i, $j, $denominator, @temp, @string, $frag_str_size,
$remaining, $frag);
for($i=0; $i< @_; $i++){
if(ref($_[$i]) eq 'ARRAY'){
push(@array, @{$_[$i]});
}elsif(ref($_[$i]) eq 'SCALAR'){
if(${$_[$i]} =~/^(\d+)$/){
$denominator = $1;
}else{
push(@array, ${$_[$i]});
}
}elsif($_[$i]=~/^(\d+\.?\d*)$/){ # this can handle fraction number!!
$denominator = int($1); # like 9.5
}else{
push(@array, $_[$i]);
}
}
if($denominator == 0){
print "\n Denominator is 0, error, setting to 1\n";
$denominator = 1;
}
for($i=0; $i< @array; $i++){
my @temp;
my @string = split(//, $array[$i]);
$frag_str_size = int(@string/$denominator);
if($debug eq 1){ print "\n Frag str size is : $frag_str_size \n" }
$remaining = @string % $denominator;
for($j=0; $j < $denominator; $j++){
if($remaining > 0){
$frag=join('', splice(@string, 0, ($frag_str_size+1) ) );
push(@temp, $frag);
$remaining --;
}elsif( $remaining == 0 and @string>0){
$frag=join('', splice(@string, 0, $frag_str_size,) );
push(@temp, $frag);
}
}
push(@final_array_ref, \@temp);
}
wantarray? \@final_array_ref : $final_array_ref[0];
}
#____________________________________________________________
# Title : write_html_headbox
# Usage : &write_html_headbox($outfilename, \%entries);
# Function : write html format headbox explanation with
# given hashes of headbox content.
# Example :
# Warning : It takes off the last '/' when $URL has it
#
# Keywords : write_headbox_html, write headbox in html,
# write_headbox_files
# Options : 'd' for date inclusion at the top of the page
# f= for default ftp dir name
#
# Category :
# Version : 1.7
#-----------------------------------------------------------
sub write_html_headbox{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my(%in, $name, $name1);
my $URL='ftp://cyrah.med.harvard.edu/pub/Perl/';
if($vars{'f'}=~/(\S+)\/$/){ ## checking if 'f' opt is set(for default ftp dir name)
$URL=$1;
}
for ($f=0; $f<@file; $f++){
$output_file=$file[$f];
open(FILE, ">$output_file");
print FILE "\<html\>\n";
if($char_opt =~/d/i){
print FILE "\<H2 ALIGN=CENTER\>";
print FILE "\<FONT COLOR=\"#DC143C\"\>$file[$f]\<\/H2\>\<\/FONT\>\n";
print FILE "\<br\>\n";
print FILE "\<H6 ALIGN=RIGHT\>", `date`, "Created by \<A href=\"$0\.html\"\>$0\<\/A\>\<\/H6\> \<hr\>";
}
for($i=0; $i< @hash; $i++){
my %in = %{$hash[$i]};
my @keys = sort keys %in;
$name= $in{'Title'};
if($name=~/(\S+)\.pl/){ ## to prevent making 'xxxx.pl.pl'
$name1=$1;
}else{
$name1=$name;
.
format STDOUT =
>P1;@<<<<<<<<<<<<<<<
$name
structureX: @<<<<<<<<<<<<<<: @<<<: : @<<<<<: : @<<<<<<<<<<<<<: :
$name $start_seq $seq_leng $name
@*
$seq
*
.
########## Program starts ####################
if($file[0]){
open (MODELLER_ALI_FORMAT, ">$file[0]");
}else{
@seq = %{$hash[0]};
$name = $seq[0]; print "\n\$name is $name" if $debug ==1;
open (MODELLER_ALI_FORMAT, ">$name.ali");
}
for($i=0; $i<@hash;$i++){
@seq = %{$hash[$i]};
$name = $seq[0]; print "\n\$name is $name" if $debug ==1;
$seq = $seq[1]; print "\n\$seq is $seq" if $debug ==1;
$seq_leng = length($seq);
## checking PDB entry of the template ##
$pdb1 = "$ENV{'PDB'}\/$name.brk";
$pdb2 = "$ENV{'PDB'}\/$name.pdb";
if( (-e $pdb1 ) || ( -e $pdb2 ) ){
$pdb_file = $pdb1;
$start_seq = ${&get_pdb_file_start_number( $pdb_file )};
}
## Handling the starting and ending seq points of template
if( defined($seq[2])&&($seq[2]=~/^\d+/) ){
$start_seq=$seq[2]; $seq_leng+=$start_seq; }
elsif( defined($start_seq) ){
$seq_leng+=$start_seq; }
else{ $start_seq =1; }
if( (defined($seq[2])) &&($end_seq != $seq_leng ) ){
print "\n Your template seq length does not match with actual seq size
\n I will put the calculated value \"$seq_leng\" as the template length\n\n";
}
print "\n\$seq_leng is $seq_leng\n" if $debug ==1;
write MODELLER_ALI_FORMAT;
if( $char_opt=~/v/i){ write STDOUT; }
}
}
#_______________________________________________________________
# Title : make_template_from_sec_str
# Usage : %target = %{&make_template_from_sec_str(\%seq)};
# Function : makes template of sec. str. like: 'H5 E4 E2' out of '__HHHHH__EEEE__EE__'
# Example :
# Warning :
# Keywords :
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.1
#-----------------------------------------------------------
sub make_template_from_sec_str{
my(%out, @name, @fragments, %in, $name, $leng, $frag_seq, $name2);
#""""""""""""""""""""""< handle_arguments{ head Ver 1.1 >""""""""""""""""""""""""""""""""
my(@A ) = &handle_arguments( @_ ); my( $num_opt )=${$A[7]};my( $char_opt )=${$A[8]};
my(@hash) =@{$A[0]};my(@file) =@{$A[4]};my(@dir ) =@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};my(@raw_string)=@{$A[9]};
my($i, $j, $c, $d, $e, $f, $g, $h, $k, $l, $p, $q, $r, $s, $t, $u, $v, $w, $x,$y,$z);
#"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
for($j = 0; $j < @hash ; $j ++){
%in = %{$hash[$j]};
@name = keys %in;
print "\@name is @name\n", if($debug eq 1);
for($t=0; $t < @name; $t++){
$name = $name[$t];
print "\$name is $name\n", if($debug eq 1);
$leng=length($in{$name});
print "$leng\n", if($debug eq 1);
$name2 ="$name"."$leng"; # to attach sequence length
@fragments = split(/_+/, $in{$name});
print "\@fragments is @fragments\n", if($debug eq 1);
for($i = 0; $i < @fragments; $i++){
if($fragments[$i] =~/(\w)\w+/){
$fraglength = length($fragments[$i]);
$frag_seq .= "$1"."$fraglength "; # space is delimiter 'H5 E3 E5 E4'
print "\$frag_seq is $frag_seq\n", if($debug eq 1);
}
$out{$name2}=$frag_seq;
}
}
}
return(\%out);
}
#_______________________________________________________________
# Title : calculate_protein_volume
# Usage : %volumes=%{&calculate_protein_volume(\%seq)}
# Function :
# Example :
# Warning :
# Keywords :
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.0
#-----------------------------------------------------------
sub calculate_protein_volume{
my %final_volume;
my %volume=("A","88.6" ,"C","108.5","D","111.1","E","138.4",
"F","189.9","G","60.1" ,"H","153.2","I","166.7",
"K","168.6","L","166.7","M","162.9","N","117.7",
"P","122.7","Q","143.9","R","173.4","S","89.0",
"T","116.1","V","140.0","W","227.8","Y","193.6",
"a","88.6" ,"c","108.5","d","111.1","e","138.4",
"f","189.9","g","60.1" ,"h","153.2","i","166.7",
"k","168.6","l","166.7","m","162.9","n","117.7",
"p","122.7","q","143.9","r","173.4","s","89.0",
"t","116.1","v","140.0","w","227.8","y","193.6");
#""""""""""""""""""""""< handle_arguments{ head Ver 1.1 >""""""""""""""""""""""""""""""""""""""
my(@A ) = &handle_arguments( @_ ); my( $num_opt )=${$A[7]};my( $char_opt )=${$A[8]};
my(@hash) =@{$A[0]};my(@file) =@{$A[4]};my(@dir ) =@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};my(@raw_string)=@{$A[9]};
my($i, $j, $c, $d, $e, $f, $g, $h, $k, $l, $p, $q, $r, $s, $t, $u, $v, $w, $x,$y,$z);
#"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my %seq_hash = %{$hash[0]};
my @names = keys %seq_hash;
#if($debug == 1){
print "\n ",__LINE__, " Seq names are @names\n";
# }
for( $i=0 ; $i < @names; $i++){
my @seq=split(//, $seq_hash{$names[$i]});
for( $j=0; $j < @seq; $j ++){
$final_volume{$names[$i]} += $volume{$seq[$j]};
}
}
\%final_volume;
}
#_______________________________________________________________
# Title : extract_words
# Usage : @words = @{&extract_words(\$string)};
$scop_class=$2;
if($char_opt=~/R/ and $range !~/\d+\-\d+/){ $range='' }else{ $range="_$range" }
splice(@string, $j, 1); ## for optimization
$sequence{"${pdb_name}${range}"} = $scop_class;
return(\%sequence) if(@string < 1);
}
}
}
}else{
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Simplified style reading
#_______________________________________________
while(<FILE_1>){
if(/^ *\> *($string[$j]) +([\d\.]+) *\(?(\d*\-?\d*)\)?/){
$pdb_name=$1;
$range=$3;
$scop_class=$2;
if($NO_range){ $range='' }elsif($range=~/\d+\-\d+/){ $range="_$range" }
$sequence{"${pdb_name}${range}"} = $scop_class;
}
}
}
}else{ # getting all seq in the given file(s)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Simplified style reading
#_______________________________________________
if($pdb_style_name){
while(<FILE_1>){
if(/^ *\> *[cd](\S\S\S\S)(\S)(\S) +([\d\.]+) *\(?(\d*\-?\d*)\)?/){
$pdb_name=$1;
$range=$5;
$domain=$2; #<-- like A or B or C
$domain_num=$3;
$scop_class=$4;
if($range !~/\d+\-\d+/){ $range='' }else{ $range="_$range" }
if($domain !~ /^_$/){
$domain="\U$domain";
}else{ undef($domain) }
if($domain_num !~/^_$/){
$domain_num=$3;
}else{ undef($domain_num) }
$sequence{"${pdb_name}${domain}${range}"} = $scop_class;
}else{ next }
}
}else{ #### This takes only the real names rather than dxxxx__ style
while(<FILE_1>){
if(/^ *\> *(\S+) +([\d\.]+)/){
$pdb_name=$1;
$scop_class=$2;
$sequence{"${pdb_name}"} = $scop_class;
}
}
}
}
close FILE_1;
}
if(defined(@range)){
%seq_fragments=%{&get_seq_fragments(\%sequence, \@range)};
return(\%seq_fragments);
}elsif($char_opt=~/s/i){ # when SIZE return only option is set
@seq_Names=keys %sequence;
for($i=0; $i < @seq_Names; $i++){
$Sizes{$seq_Names[$i]}=length($sequence{$seq_Names[$i]});
}
return(\%Sizes);
}elsif(@hash >=1){
for($h=0; $h< @hash; $h++){
my %hash=%{$hash[$h]};
my @Keys=keys %hash;
for($k=0; $k<@Keys; $k++){
if(defined($hash{$Keys[$k]})){
($S_start, $S_end)=($hash{$Keys[$k]}=~/(\d+)\-(\d+)/);
$sequence{$Keys[$k]}=substr($sequence{$Keys[$k]}, ($S_start-1), ($S_end-$S_start));
}
}
}
return(\%sequence);
}else{
return(\%sequence);
}
}
#________________________________________________________________________
# Title : open_cel_files
# Usage :
# Function :
# Example :
# Example INPUT file:
# # Tabl delimited data file
# X Y MeanIn STDEV Pixel num
#
# 0 0 200.0 15 16 1.24343
# 1 0 200.0 15 16 1.24343
# 2 0 200.0 15 16 1.24343
# 3 0 200.0 15 16 1.24343
#
# Warning :
# Keywords :
# Options :
# Returns : 3D array =( [X][Y][0], [X][Y][1], [X][Y][2], [X][Y][3])
# Argument :
# Version : 1.1
#--------------------------------------------------------------------
sub open_cel_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my($X_axis, $Y_axis, $intensity, $stdev_of_intensity, $pixel_number,
$temp, @out_3D_array, @out_3D_array_ref_list);
# If the seq name has ranges like XXXXXX_1-30, it will only
# return 1-30 of XXXXXX sequence. EMBL file is from SANGER
# Centre for genome seq.
#
# FASTA sequence file format is like this;
#
# FT /db_xref="PID:e322131"
# FT /db_xref="SPTREMBL:O00531"
# FT /product="tenascin-R (restrictin)"
# FT /translation="CLFLHPLIDLVPGTEYGVGISAVMNSQQSVPATMNARTELDSPRD
# FT IISVTAERGRQQSLESTVDAFTGFRPISHLHFSHVTSSSVNITWSDPSPPADRLILNYS
# FT PRDEEEEMMEVSLDATKRHAVLMGLQPATEYIVNLVAVHGTVTSEPIVGSITTGIDPP"
#
# SQ Sequence 106319 BP; 33624 A; 23004 C; 22373 G; 27318 T; 0 other;
# gatcatttta ccttgtacgt atcaaaacat cactgtgtac cccatacata tgtacaatta 60
# ttattcatca tttagaaaat aaaattaaaa aataaataaa tctctactca cccaggccct 120
# tacacagtgg catttctaga accttatttc cttttgccaa cggagcctcc attccaggct 180
#
# This can also return the sizes of sequences rather than seqs.
#
# This ignores any dup entrynames coming later.
#
# Example : %out = %{&open_embl_files(@ARGV)};
# %out2=%{&open_embl_files('seq.fa', \%index)};
# %out3=%{&open_embl_files('seq.fa', \%range)};
# %seq=%{&open_embl_files($PDB40_FASTA, \@seq_to_fetch)};
#
# while @ARGV at prompt was: 'GMJ.pep MJ0084'
#
# Keywords : open_fasta, open_fa_files, open_FASTA_files,
# Options : Seq name to fetch the specified seq only.
# as open_embl_files.pl MY_SEQ_NAME Swissprot.fasta
# -d for giving back desc as well as the name. so it
# gives 'HI0002 This is the description part'
# as the key
# If you put hash which is like('seq_name', ['20-30', '30-44',..])
# it will produce hash which has got:
# ( seq_name_20-30 'asdfasdfasdfasdfasd',
# seq_name_30-44 'kljkljkjkjljkjljkll',
# .... .... )
# -s for returning sequence size only
# $reverse_seq=r by r ## to reverse seq.
# Category :
# Version : 4.1
#--------------------------------------------------------------------
sub open_embl_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my (%sequence, %na_sequence, %HASH, @Keys, $seq_found1,
$S_start, $S_end, $seq_found,
$present_seq, @seq_Names, %Sizes, $bare_seq_name, $fasta_seq_idx_file,
%seq_fragments, $get_protein_seq, $reverse_seq, $translation_found);
$get_nucleic_acid_seq='n';
if($char_opt=~/r/){ $reverse_seq='r' }
if(@file<1){
print "\n \@file has less than 1 elem. There is no fileinput for open_embl_files\n";
exit
}
print "\n# \@file are @file\n" if $verbose;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (1) opening FASTA files (NORmal, when no perl index pos number is given)
#_______________________________________________________________________
for($i=0; $i< @file; $i++){
unless(-s $file[$i]){ print "\n# (E) File not found $file[$i]\n"; next; } ## this is essential as handle_arguments has a problem
my($entry_found, $name, $matched);
my($input_file) = ${$file[$i]} || $file[$i];
if($debug eq 1){ print "\n open_embl_files: Inputfile is $input_file\n" };
unless (-e $input_file){
print chr(7);
print "\n\n\t This is sub open_fas_files in $0 \n\n";
print "\t Fatal: The input file $input_file is not in the directory \n";
}
open(FILE_1,"$input_file");
if(@hash >=1){ ## if seq names are given in hash
for($h=0; $h< @hash; $h++){
@string=(@string, keys %{$hash[$h]});
}
}
@string=sort @string;
$num_of_seq_to_fetch=@string;
if(@string > 0){
print "\n# open_embl_files(normal fasta fetch): \$num_of_seq_to_fetch is $num_of_seq_to_fetch\n";
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (2.1) when seq to fetch is given by \@sequences format
#_______________________________________________________________________
if( @_ > 1 and @string > 0 ){
print "\n# open_embl_files is fetching sequences from \$input_file= $input_file\n";
%sequence=%{&fetch_sequence_from_db($input_file, \@string)};
print "\n# $fasta_seq_idx_file file is made by open_embl_files(fetch_sequence_from_db), you may remove it\n";
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (2.2) When seq names NOT given, fetches all (THE DEFAULT)
#____________________________________________________________
else{
$translation_found=0;
while(<FILE_1>){
if(/^ *ID[\t ]+(\S+)/){
$seq_name=$1;
}elsif($get_protein_seq and /^ *FT[\t ]+\/translation\=\"(\w+)/){
$translation_found=1;
$seq .=$1;
}elsif($get_protein_seq and $translation_found and /^ *FT[\t ]+(\w+) *$/){
$seq .=$1;
}elsif($get_protein_seq and $translation_found and /^ *FT[\t ]+(\w*)\"/){
$seq .=$1;
$translation_found=0;
}elsif($get_nucleic_acid_seq and /^ *SQ[\t ]+/){
$nucleic_acid_found=1;
}elsif($nucleic_acid_found and /^[\t ]+([\S+ ]+)[\t ]+\d*/){ # \d* NOT \d+ as some don't have digits
$na_seq.=$1;
}elsif($nucleic_acid_found and /^ *\/\//){
last;
}
}# end of while
$sequence{$seq_name}=$seq;
$seq=~s/[\t ]+//g;
$na_seq=~s/[\t ]+//g if $get_nucleic_acid_seq;
$na_seq=~tr/[a-z]/[A-Z]/;
$na_sequence{$seq_name}=$na_seq if $get_nucleic_acid_seq;
$seq=$na_seq='';
}
close FILE_1;
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`~~~~~~~~~~~~~~~~~~~~~`
# (3) When ranges information is given(via \@range), seq in those ranges are returned
#______________________________________________________________________________________
if(defined(@range)){
%seq_fragments=%{&get_seq_fragments(\%sequence, \@range)};
if($reverse_seq){
%seq_fragments=%{&reverse_sequences(\%seq_fragments)};
}
return(\%seq_fragments);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (4) When only size is asked with -s option
#_____________________________________________________________________________
elsif($char_opt=~/s/){ # when SIZE(length of seq) return only option is set
#@seq_Names=keys %sequence;
#for($i=0; $i<@seq_Names; $i++){
# $Sizes{$seq_Names[$i]}=length($sequence{$seq_Names[$i]});
# }
return(\%Sizes);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (5) when hash which has range info is given(@range should not be defined)
#_____________________________________________________________________________
elsif(@hash >=1){
for($h=0; $h< @hash; $h++){
my %hash=%{$hash[$h]};
my @Keys=keys %hash;
for($k=0; $k<@Keys; $k++){
if(defined($hash{$Keys[$k]})){
($S_start, $S_end)=$hash{$Keys[$k]}=~/(\d+)\-(\d+)/;
$sequence{$Keys[$k]}=substr($sequence{$Keys[$k]}, ($S_start-1), ($S_end-$S_start));
}
}
}
if($reverse_seq){ %sequence=%{&reverse_sequences(\%sequence)}; }
if($get_nucleic_acid_seq){
return(\%na_sequence);
}else{
return(\%sequence);
}
}else{
if($reverse_seq){ %sequence=%{&reverse_sequences(\%sequence)}; }
if($get_nucleic_acid_seq){
return(\%na_sequence);
}else{
return(\%sequence);
}
}
}
#________________________________________________________________________
# Title : open_fasta_files
# Usage : %fasta_seq=%{&open_fasta_files($fasta_file, ['MJ0084'])};
# if you put additional seq name as MJ0084 it will
# fetch that sequence only in the database file.
#
# %out=%{&open_fasta_files(@ARGV, \%index)};
# while %index has (seq indexpos seq2 indexpos2,,,)
# In this case, the fasta file should have xxxx.fa format
#
# Function : open fasta files and put sequences in a hash
# If hash(es) is put which has sequence names and seek position
# of the index file, it searches the input FASTA file to
# fetch at that seek position. This is useful for Big fasta DBs
# If the seq name has ranges like XXXXXX_1-30, it will only
# return 1-30 of XXXXXX sequence.
#
# FASTA sequence file format is like this;
#
# > 1st-seq
# ABCDEFGHIJKLMOPABCDEFGHIJKLMOPABCDEFGHIJKLMOPABCDEFG
# > 2nd.sequ
# ABCDEFGHIJKLMOYYUIUUIUIYIKLMOPABCDEFGHIJKLMOPABCDEFG
# >owl|P04439|1A03_HUMAN HLA CLASS I HISTOCOMPATIBILITY ANTIGEN, A-3 ALPHA CHAIN PRECURSOR....
# MARGDQAVMAPRTLLLLLSGALALTQTWAGSHSMRYFFTSVSRPGRGEPRFIAVGYVDDT
#
# This can also return the sizes of sequences rather than seqs.
#
# This ignores any dup entrynames coming later.
#
# Example : %out = %{&open_fasta_files(@ARGV)};
# %out2=%{&open_fasta_files('seq.fa', \%index)};
# %out3=%{&open_fasta_files('seq.fa', \%range)};
# %seq=%{&open_fasta_files($PDB40_FASTA, \@seq_to_fetch)};
#
# while @ARGV at prompt was: 'GMJ.pep MJ0084'
#
# Keywords : open_fasta, open_fa_files, open_FASTA_files,
# Options : Seq name to fetch the specified seq only.
# as open_fasta_files.pl MY_SEQ_NAME Swissprot.fasta
# -d for giving back desc as well as the name. so it
# gives 'HI0002 This is the description part'
# as the key
# If you put hash which is like('seq_name', ['20-30', '30-44',..])
# it will produce hash which has got:
# ( seq_name_20-30 'asdfasdfasdfasdfasd',
# seq_name_30-44 'kljkljkjkjljkjljkll',
# .... .... )
# -s for returning sequence size only
# $reverse_seq=r by r ## to reverse seq.
# Category :
# Version : 4.1
#--------------------------------------------------------------------
sub open_fasta_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my (%sequence, %HASH, @Keys, $seq_found1, $S_start, $S_end, $seq_found,
$present_seq, @seq_Names, %Sizes, $bare_seq_name, $fasta_seq_idx_file,
%seq_fragments, $reverse_seq);
if($char_opt=~/r/){ $reverse_seq='r' }
if(@file<1){
print "\n \@file has less than 1 elem. There is no fileinput for open_fasta_files\n";
exit
}
print "\n# \@file are @file\n" if $verbose;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~``~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# (1) When perl file INDEX pos info is given in hash, this speeds up
#_________________________________________________________________________________
for($d=0; $d < @hash; $d++){
my ($sequence, $NAME, $range_start, $range_leng);
%HASH=%{$hash[$d]};
my @Keys=keys %HASH; ## <<< NOTE it is @Keys, not @keys
for($f=0; $f< @file; $f++){
#====== It must be xxxx.fa format =======
unless($file[$f]=~/\S\.fa[sta]?$/){
print "\n# open_fasta_files: \$file\[\$f\] does not have fasta extension, skipping"; next; }
open(FASTA, $file[$f]);
F0: for($e=0; $e< @Keys; $e++){
my $sequence;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# When seq name has range attachment, it handles
#________________________________________________
if($Keys[$e]=~/^(\S+)_(\d+)\-(\d+)/){
$NAME=$1;
$range_start=$2-1; ## to fit in substr function
$range_leng =$3-$2+1; ## to fit in substr
}else{
$NAME=$Keys[$e];
}
if($HASH{$Keys[$e]}=~/^(\d+)$/){
splice(@hash, $d, 1);
$d--;
splice(@file, $f, 1);
$f--;
seek(FASTA, $1-220, 0); # -220 is necessary
while(<FASTA>){
if( /^\> *$NAME/ or
/^\> *owl\|\S+\|$NAME/){ # to handle ">owl|P04439|1A03_HUMAN HLA CLASS I HISTOCOMPATIBILITY
$seq_found1=1;
}elsif(/^(\w+)$/ and $seq_found1==1){ $sequence .=$1;
}elsif(/^\> *\S+/ and $seq_found1==1){
#======= When range is defined, take only the ranged part==================
if($range_start =~/\d+/){
$sequence{$Keys[$e]}=substr($sequence, $range_start, $range_leng);
}else{ $sequence{$Keys[$e]}=$sequence; }
$range_start='';
$sequence='';
$seq_found1=0; next F0;
}
}
}
}
}
}
# (2.1) when seq to fetch is given by \@sequences format
#_______________________________________________________________________
if( @_ > 1 and @string > 0 ){
print "\n# open_fasta_files is fetching sequences from \$input_file= $input_file\n";
%sequence=%{&fetch_sequence_from_db($input_file, \@string)};
print "\n# $fasta_seq_idx_file file is made by open_fasta_files(fetch_sequence_from_db), you may remove it\n";
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (2.2) When seq names NOT given, fetches all (THE DEFAULT)
#____________________________________________________________
else{
while(<FILE_1>){ # file1 needs to be xxxx.fasta for the moment, automatic later
if(/^> *gi\|\d+\|\S+\|(\S+)\|.*/){ ## for >gi|1669546|dbj|D84107|D84107 Human mRNA for Werner syndrome-1/type 1, complete cds
if($char_opt=~/[\-]?d/i){ # To add the description
$name=$_; # entire line becomes the name of the seque.
}else{
if( $sequence{$1} ){
#------- To avoid identical entry reading repeatedly -----
print "\n# I am open_fasta_files: $1 seems to be the same as previous entry, ERROR??\n";
$entry_found=0;
}else{ $name=$1; $entry_found=1; }
}
}elsif(/^> *owl\|\S*\|(\S+)/ or /^> *pir\|\S*\|(\S+)/){ ## for ">owl|P04439|1A03_HUMAN HLA CLASS I HISTOCOMPATIB
if($char_opt=~/[\-]?d/i){ # To add the description
$name=$_; # entire line becomes the name of the seque.
}else{
if( $sequence{$1} ){
#------- To avoid identical entry reading repeatedly -----
print "\n# I am open_fasta_files: $1 seems to be the same as previous entry, ERROR??\n";
$entry_found=0;
}else{ $name=$1; $entry_found=1; }
}
}elsif(/^> {0,5}([\w\-\.]+) *.*$/){
if($char_opt=~/[\-]?d/i){ $name=$_; # To add the description
}else{
if( $sequence{$1} ){ # check if the entry already exists
print "\n# $1 seems to be the same as previous entry, ERROR??\n";
$entry_found=0;
}else{ $name=$1; $entry_found=1; }
}
}elsif(/^([\w\.\- ]+)$/ and $entry_found == 1){
$matched=$1; $matched=~s/ //g;
if($char_opt=~/s/){ # when SIZE(length of seq) return only option is set
$Sizes{$name}+=length($matched);
#print $Sizes{$name}, "\n";
}else{
$sequence{$name}.= $matched if defined($name);
}
}elsif(/^$/){ next;
}else{ $entry_found=0; } ## this is when rubbish is matched
}# end of while
}
close FILE_1;
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`~~~~~~~~~~~~~~~~~~~~~`
# (3) When ranges information is given(via \@range), seq in those ranges are returned
#______________________________________________________________________________________
if(defined(@range)){
%seq_fragments=%{&get_seq_fragments(\%sequence, \@range)};
if($reverse_seq){
%seq_fragments=%{&reverse_sequences(\%seq_fragments)};
}
return(\%seq_fragments);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (4) When only size is asked with -s option
#_____________________________________________________________________________
elsif($char_opt=~/s/){ # when SIZE(length of seq) return only option is set
#@seq_Names=keys %sequence;
#for($i=0; $i<@seq_Names; $i++){
# $Sizes{$seq_Names[$i]}=length($sequence{$seq_Names[$i]});
# }
return(\%Sizes);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
# (5) when hash which has range info is given(@range should not be defined)
#_____________________________________________________________________________
elsif(@hash >=1){
for($h=0; $h< @hash; $h++){
my %hash=%{$hash[$h]};
my @Keys=keys %hash;
for($k=0; $k<@Keys; $k++){
if(defined($hash{$Keys[$k]})){
($S_start, $S_end)=$hash{$Keys[$k]}=~/(\d+)\-(\d+)/;
$sequence{$Keys[$k]}=substr($sequence{$Keys[$k]}, ($S_start-1), ($S_end-$S_start));
}
}
}
if($reverse_seq){ %sequence=%{&reverse_sequences(\%sequence)}; }
return(\%sequence);
}else{
if($reverse_seq){ %sequence=%{&reverse_sequences(\%sequence)}; }
return(\%sequence);
}
}
#________________________________________________________________________
# Title : msf_permute_hash_write
# Usage : &msf_permute_hash_write(\%hash, $group_name); # void
# Function : gets 2 references (one for %hash the other for group $name)
# uses &msf_permute_array_write(\%hash, \$group_name)
# the second arg is for output file name. can be anything.
# Example :
# Warning :
# Keywords :
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.0
#--------------------------------------------------------------------
sub msf_permute_hash_write{ my(%hash)=%{$_[0]}; my($group_name)=${$_[1]};
my(@array)= keys %hash; &msf_permute_array_write (\@array, $group_name );
}
#________________________________________________________________________
# Title : msf_permute_array_write
# Usage : &msf_permu_array_write(\%hash, \$group_name); # void
# Function :
# the second arg is for output file name. can be anything.
print C3SS_FILE "PFRMAT SS\n";
print C3SS_FILE "TARGET $base\n";
print C3SS_FILE "AUTHOR 1143-9512-5546\n";
print C3SS_FILE "REMARK secondary struc prediction, _ means NO prediction\n";
print C3SS_FILE "METHOD FBRV postprocessing method\n";
print C3SS_FILE "MODEL 1\n";
@position=sort{$a<=>$b} keys %prdl;
for($i=0; $i<@position; $i++){
$residue=$prdl{$position[$i]}->[0];
if($residue =~/\W/){ next } # skipping non residue(gap)
$reliability=substr($prdl{$position[$i]}->[2], 0, 5);
$sec_str="\U$prdl{$position[$i]}->[1]";
if($reliability > 0.999){
$reliability=0.999;
}
if($reliability <=0){ $sec_str='_' }
print C3SS_FILE "$residue $sec_str $reliability\n"
}
print C3SS_FILE "END\n";
return(\$out_c3ss_file);
}
#________________________________________________________________________
# Title : write_pred_files
# Usage :
# Function :
# Example : &write_pred_files(\%gapped_av_for_back_pred, $final_output_pred_name,
# $graphical_rep_of_str, "$put_reliability_line");
#
# Warning :
# Keywords : write_predator_short_out_file, write_pred_file, write_prd_file
# Options :
# $put_reliability_line=r by r
# $omit_coil_region=c by c
# $protein_name= by n=
# $graphical_rep_of_str=g by g
# $show_on_screen_only=s by s
# $seq_block_size= by b=
# Returns :
# Argument :
# Category :
# Version : 1.4
#--------------------------------------------------------------------
sub write_pred_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#$| = 1;
my($sort_seq_names, $string, $string_leng, $na,$out_file_name_provided,
$write_pure_seq_only, $show_coil_region, %one_to_three, $each_reliability,
$each_seq_fragment, $each_str_fragment, $end_seq_point, $each_reliability_frag,
$protein_name, $omit_coil_region, $ext, $graphical_rep_of_str,
$show_on_screen_only, $seq_block_size);
my($output_file) ='default_out.pred'; ### when no output file name is given, this is used
$protein_name='pred_temp_name';
my $seq_block_size=80;
if($char_opt=~/r/){ $put_reliability_line='r' } ## putting additional reliability line
if($char_opt=~/c/){ $omit_coil_region='c' } ## putting additional reliability line
if($char_opt=~/g/){ $graphical_rep_of_str='g' } ## putting additional reliability line
if($char_opt=~/s/){ $show_on_screen_only='s' } ## putting additional reliability line
if($vars{'n'}=~/(\S+)/){ $protein_name=$1; }
if($vars{'b'}=~/(\d+)/){ $seq_block_size=$1; } ## final seq block size
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If any input file has the prdl extension, I assume it is the output name
#__________________________________________________________________________
$ext='pred';
for($i=0; $i< @file; $i++){
if($file[$i]=~/$ext/){
$output_file=$file[$i];
}
}
%one_to_three = (
"A", "ALA", "C", "CYS", "D", "ASP", "E", "GLU", "F", "PHE",
"G", "GLY", "H", "HIS", "I", "ILE", "K", "LYS", "L", "LEU",
"M", "MET", "N", "ASN", "P", "PRO", "Q", "GLN", "R", "ARG",
"S", "SER", "T", "THR", "V", "VAL", "W", "TRP", "Y", "TYR"
);
if(@file>0){
$output_file = $file[0];
$out_file_name_provided=1;
}elsif($protein_name){ $output_file="$protein_name\.pred"
}else{ $output_file='default_out.pred'; }
for ($n=0 ; $n < @hash; $n ++){
my($reliability, $position, @keys, $residue, $sec_str, $protein_name);
my %hash=%{$hash[$n]};
@keys= sort {$a <=> $b} keys %hash; # <------- old style
$protein_name=$hash{0}->[3];
unless($show_on_screen_only){
open (PRED_WRITE,">$output_file"); # $string is the seq string.
print PRED_WRITE "# written by $0 . predator sec. str prediction based\n";
print PRED_WRITE "\n\> $protein_name\n";
}else{
print "\n# written by $0 . predator sec. str prediction based\n";
print "\n\> $protein_name\n";
}
for ($i=0; $i < @keys; $i++){
$position= $keys[$i];
if($put_reliability_line){
$each_reliability=int($hash{$position}->[2]*10);
if($each_reliability > 9){
$each_reliability= 9;
$sec_str .= "\U$hash{$position}->[1]";
}elsif($each_reliability < 4 and $each_reliability > 0){
$sec_str .= $hash{$position}->[1];
}elsif($each_reliability >= 4){
$sec_str .= "\U$hash{$position}->[1]";
}elsif($each_reliability <= 0){ ## If the reliability is 0, Don't show the sec. str
$sec_str .='_';
}
$residue .= $hash{$position}->[0];
$reliability .= $each_reliability;
}else{
$residue .= $hash{$position}->[0];
$sec_str .= "\U$hash{$position}->[1]"; ## usually one of 'c h e'
}
}
if( $omit_coil_region){
$sec_str=~s/[Cc]/_/g;
}
if($graphical_rep_of_str){
$sec_str=~s/[Hh]/>/g;
$sec_str=~s/[Ee]/=/g;
$sec_str=~s/[Cc]/~/g;
$sec_str=~s/_/ /g;
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~``
# writing down the 80 residue seq fragment
#______________________________________________________
$seq_ruler_line=' . . . . . . . .';
for($i=0; $i< @keys; $i+= $seq_block_size){
$each_seq_ruler_line=substr($seq_ruler_line, 0, @keys-$i-1);
$each_seq_fragment=substr($residue, $i, $seq_block_size);
$each_str_fragment=substr($sec_str, $i, $seq_block_size);
$end_seq_point=$i+length($each_seq_fragment);
unless($show_on_screen_only){
printf PRED_WRITE ("%5s %-${seq_block_size}s %-5s\n", ' ', $each_seq_ruler_line, ' ');
printf PRED_WRITE ("%5s %-${seq_block_size}s %-5s\n", $i+1, $each_seq_fragment, $end_seq_point);
}else{
printf ("%5s %-${seq_block_size}s %-5s\n", ' ', $each_seq_ruler_line, ' ');
printf ("%5s %-${seq_block_size}s %-5s\n", $i+1, $each_seq_fragment, $end_seq_point);
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# if you want to put reliability index line
#____________________________________________
if($put_reliability_line){
$each_reliability_frag=substr($reliability, $i, $seq_block_size);
if(! $show_on_screen_only){
printf PRED_WRITE ("%5s %-${seq_block_size}s %-5s\n", ' ', $each_str_fragment, ' ');
printf PRED_WRITE ("%5s %-${seq_block_size}s %-5s\n\n", ' ', $each_reliability_frag, ' ');
}else{
printf ("%5s %-${seq_block_size}s %-5s\n", ' ', $each_str_fragment, ' ');
printf ("%5s %-${seq_block_size}s %-5s\n\n", ' ', $each_reliability_frag, ' ');
}
}else{
if(! $show_on_screen_only){
printf PRED_WRITE ("%5s %-${seq_block_size}s %-5s\n\n", ' ', $each_str_fragment, ' ');
}else{
printf ("%5s %-${seq_block_size}s %-5s\n\n", ' ', $each_str_fragment, ' ');
}
}
}
close PRED_WRITE unless $show_on_screen_only;
splice(@hash, $n, 1); $n--;
}
if( $out_file_name_provided != 1){
print "\n\n# You didnt give out file name, $output_file used\n" unless $show_on_screen_only;
}
if( -s $output_file ){
if($verbose=~/\S/){ ## if v option is given, mesg is omitted to prevent comments to a redirected output
print "\n# Sequences were written in $output_file " unless $show_on_screen_only;
}
}else{
print "\n# The size of written outfile \"$output_file\" is 0, error \n\n" unless $show_on_screen_only;
}
}
#________________________________________________________________________
# Title : write_prdl_files
# Usage :
#
# Function :
# Example :
#
# Warning :
# Keywords : write_predator_long_out_file
# Options :
# $protein_name= by n=
# Returns :
# Argument :
# Category :
# Version : 1.1
#--------------------------------------------------------------------
sub write_prdl_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#$| = 1;
my($sort_seq_names, $string, $string_leng, $na,$out_file_name_provided,
$write_pure_seq_only, %one_to_three, $protein_name, $ext);
my($output_file) ='default_out.prdl'; ### when no output file name is given, this is used
$protein_name='temp_prdl';
if($vars{'n'}=~/(\S+)/){ $protein_name=$1; }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If any input file has the prdl extension, I assume it is the output name
#__________________________________________________________________________
$ext='prdl';
for($i=0; $i< @file; $i++){
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my(@names, %hash, %hash_rv, @out_hash_ref_list, $no_gap_char_included,
@seq_order, $produce_seq_oder_info, $reverse_seq,
%redundency_check);
if($char_opt=~/n/){ $no_gap_char_included='n' }
if($char_opt=~/r/){ $reverse_seq='r' }
if($char_opt=~/o/){ $produce_seq_oder_info='o'; print "\n# (INFO) open_msf_files: \$produce_seq_oder_info is set\n"; }
for($i=0; $i< @file; $i++){
my(%redundency_check, %hash_rv);
print "\n# (INFO) Opening $file[$i] in @file\n";
open(FILE_1, "$file[$i]");
while(<FILE_1>){ # file1 needs to be xxxx.msf
if(/^\S+[\t ]*$/ or /^\#/||/^\-+/ ){ next; }
if(/^(\S+)[\t ]+([\.\-\w ]+) *$/){
$n=$1; $s=$2; $s=~s/ //g;
if($no_gap_char_included){ $s=~s/\.//g; }
$hash{$n}.= $s;
if($produce_seq_oder_info){
push(@seq_order, $n) unless $redundency_check{$n};
$redundency_check{$n}=$n;
}
}
}
if($reverse_seq){ print "\n# (INFO) I am reversing seq \n";
%hash_rv=%{&reverse_sequences(\%hash)};
push(@out_hash_ref_list, \%hash_rv);
}else{
push(@out_hash_ref_list, \%hash);
}
close(FILE_1);
}
if(@out_hash_ref_list == 1 ){
if($produce_seq_oder_info){
return($out_hash_ref_list[0], \@seq_order);
}else{ return($out_hash_ref_list[0]) }
}elsif(@out_hash_ref_list > 1){
if($produce_seq_oder_info){
return(\@out_hash_ref_list, \@seq_order);
}else{
return(@out_hash_ref_list);
}
}
}
#__________________________________________________________________________
# Title : open_hmmls_files
# Usage : %out=%{&open_hmmls_files(\@file)};
# Function : hmmls matches the full length model to target seq. while, hmmfs
# does for fragments as well.
# Example :
# Keywords :
# Options :
# t=$thresh for bits score threshold
# e=$evalue_thresh for bits score threshold
# r for adding ranges
# m for making MSP file format output
# E=Enguiry_name for specifying enquiry seq name rather than 'HMM', the default
# Returns :
# Argument :
# Version : 1.5
#----------------------------------------------------------------------------
sub open_hmmls_files{
#"""""""""""""""""< handle_arguments{ head Ver 4.1 >"""""""""""""""""""
my(@A)=&handle_arguments(@_);my($num_opt)=${$A[7]};my($char_opt)=${$A[8]};
my(@hash)=@{$A[0]};my(@file)=@{$A[4]};my(@dir)=@{$A[3]};my(@array)=@{$A[1]};
my(@string)=@{$A[2]};my(@num_opt)=@{$A[5]};my(@char_opt)=@{$A[6]};
my(@raw_string)=@{$A[9]};my(%vars)=%{$A[10]};my(@range)=@{$A[11]};
my($i,$j,$c,$d,$e,$f,$g,$h,$k,$l,$m,$n,$o,$p,$q,$r,$s,$t,$u,$v,$w,$x,$y,$z);
if($debug==1){print "\n\t\@hash=\"@hash\"
\@raw_string=\"@raw_string\"\n\t\@array=\"@array\"\n\t\@num_opt=\"@num_opt\"
\@char_opt=\"@char_opt\"\n\t\@file=\"@file\"\n\t\@string=\"@string\"\n" }
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
my($query, $thresh, $evalue_thresh, $read_point_found); # $read_point_found is not used in this sub
$thresh = 10; # default threshold
$evalue_thresh=5;
if($vars{'t'}=~/\S/){ $thresh = $vars{'t'} }
if($vars{'e'}=~/\S/){ $thresh = $vars{'e'} }
if($vars{'E'}=~/\S/){ $enquiry_name = $vars{'E'} }
if(@file<1){
print "\n# open_hmmls_files needs one or more files !!\n";
}else{
print "\n# open_hmmls_files: threshold used = $thresh \n";
}
print "\n# (INFO) open_hmmls_files: score \$thresh is $thresh\n";
for($i=0; $i< @file; $i++){
$input_file=$file[$i];
open(HMMLS, "$input_file");
while(<HMMLS>){
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# HMMER version 1
#_____________________________________________________
if(/^ {0,2}(\S+) +\S+ +\S+ +(\d+) +\S+ +(\d+) +\S+ +(\S+) *$/){
if($1 > $thresh){
$score=$1;
$matched=$4;
$from=$2;
$to=$3;
if($char_opt=~/r/){
if($char_opt=~/m/){ ## convert to msp file format
$hash{"$matched"}="$score $score $from $to HMM $from $to $matched\n";
}else{
$hash{"$matched\_$from\-$to"}=$score;
}
}else{
if($char_opt=~/m/){ ## convert to msp file format
$hash{"$matched"}="$score $score $from $to HMM $from $to $matched\n";
my($input_file) = ${$file[$i]} || $file[$i];
if($debug eq 1){ print "\n open_seq_files: Inputfile is $input_file\n" };
open(FILE_1,"$input_file");
@lines=<FILE_1>;
close FILE_1;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# when hash is directly given for selected sequence names
#_____________________________________________________________
if(@hash >=1){
for($h=0; $h< @hash; $h++){
@string=(@string, keys %{$hash[$h]});
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# When there are asked seqs to fetch
#_______________________________________________________
if( (@_ > 1)&&(@string > 0) ){
for $lines (@lines){
for($j=0; $j< @string; $j++){
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If the seq name has ranges -> 1xxx_10-20
#___________________________________________
if($string[$j]=~/^(\S+)_(\d+)\-(\d+)$/){
$name_stem=$1;
$seq_start=$2-1;
$seq_leng=$3-$2+1;
if($lines=~/^\>?($name_stem) +(\d+) +(\w+)$/i){
$seq_in_range=substr($3, $seq_start, $seq_leng);
$sequence{$string[$j]}=$seq_in_range;
splice(@string, $j, 1); ## for optimization
$j--;
}
}else{
if($lines=~/^\>?($string[$j]) +(\d+) +(\w+)/i){
if($min_size < $2){
splice(@string, $j, 1);
$j--; ## for optimization
$sequence{$1} = $3;
}
}
}
}
}
return(\%sequence);
}else{ # getting all seq in the given file(s)
for(@lines){ # file1 needs to be xxxx.fasta for the moment, automatic later
if(/^\>?(\S+) +(\d+) +(\w+)$/){
if($min_size < $2){
$sequence{$1} = $3;
}
}elsif(/^\>?(\S+)[\t ]+(\D\D+) *$/){ ## plain 'seqname ASDASDFAASSEQUECN'
$sequence{$1} = $2;
}else{ next }
}
}
close FILE_1;
}
if(defined(@range)){
%seq_fragments=%{&get_seq_fragments(\%sequence, \@range)};
return(\%seq_fragments);
}elsif($char_opt=~/s/i){ # when SIZE return only option is set
@seq_Names=keys %sequence;
for($i=0; $i < @seq_Names; $i++){
$Sizes{$seq_Names[$i]}=length($sequence{$seq_Names[$i]});
}
return(\%Sizes);
}elsif(@hash >=1){
for($h=0; $h< @hash; $h++){
my %hash=%{$hash[$h]};
my @Keys=keys %hash;
for($k=0; $k<@Keys; $k++){
if(defined($hash{$Keys[$k]})){
($S_start, $S_end)=$hash{$Keys[$k]}=~/(\d+)\-(\d+)/;
$sequence{$Keys[$k]}=substr($sequence{$Keys[$k]}, ($S_start-1), ($S_end-$S_start));
}
}
}
return(\%sequence);
}else{
return(\%sequence);
}
}
#_________________________________________________________________________________
# Title : open_sso_files
# Usage : @sso=@{&open_sso_files(@file, $add_range, $add_range2, "u=$upper_expect_limit",
# "l=$lower_expect_limit", "m=$margin", $new_format)};
# Function : This reads the parseable( -m 10 option)
# and non-parseable form of ssearch program output
# If you give 5 files, it produces 5 hashes as a ref of array.
# This understands xxxx.gz files.
# This reads FASTA -m 10 output, too.
# Example :
# 717 0 0.343 16 373 EC1260_16-373 74 434 YBL6_YEAST_74-434
# 348 9e-16 0.500 113 233 EC1260_113-233 27 146 YDBG_ECOLI_27-146
# 472 2.9e-08 0.271 13 407 EC1260_13-407 148 567 YHJ9_YEAST_148-567
# 459 1.9e-22 0.260 1 407 EC1260_1-407 65 477 YLQ6_CAEEL_65-477
# 452 4.5e-14 0.275 1 407 EC1260_1-407 103 537 YSCPUT2_103-537
# 1131 0 0.433 1 407 EC1260_1-407 112 519 ZMU43082_112-519
#
# Input SSO file example)-> below
#
# >>MG032 ATP-dependent nuclease (addA) {Bacillus subtilis (666 aa)
# Z-score: 88.3 expect() 1.9
# Smith-Waterman score: 77; 27.143% identity in 70 aa overlap
#
# 30 40 50 60 70 80
# MJ0497 RSAGSKGVDLIAGRKGEVLIFECKTSSKTKFYINKEDIEKLISFSEIFGGKPYLAIKFNG
# : .. ... . .:.:::. :: : ..:
# MG032 HDKVRYAFEVKFNIALVLSINKSNVDFDFDFILKTDNFSDIENFNEIFNRKPALQFRFYT
# 200 210 220 230 240 250
#
# 90 100 110 120 130
# MJ0497 EMLFINPFLLSTNGK------NYVIDERIKAIAIDFYEVIGRGKQLKIDDLI
# . :: :: ::. : ....... . ::. . :
# MG032 K---INVHKLSFNGSDSTYIANILLQDQFNLLEIDLNKSIYALDLENAKERFDKEFVQPL
# 260 270 280 290 300 310
#
for($k=0; $k<$query_leng; $k++){
$query[$k]=0;
}
$query_fasta="$query_name\.fa";
open(QUERY_FASTA, "$query_fasta") || die "\n# open_msp_files, die: No $matched_fasta. Didn't you want \'n\' opt?\n";
my @QUERY_FASTA=<QUERY_FASTA>;
for(@QUERY_FASTA){
if(/\>$query_name/){
$query_found=1; next;
}elsif(($query_found==1)&&(/^(\w{10,})$/)){
$matched_seq .=$1;
}
}
$matched_seq{$query_name}=$matched_seq;
close(QUERY_FASTA);
}elsif($6 ne $9){
if($match_name ne $9){
$match_name=$9;
$query_leng=($5 - $4 + 1);
if($query_leng < $query_length_thresh){ next }
my $matched_seq;
push(@matched_members, $9);
$matched_fasta="$match_name\.fa";
open(MATCHED_FASTA, "$matched_fasta") || die "\n# open_msp_files, die: No $matched_fasta. Didn't you want \'n\' opt?\n";
my @MATCHED_FASTA=<MATCHED_FASTA>;
for(@MATCHED_FASTA){
if(/\>$match_name/){
$found=1; next;
}elsif(($found==1)&&(/^(\w{10,})$/)){
$matched_seq .=$1;
}
}
$matched_seq{$match_name}=$matched_seq;
}
}
}
}
}else{ # default output. This filters with given E and Score thresholds
$simple_default_out_flag++;
F4: for($j=0; $j<@MSP; $j++){
#""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# $1 $2 $3 $4 $5 $6 $7 $8
# 171 41.18 6 73 HI1690 9 76 HI0736 sodium...
#,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
if($MSP[$j]=~/^ *(\d+) +(\S+) +\S* *(\d+) +(\d+) +(\S+) +(\d+) +(\d+) +(\S+) *(.*)/i){
my $query_length = $4-$3;
if($2 > $Evalue_thresh or $1 < $Score_thresh
or $query_length < $query_length_thresh ){ next F4;
}else{
$simple_default_output{$input_file} .= $MSP[$j];
}
}
}
}
} #<<<<<-- for($i=0; $i< @file; $i++){
if($simple_default_out_flag >=1){ %matched_seq=%simple_default_output; }
if($simple_default_out_flag > 1){ print "\n# You have asked me to open more than one file. I merge output into one hash\n"; }
if(defined(@range) && ($char_opt !~ /n/)){ ## m opt is for matched seq names only
%seq_fragments=%{&get_seq_fragments(\%matched_seq, \@range)};
return(\%seq_fragments);
}elsif($char_opt=~/s/){ # when SIZE return only option is set
@seq_Names=keys %matched_seq;
for($i=0; $i<@seq_Names; $i++){
$Sizes{$seq_Names[$i]}=length($matched_seq{$seq_Names[$i]});
}
return(\%Sizes);
}else{
if($char_opt=~/n/){
return(\%matched_seq_names);
}else{
return(\%matched_seq);
}
}
}
#________________________________________________________________________
# Title : open_dssp_files
# Usage : (*out, *out2) = @{&open_dssp_files(\$inputfile1, \$inputfile2, \$H, \$S,,,,)};
# (@out) = @{&open_dssp_files(\$inputfile1, \$inputfile2, \$H, \$S,,,,)};
# Function : open dssp files and put sequences in a hash(s)
# It can take options for specific secondary structure types. For example,
# if you put an option $H in the args of the sub with the value of 'H'
# open_dssp_files will only read secondary structure whenever it sees 'H'
# in xxx.dssp file ignoring any other sec. str. types.
# If you combine the options of 'H' and 'E', you can get only Helix and long
# beta strand sections defined as segments. This is handy to get sec. str. segments
# from any dssp files to compare with pdb files etc.
# With 'simplify' option, you can convert only all the 'T', 'G' and 'I' sec. to
# 'H' and 'E'.
# Example :
# Warning : 6taa.dssp and 6taa are regarded as the same.
# Keywords :
# Options : H, S, E, T, I, G, B, P, C, -help
# $H = 'H' by -H or -h or H or h # to retrieve 4-helix (alpha helical)
# $S becomes 'S' by -S or -s or S or s # to retrieve Extended strand, participates in B-ladder
# $E becomes 'E' by -E or -e or E or e # to retrieve residue in isolated Beta-bridge
# $T becomes 'T' by -T or -t or T or t # to retrieve H-bonded turn
# $I becomes 'I' by -I or -i or I or i # to retrieve 5-helix (Pi helical) segment output
# $G becomes 'G' by -G or -g or G or g # to retrieve 3-helix (3-10 helical)
# $B becomes 'B' by -B or -b or B or b # to retrieve only B segment
# $simplify becomes 1 by -p or P or -P, p
# $comm_col becomes 'c' by -c or c or C or -C or common
# $HELP becomes 1 by -help # for showing help
#
# Returns : (*out, *out2) or (@out_array_of_refs)
# Argument : files names like (6taa, 6taa.dssp) If you put just '6taa' without extension, it
# searches if there is a '6taa.dssp' in both PWD and $DSSP env. set directory.
# ---------- Example of dssp ---
# **** SECONDARY STRUCTURE DEFINITION BY THE PROGRAM DSSP, VERSION JUL
# REFERENCE W
# HEADER RIBOSOME-INACTIVATING PROTEIN 01-JUL-94 1MRG
# COMPND ALPHA-MOMORCHARIN COMPLEXED WITH ADENINE
# SOURCE BITTER GOURD (CUCURBITACEAE MOMORDICA CHARANTIA) SEEDS
# AUTHOR Q
# 246 1 0 0 0 TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN) .
# 112 95.0 ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2) .
( run in 1.121 second using v1.01-cache-2.11-cpan-364913b4093 )