2179 lines
81 KiB
Perl
Executable File
2179 lines
81 KiB
Perl
Executable File
#!/usr/bin/env perl
|
|
|
|
# ps_scan - a PROSITE scanning program
|
|
#
|
|
# Revision: 1.91
|
|
#
|
|
# Copyright (C) 2001-2020 SIB Swiss Institute of Bioinformatics
|
|
# Authors:
|
|
# edouard.decastro@sib.swiss
|
|
# Alexandre Gattiker
|
|
# Beatrice Cuche (evaluated_by post-processing)
|
|
# Lorenzo Cerutti (FTREP method)
|
|
#
|
|
# Contributions:
|
|
# Lorenza Bordoli (repeat method)
|
|
|
|
# This program is free software; you can redistribute it and/or modify it under
|
|
# the terms of the GNU General Public License as published by the Free Software
|
|
# Foundation; either version 2 of the License, or (at your option) any later
|
|
# version.
|
|
#
|
|
# This program is distributed in the hope that it will be useful, but WITHOUT
|
|
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
|
# details.
|
|
#
|
|
# You should have received a copy of the GNU General Public License along with
|
|
# this program; if not, write to the Free Software Foundation, Inc., 59 Temple
|
|
# Place - Suite 330, Boston, MA 02111-1307, USA.
|
|
|
|
|
|
use 5.005_03; # uses perl >= v5.5
|
|
use IO::File;
|
|
use Carp qw(confess cluck);
|
|
use vars qw(@ISA $VERSION $errpos $errstr);
|
|
use strict;
|
|
use Data::Dumper;
|
|
|
|
|
|
################################################################################
|
|
# integrated subs taken from Prosite.pm module:
|
|
|
|
|
|
# scan a sequence with a perl pattern
|
|
sub scanPattern {
|
|
my ($pattern, $sequence, $behavior, $max_x, $opt_miniprofiles) = @_;
|
|
$behavior ||= 0;
|
|
my $allowOverlap = !($behavior & 1);
|
|
my $allowInclude = $behavior & 2 ;
|
|
$max_x ||= 0;
|
|
my @hits;
|
|
my $pos = 0;
|
|
my @comb = $pattern =~ /\((.*?)\)/g;
|
|
if ($pattern) {
|
|
my $prevstop = -1;
|
|
my @tok;
|
|
my $nter_anchor = $pattern=~/^\^/ ? 1 : 0;
|
|
while (@tok = (substr($sequence, $pos) =~ /^(.*?)($pattern)/)) {
|
|
my $prematch = shift @tok;
|
|
my $subseq = shift @tok;
|
|
my $length = length $subseq;
|
|
my $number_x = 0;
|
|
if (@tok == @comb and @tok) {
|
|
$subseq = "";
|
|
for (my $i = 0; $i < @tok; $i++) {
|
|
if ($comb[$i] =~ /\./) {
|
|
$tok[$i] =~ tr/A-Z/a-z/;
|
|
}
|
|
# don't count X's that match pattern elements which exclude
|
|
# certain AA's
|
|
elsif ($comb[$i] =~ /^\[\^/) {
|
|
}
|
|
# count number of X's matching a non-x position in the
|
|
# pattern
|
|
elsif (my $x_count = $tok[$i] =~ tr/Xx/Xx/) {
|
|
$number_x += $x_count;
|
|
}
|
|
my $tok = $comb[$i] =~ /\./ ? lc($tok[$i]) : $tok[$i];
|
|
$subseq .= $tok;
|
|
if (my @numbers = $comb[$i] =~ /(\d+)/g) {# insertion
|
|
my $biggest = pop @numbers;
|
|
$subseq .= "." x ($biggest - length $tok);
|
|
}
|
|
}
|
|
}
|
|
elsif (@tok != @comb) {
|
|
cluck "Internal error with regular expression $pattern\n";
|
|
}
|
|
my $shift = ($length || 1) - 1;
|
|
my $stop = $pos + $length + length $prematch;
|
|
$pos = $stop;
|
|
$pos -= $shift if $allowOverlap;
|
|
# may happen if empty pattern match
|
|
last if $pos > length $sequence;
|
|
if ($length) {
|
|
if ($allowInclude or $stop > $prevstop) {
|
|
if ($number_x <= $max_x or $max_x<0) {
|
|
push @hits, [$subseq, $stop - $length + 1,
|
|
$stop, undef,
|
|
($opt_miniprofiles && !@hits ? $sequence : undef) ];
|
|
# also push sequence for ps_scan option -b
|
|
$prevstop = $stop;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
$pos++;# empty pattern match
|
|
}
|
|
last if $nter_anchor and $pos;
|
|
}
|
|
}
|
|
return \@hits;
|
|
}
|
|
|
|
sub scanProfiles {
|
|
my ($file, $level_min) = @_;
|
|
$level_min = -99 unless defined $level_min;
|
|
my $results = {};
|
|
local $/ = "\n";
|
|
my $must_open_file = !UNIVERSAL::isa($file, "GLOB");
|
|
my $pfscan_h;
|
|
if ($must_open_file) {
|
|
$pfscan_h = new IO::File $file or confess "Could not open $file : $!";
|
|
}
|
|
else {
|
|
$pfscan_h = $file;
|
|
}
|
|
my $lastac = "";
|
|
my $last_entry;
|
|
while (defined (local $_=<$pfscan_h>)) {
|
|
if (my ($id1, $level, $levelna, $nscore, $rawscore, $from, $to,
|
|
$pffrom, $pfto, $repregion, $repnumber, $ac, $de) = m/
|
|
(
|
|
>\S*
|
|
)? #ID of the profile, if option -x or -s
|
|
(?<!L=\d) #These three lines are to fix a bug in the output of pfscan if >99 matches are found: the id and the level are then pasted together with an intervening "1" (e.g. Q9LGZ9 vs. preprofile PS50079). FIXME: This still doesn't work with pfsearch -xlr
|
|
(?<!L=[-\d]\d)
|
|
(?<!L=[-\d]\d\d)
|
|
\s*
|
|
(?:L=(\S+)|(NA))? #level, if option -l
|
|
\s*
|
|
(?:(\S+)\s+)? #normalized score, unless option -r
|
|
(\S+) #raw score
|
|
\s+
|
|
pos\.
|
|
\s*
|
|
(\d+) #seq from
|
|
\s*-\s*
|
|
(\d+) #seq to
|
|
\s*
|
|
(?:
|
|
\[\s*
|
|
(\d+) #profile from
|
|
,\s*
|
|
(-\d+) #profile to
|
|
\]
|
|
\s*
|
|
)? #if option -z
|
|
(REGION\d+\s)? #repeat region, if option -m
|
|
(\d+REP\d+\s)? #repeat number, if option -m
|
|
(\S+) #profile or sequence AC and-or ID
|
|
(?:\s*(.+))? #sequence description when running pfsearch, may be absent
|
|
/x) {
|
|
# fix bug in pfsearch/pfscan which report "******" if nscore>999.999
|
|
$nscore = "999.999" if $nscore eq "*******";
|
|
$level = $level_min if !defined($level) and defined $levelna;
|
|
# 0 1 2 3 4 5 6 7 8 9 10
|
|
my $entry = ["", $from, $to, $ac, $pffrom, $pfto, $rawscore, $nscore, $level, undef, $de, []];
|
|
if ($repnumber) { push @{$results->{$lastac}->[-1]->[11]}, $entry }
|
|
else { push @{$results->{$ac}}, $entry }
|
|
$lastac = $ac;
|
|
$last_entry = $entry;
|
|
}
|
|
else {
|
|
# sequence
|
|
next unless $last_entry;
|
|
$last_entry->[0] .= $_;
|
|
}
|
|
}
|
|
if ($must_open_file) {
|
|
close $pfscan_h or confess "Error $? closing $file";
|
|
}
|
|
return $results;
|
|
}
|
|
|
|
|
|
sub prositeToRegexp_old {
|
|
local $_ = shift;
|
|
my $notGreedy = shift;
|
|
my $ungreed = $notGreedy ? "?" : "";
|
|
s/\.$//;# possible end dot from parsing
|
|
s/-//g;
|
|
1 while s/\{([^\}^\^]*)([^\^^\}][^\}]*\})/\{^$1$2/; #insert ^
|
|
s/\}/]/g;
|
|
s/\{/[/g;
|
|
s/\(/{/g;
|
|
s/\)/}$ungreed/g;
|
|
s/x/./ig;
|
|
s/B/[ND]/g;
|
|
s/Z/[QE]/g;
|
|
# PS00539 is P-R-L-[G>] which would be converted to PRL[G$], but "$" is
|
|
# not valid in a character range, we have to convert this to (?:[G]|$)
|
|
s/\[([^\[\]]*)([<>])([^\[\]]*)\]/(?:[$1$3]|$2)/g;
|
|
# do this after the previous step, so as not to mix ^ for excluding
|
|
# character classes with ^ to indicate anchoring
|
|
s/</^/g; # anchors
|
|
s/>/\$/g;
|
|
# insert () around match states and insertions
|
|
s/ (\[[^\]]*\]|[\w.]) ( \{ \d+(,\d+)? \} )/)($1$2)(/xg;
|
|
return "($_)";
|
|
}
|
|
|
|
# same, using tokenizing parser
|
|
sub prositeToRegexp {
|
|
my $string = shift;
|
|
my $notGreedy = shift;
|
|
my $ungreed = $notGreedy ? "?" : "";
|
|
my $preventX = shift;
|
|
$errstr = undef;
|
|
$errpos = undef;
|
|
my $pushback = "";
|
|
my $ntok = 0;
|
|
my $regexp = "";
|
|
my $get = sub {
|
|
$ntok++;
|
|
return ($pushback =~ s/(.)// ? $1 : ($string =~ s/(.)// ? $1 : undef));
|
|
};
|
|
while (defined (my $tok = &$get)) {
|
|
my $state;
|
|
my $not;
|
|
if ($tok eq "-") {
|
|
# ignore
|
|
}
|
|
elsif ($tok eq "[") {# RANGE
|
|
while(defined (my $tok = &$get)) {
|
|
last if $tok eq "]";
|
|
$state .= $tok;
|
|
}
|
|
}
|
|
elsif ($tok eq "{") {# negative RANGE
|
|
$not = 1;
|
|
while(defined(my $tok = &$get)) {
|
|
last if $tok eq "}";
|
|
$state .= $tok;
|
|
}
|
|
}
|
|
elsif ($tok =~ /[A-Za-z]/) {# single char
|
|
$state = $tok;
|
|
}
|
|
elsif ($tok eq "<") {
|
|
$regexp .= "^";
|
|
}
|
|
elsif ($tok eq ">") {
|
|
$regexp .= '$';
|
|
}
|
|
else {
|
|
$errstr = "Parsing error";
|
|
$errpos = $ntok-1;
|
|
return undef;
|
|
}
|
|
if (defined $state) {
|
|
# read range, e.g. "x(2,5)"
|
|
my $range;
|
|
my $range_char;
|
|
if (defined(my $tok = &$get)) {
|
|
if ($tok eq "(") {
|
|
while (defined(my $tok = &$get)) {
|
|
last if $tok eq ")";
|
|
$range .= $tok;
|
|
}
|
|
}
|
|
elsif ($tok eq "*") {# support e.g. "<{C}*>"
|
|
$range_char = $tok;
|
|
}
|
|
else {
|
|
$pushback .= $tok;
|
|
$ntok--;
|
|
}
|
|
}
|
|
|
|
if ($state =~ /x/i) {
|
|
$state = ".";
|
|
}
|
|
else {
|
|
# handle B/Z unsure amino acids both in pattern and sequence
|
|
if ($not) {
|
|
$state =~ s/B/NDB/g;
|
|
$state =~ s/Z/QEZ/g;
|
|
}
|
|
else {
|
|
$state =~ s/B/NDB/g or $state =~ s/([ND])/$1B/g;
|
|
$state =~ s/Z/QEZ/g or $state =~ s/([QE])/$1Z/g;
|
|
$state .= "X" unless $preventX;
|
|
}
|
|
}
|
|
my $mod = $1 if $state =~ s/([<>])//g;
|
|
# "$" is not valid in a character range, we have to convert this
|
|
# to (?:[GH]|$)
|
|
$regexp .= "(";
|
|
$regexp .= "(?:" if $mod;
|
|
$regexp .= "[" if length($state) > 1 or $not;
|
|
$regexp .= "^" if $not;
|
|
$regexp .= $state;
|
|
$regexp .= "]" if length($state) > 1 or $not;
|
|
$regexp .= "|" . ($mod eq "<" ? "^" : '$') . ")" if $mod;
|
|
$regexp .= "{$range}$ungreed" if defined $range;
|
|
$regexp .= "$range_char" if defined $range_char;
|
|
$regexp .= ")";
|
|
}
|
|
}
|
|
return $regexp;
|
|
}
|
|
|
|
sub parseProsite {
|
|
local $_ = shift;
|
|
my $ac = $1 if /^AC (\w+)/m;
|
|
my ($id, $type) = ($1, $2) if /^ID (\w+); (\w+)/m;
|
|
my $de_line = join "\n", /^(DE.*\S)/mg;
|
|
my $de = join " ", $de_line =~ /^DE (.*\S)/mg;
|
|
my $pa = join "", /^PA (.*\S)/mg;
|
|
$pa =~ s/\.$//;
|
|
my $rule = join "\\\n", /^RU (.*\S)/mg;# obsolete
|
|
my $skip = /^CC \/SKIP-FLAG=TRUE/m;
|
|
my $pdoc = $1 if /^DO (\w+)/m;
|
|
my $tax = $1 if /^CC \/TAXO-RANGE=(.*?);/m;
|
|
my $rep = $1 if /^CC.*\/MAX-REPEAT=(\d+)/m;
|
|
my @sites;push (@sites, [$1, $2])
|
|
while(m/\/SITE=(\d+),(.*?);/g);
|
|
my @cutoffs = map {my %a=map {split /=/, $_, 2}
|
|
split /; */, $_; \%a } /^MA \/CUT_OFF: (.*)/mg;
|
|
my @pps;push (@pps, {'type'=>$1,'discriminators'=>[split(';\s*',$2)]})
|
|
while(m/^PP \/?(\w+):\s*([^\r\n]+)/mg);
|
|
return ($ac, $id, $type, $de, $pa, $rule,
|
|
$pdoc, $skip, $tax, $rep, \@sites, \@cutoffs, \@pps);
|
|
};
|
|
|
|
|
|
|
|
################################################################################
|
|
# integrated subs taken from FTRep.pm module:
|
|
|
|
use constant PI => 3.14159265;
|
|
use constant PI2 => 2*PI;
|
|
|
|
#-------------------------------------------------------------#
|
|
# sinus cardinal
|
|
sub sinc
|
|
{
|
|
my $x = shift;
|
|
return ($x == 0)? 1 : sin($x)/$x;
|
|
}
|
|
|
|
#-------------------------------------------------------------#
|
|
# Fourier transform
|
|
sub ft
|
|
{
|
|
my ($x,$score,$d,$size) = @_;
|
|
return 0 if ($size < 1);
|
|
my $t = PI2*($x-$d)/$size;
|
|
my $sint = sinc($t);
|
|
my $sint2 = $sint*$sint;
|
|
return $score * $sint2;
|
|
}
|
|
|
|
sub FT
|
|
{
|
|
my ($matchlist,$groups,$motif_param,$amplitude) = @_;
|
|
|
|
$amplitude = 0; # Turn it off
|
|
|
|
# ft and autocorrelation
|
|
my $ft = []; # for plot
|
|
my $fa = []; # for plot and integral
|
|
|
|
for (my $i = 0; $i < scalar(@$groups); $i++)
|
|
{
|
|
for (my $x = $groups->[$i]->{start}; $x <= $groups->[$i]->{end}; $x++) # walk around the group region in the sequence
|
|
{
|
|
my $y = 0; # store ft function result
|
|
my $yk = 0; # store ft with lag function
|
|
my $n_members = scalar(@{$groups->[$i]->{members}});
|
|
if ($n_members > 1)
|
|
{
|
|
for (my $g = 0 ; $g < $n_members; $g++)
|
|
{
|
|
my $k = ($g == 0)? 1 : -1; # deal with left border g-1,g becomes g,g+1
|
|
|
|
my $m = $groups->[$i]->{members}->[$g];
|
|
my $n = $groups->[$i]->{members}->[$g+$k];
|
|
my $m_amp = ($amplitude)? $amplitude : $matchlist->[$m]->{match_rscore}/$motif_param->{maxscore};
|
|
my $n_amp = ($amplitude)? $amplitude : $matchlist->[$n]->{match_rscore}/$motif_param->{maxscore};
|
|
|
|
$y += ft($x,$m_amp,$matchlist->[$m]->{match_pos},$matchlist->[$m]->{match_size});
|
|
$yk += ft($x+$motif_param->{motif_len}*$k,$n_amp,$matchlist->[$n]->{match_pos},$matchlist->[$n]->{match_size});
|
|
}
|
|
}
|
|
$ft->[$i][$x] += $y;
|
|
$fa->[$i][$x] += $y * $yk;
|
|
}
|
|
}
|
|
return ($ft,$fa);
|
|
}
|
|
|
|
#-------------------------------------------------------------#
|
|
# autocorrelation
|
|
sub AC
|
|
{
|
|
my ($matchlist,$groups,$fa) = @_;
|
|
my $ac = [];
|
|
|
|
# integral of autocorrelation for each group
|
|
for (my $i = 0; $i < scalar(@$groups); $i++)
|
|
{
|
|
$ac->[$i] = 0;
|
|
for (my $x = $groups->[$i]->{start}+1; $x <= $groups->[$i]->{end}; $x++)
|
|
{
|
|
$ac->[$i] += ($fa->[$i][$x] + $fa->[$i][$x-1]) / 2;
|
|
}
|
|
}
|
|
return $ac;
|
|
}
|
|
|
|
#-------------------------------------------------------------#
|
|
# select match
|
|
sub select_match
|
|
{
|
|
my ($matchlist,$groups,$ac,$cutoff,$filter) = @_;
|
|
my $max_prot_score = -999999;
|
|
|
|
for (my $i = 0; $i < scalar(@$groups); $i++)
|
|
{
|
|
foreach my $g (@{$groups->[$i]->{members}})
|
|
{
|
|
if ( $cutoff < $ac->[$i] )
|
|
{
|
|
$matchlist->[$g]->{selected} = 1;
|
|
}
|
|
|
|
if ($matchlist->[$g]->{match_nscore} > $max_prot_score)
|
|
{
|
|
$max_prot_score = $matchlist->[$g]->{match_nscore}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (defined($filter) && $max_prot_score < $filter)
|
|
{
|
|
for (my $i = 0; $i < scalar(@$groups); $i++)
|
|
{
|
|
foreach my $g (@{$groups->[$i]->{members}})
|
|
{
|
|
$matchlist->[$g]->{selected} = 0;
|
|
}
|
|
}
|
|
}
|
|
|
|
return (1);
|
|
}
|
|
|
|
#-------------------------------------------------------------#
|
|
# Sort positions
|
|
sub sortlist
|
|
{
|
|
my ($matchlist) = shift;
|
|
my %sortindex;
|
|
for (my $i = 0; $i < scalar(@$matchlist); $i++)
|
|
{
|
|
$sortindex{$matchlist->[$i]->{match_start}} = $i;
|
|
}
|
|
|
|
my @sortindex;
|
|
foreach (sort {$a<=>$b} keys %sortindex)
|
|
{
|
|
push(@sortindex,$sortindex{$_});
|
|
}
|
|
|
|
return @sortindex;
|
|
}
|
|
|
|
#-------------------------------------------------------------#
|
|
# Group matches
|
|
sub group_matches
|
|
{
|
|
my ($matchlist,$grouping_size) = @_;
|
|
|
|
# sort matches by their position
|
|
my @index = sortlist($matchlist);
|
|
|
|
my $keep_prot = 0; # used for filter
|
|
my ($start,$end);
|
|
my $group_index = 0;
|
|
my $groups = [];
|
|
|
|
# First match
|
|
push(@{$groups->[$group_index]->{members}},$index[0]);
|
|
$start = $matchlist->[$index[0]]->{match_start};
|
|
|
|
for (my $i = 1; $i < scalar(@index); $i++)
|
|
{
|
|
if ($matchlist->[$index[$i]]->{match_pos} - $matchlist->[$index[$i-1]]->{match_pos} > $grouping_size)
|
|
{
|
|
# Set limits for current group
|
|
$end = $matchlist->[$index[$i-1]]->{match_end};
|
|
$groups->[$group_index]{start} = $start;
|
|
$groups->[$group_index]{end} = $end;
|
|
$start = $matchlist->[$index[$i]]->{match_start};
|
|
|
|
# Create new group
|
|
$group_index++;
|
|
}
|
|
elsif (!$keep_prot)
|
|
{
|
|
$keep_prot = 1; # keep the protein if at least 2 match have distance < grouping_size
|
|
}
|
|
|
|
$groups->[$group_index]->{match_index} = 0;
|
|
push(@{$groups->[$group_index]->{members}},$index[$i]);
|
|
}
|
|
|
|
# Deal with last match
|
|
$end = $matchlist->[$index[scalar(@index)-1]]->{match_end};
|
|
$groups->[$group_index]->{start} = $start;
|
|
$groups->[$group_index]->{end} = $end;
|
|
|
|
if (!$keep_prot)
|
|
{
|
|
#$groups = []
|
|
}
|
|
return $groups;
|
|
}
|
|
|
|
|
|
|
|
|
|
################################################################################
|
|
# ps_scan_source.pl core:
|
|
|
|
|
|
# initializations & parameters processing
|
|
|
|
BEGIN {
|
|
$VERSION = '1.90';
|
|
}
|
|
|
|
# Can we use the IPC::Open2 module to communicate with
|
|
# pfscan via pipes (instead of temp files) ?S
|
|
eval { require IPC::Open2 };
|
|
my $NO_DIRECT_PIPE=$? || $^O eq "MSWin32" ? 1 : 0;
|
|
# change this to the absolute path to the programs,
|
|
# unless they are located in a directory in your PATH
|
|
# or use option --pfscan and/or --psa2msa to give the full path.
|
|
my $PFSCAN = 'pfscan';
|
|
my $PSA2MSA = 'psa2msa';
|
|
my $errcode = 0;
|
|
my $MOTIF_AC_REGEXP = '\w+\d+';#'PS\d{5}';
|
|
$|= 1;
|
|
use Getopt::Long;
|
|
use Data::Dumper;
|
|
use IO::File;
|
|
my @formats = qw(scan fasta psa msa gff pff epff sequence matchlist ipro);
|
|
my $formats_string = join " | ", @formats;
|
|
|
|
sub usage {
|
|
my $progname = $0;
|
|
$progname =~ s/.*[\\\/]//;
|
|
print <<EOF;
|
|
$progname [options] sequence-file(s)
|
|
ps_scan version $VERSION options:
|
|
-h : this help screen
|
|
|
|
Input/Output:
|
|
-e <string> : specify the ID or AC of an entry in sequence-file
|
|
-o <string> : specify an output format : scan | fasta | psa | msa |
|
|
gff | pff | epff | sequence | matchlist
|
|
-d <file> : specify a prosite motif file
|
|
-p <string> : specify a pattern or the AC of a prosite motif
|
|
-f <string> : specify a motif AC to scan against together with all its
|
|
related post-processing motifs (but show only specified
|
|
motif hits)
|
|
|
|
Selection:
|
|
-r : do not scan profiles
|
|
-m : only scan profiles
|
|
-s : skip frequently matching (unspecific) patterns and
|
|
profiles
|
|
-l <number> : cut-off level for profiles (default : 0)
|
|
|
|
Pattern option:
|
|
-x <number> : specify maximum number of accepted matches of X's in
|
|
sequence (default=0)
|
|
-g : Turn greediness off
|
|
-v : Turn overlaps off
|
|
-i : Allow included matches
|
|
-b <file> : use profiles from <f> to evaluate pattern matches
|
|
(a LevelTag is assigned to each pattern match). If no
|
|
file is specified, evaluator.dat will be used (searched
|
|
in the paths \$PROSITE/ and \$PROSITE/prosite/)
|
|
|
|
profile options:
|
|
-w pfsearch : Compares a query profile against a protein sequence
|
|
library
|
|
|
|
Other expert options:
|
|
--nopp : do not post-process matches
|
|
--reverse : randomize the sequence database by taking the
|
|
reverse sequence of each individual entry
|
|
--shuffle : randomize the sequence database by local shuffling
|
|
of the residues in windows of 20 residues
|
|
--gff <file> : use an existing ps_scan gff result <file> as input to
|
|
e.g. post-process it and/or to reformat it into another
|
|
format (defined by -o <string>), instead of doing a
|
|
real scan.
|
|
-w pfsearch -R : use raw scores rather than normalized scores for
|
|
match selection
|
|
-w pfsearch -C <x> : report only match scores higher than the specified
|
|
parameter <x>; an integer argument is interpreted
|
|
as a raw score value, a decimal argument as a
|
|
normalized score value. -R and -C options can be
|
|
combined.
|
|
--pfscan <path> : pathname to pfscan executable (if not defined will be
|
|
'pfscan', so executable has to be found within PATH env).
|
|
This option could be used e.g. to use pftool v3 pfscanV3.
|
|
|
|
Note:
|
|
The sequence-file may be in Swiss-Prot or FASTA format.
|
|
|
|
If no prosite motif file is specified, prosite.dat will be used
|
|
(searched in the paths \$PROSITE/ and \$SPROT/prosite/).
|
|
|
|
There may be several -d, -p and -e arguments.
|
|
EOF
|
|
exit 1;
|
|
}
|
|
|
|
my $opt_noprofiles;
|
|
my $opt_onlyprofiles;
|
|
my $opt_skip;
|
|
my $opt_skiponly;
|
|
my $opt_help;
|
|
my $opt_format;
|
|
my $opt_max_x;
|
|
my $opt_nongreedy;
|
|
my $opt_nooverlaps;
|
|
my $opt_miniprofiles;
|
|
my $opt_includes;
|
|
my $opt_level = 0;
|
|
my $opt_rep_pp_4allprofiles;
|
|
my $opt_pfsearch;
|
|
my $opt_cutoff;
|
|
my $opt_raw;
|
|
my $opt_minhits;
|
|
my $opt_maxhits;
|
|
my $opt_filterheader;
|
|
my $opt_reverse;
|
|
my $opt_shuffle;
|
|
my $opt_no_postprocessing;
|
|
|
|
|
|
# list of prosite 'dat' files to be scanned
|
|
my @prosite_files;
|
|
# !can contain a list of user pattern (and/or)
|
|
# list of prosite AC (profile/pattern)!
|
|
my @motifAC_or_userpattern;
|
|
# list of ID/AC
|
|
my @entries;
|
|
# list of motif AC for which one wants to retrieve (and scan)
|
|
# all other motif linked by post-processing
|
|
my @followpp;
|
|
# list of ps_scan gff result files to be scanned
|
|
my @external_gff_files;
|
|
|
|
my $SLASH = $^O eq "MSWin32" ? "\\" : "\/";
|
|
my $TMPDIR = ".";
|
|
for my $dir ( $ENV{TMPDIR}, $ENV{SP_TEMP}, $ENV{TMP}, $ENV{TEMP},
|
|
"/tmp", "c:\\temp", "c:\\tmp" ) {
|
|
if ( defined($dir) and -d $dir ) {
|
|
$TMPDIR = $dir;
|
|
last;
|
|
}
|
|
}
|
|
my $TMP_COUNTER = 1;
|
|
my $TMP_BASE = int( rand( 1000000 ) ) + int( substr( abs( $$ ), -6 ) );
|
|
|
|
my $scan_profiles;
|
|
my $scan_pattern;
|
|
|
|
my $last_profile_tmp_filename;
|
|
|
|
Getopt::Long::Configure ("bundling", "no_ignorecase");
|
|
GetOptions (
|
|
"r" => \$opt_noprofiles,
|
|
"m" => \$opt_onlyprofiles,
|
|
"s" => \$opt_skip,
|
|
"h" => \$opt_help,
|
|
"v" => \$opt_nooverlaps,
|
|
"i" => \$opt_includes,
|
|
"g" => \$opt_nongreedy,
|
|
"b:s" => \$opt_miniprofiles,
|
|
"x=i" => \$opt_max_x,
|
|
"l=i" => \$opt_level,
|
|
"o=s" => \$opt_format,
|
|
"d=s" => \@prosite_files,
|
|
"p=s" => \@motifAC_or_userpattern,
|
|
"e=s" => \@entries,
|
|
"f=s" => \@followpp,
|
|
# a: undocumented, just for PROSITE team to test repeat post-processing
|
|
# on any normal ('!' level0) profile...
|
|
"a" => \$opt_rep_pp_4allprofiles,
|
|
"w=s" => \$opt_pfsearch,
|
|
"C=f" => \$opt_cutoff,
|
|
"R" => \$opt_raw,
|
|
"pfscan=s" => \$PFSCAN,
|
|
"psa2msa=s" => \$PSA2MSA,
|
|
"minhits=i" => \$opt_minhits,
|
|
"maxhits=i" => \$opt_maxhits,
|
|
"filterheader=s" => \$opt_filterheader,
|
|
"reverse" => \$opt_reverse,
|
|
"shuffle=i" => \$opt_shuffle,
|
|
"skipflag-only" => \$opt_skiponly,
|
|
"nopp" => \$opt_no_postprocessing,
|
|
"gff=s" => \@external_gff_files
|
|
) or &usage;
|
|
|
|
&usage if $opt_help;
|
|
&usage if !@ARGV && -t STDIN && !@external_gff_files;
|
|
if ( $opt_pfsearch ) {
|
|
die "OPTION CONFLICT: can't use option".
|
|
" -reverse together with -w option (pfsearch)"
|
|
if $opt_reverse;
|
|
die "OPTION CONFLICT: can't use option".
|
|
" -shuffle together with -w option (pfsearch)"
|
|
if $opt_shuffle;
|
|
die "OPTION CONFLICT: can't use option".
|
|
" -e together with -w option (pfsearch)"
|
|
if @entries;
|
|
&usage if !@prosite_files &&
|
|
!grep {/^$MOTIF_AC_REGEXP$/} @motifAC_or_userpattern;
|
|
# integer cutoff forces raw scores
|
|
$opt_raw = $1 if defined($opt_cutoff) and $opt_cutoff=~ /^(\d+)$/mg;
|
|
}
|
|
my $use_pfsearchV3 = index( $opt_pfsearch, 'pfsearchV3') != -1 ? 1 : 0;
|
|
|
|
|
|
my $scan_behavior;
|
|
$scan_behavior |= 1 if $opt_nooverlaps;
|
|
$scan_behavior |= 2 if $opt_includes;
|
|
|
|
$opt_format = "scan" unless defined $opt_format;
|
|
$opt_format =~ tr/A-Z/a-z/;
|
|
die "ERROR:Output format must be one of $formats_string\n"
|
|
unless grep {$_ eq $opt_format} @formats;
|
|
my $opt_psa_or_msa = $opt_format eq "msa" || $opt_format eq "psa";
|
|
|
|
# user patterns specified with -p option
|
|
my @userpat;
|
|
# hash of specified (-p option) prosite ac (without specified user patterns)
|
|
my $specifiedPrositeMotifByAc={};
|
|
map {
|
|
( /^$MOTIF_AC_REGEXP$/ ?
|
|
$specifiedPrositeMotifByAc->{$_} = 1 : push @userpat,$_ )
|
|
} @motifAC_or_userpattern;
|
|
# if no prosite ac specified: undef struct...
|
|
keys(%$specifiedPrositeMotifByAc) or $specifiedPrositeMotifByAc = undef;
|
|
|
|
# find default prosite.dat file
|
|
if ( !@prosite_files && !@userpat ) {
|
|
if (defined $ENV{PROSITE} and -e "$ENV{PROSITE}/prosite.dat") {
|
|
@prosite_files = "$ENV{PROSITE}/prosite.dat";
|
|
}
|
|
elsif ( defined $ENV{PROSITE} and -e $ENV{PROSITE} ) {
|
|
@prosite_files = $ENV{PROSITE};
|
|
}
|
|
elsif ( defined $ENV{SPROT} and -e "$ENV{SPROT}/prosite/prosite.dat" ) {
|
|
@prosite_files = "$ENV{SPROT}/prosite/prosite.dat";
|
|
}
|
|
elsif ( -e "prosite.dat" ) {
|
|
@prosite_files = "prosite.dat";
|
|
}
|
|
else {
|
|
die "prosite.dat file not found, please use the -d option";
|
|
}
|
|
}
|
|
|
|
# -b option with no pathname specified: find default evaluator.dat
|
|
if ( defined($opt_miniprofiles) && !$opt_miniprofiles ) {
|
|
if ( defined $ENV{PROSITE} and -e "$ENV{PROSITE}/evaluator.dat" ) {
|
|
$opt_miniprofiles = "$ENV{PROSITE}/evaluator.dat";
|
|
}
|
|
elsif ( defined $ENV{SPROT} and -e "$ENV{SPROT}/prosite/evaluator.dat" ) {
|
|
$opt_miniprofiles = "$ENV{SPROT}/prosite/evaluator.dat";
|
|
}
|
|
elsif ( -e "evaluator.dat" ) {
|
|
$opt_miniprofiles = "evaluator.dat";
|
|
}
|
|
else {
|
|
die "evaluator.dat file not found,".
|
|
" please specify which evaluator.dat to use after the -b option";
|
|
}
|
|
}
|
|
|
|
my %SkipFlag;
|
|
my %KnownFalsePos;
|
|
my @MotifInfo;
|
|
|
|
my $hideMotifByPSAC = {};
|
|
my $files_miniprofiles = {};
|
|
my $postProcessingByPSAC = {};# PP data by PSAC (then by PP type)
|
|
my $motifRank4PostProcessing = {};
|
|
|
|
# dispatch table for post processing of (matching) target
|
|
# dataset against effector dataset (within hits on one sequence)
|
|
my $postProcessDispatchTable = {
|
|
# match data =
|
|
# 0 1 2 3 4(if profile, sequence if 1st match on a pattern) 5 6 7 8 9 10
|
|
# [matchseq, from, to, pfac, pffrom (or full entry seq if pattern + option b), pfto, rawscore, nscore, level, leveltag, de
|
|
|
|
# promote weak match (add 1 to level)
|
|
# if effector has at least one hit at level>=$minlevel or 0
|
|
# $target = ref to array that contains match data
|
|
# for target motif on 1 sequence
|
|
# $effector = ref to array that contains match data
|
|
# for effector motif on 1 sequence
|
|
'PROMOTED_BY' => sub {
|
|
my ( $target, $effector, $minlevel ) = @_;
|
|
return unless $target and $effector;
|
|
|
|
$minlevel ||= 0;
|
|
( my $ac = $effector->[0]->[3] || '' ) =~ s/\|\w+$//;
|
|
map { $_->[8]++, $_->[9] = "PROMOTED_BY_$ac" if $_->[8] < 0 } @$target
|
|
if grep { defined( $_->[8] ) && $_->[8] >= $minlevel } @$effector;
|
|
},
|
|
|
|
# gives a "level" (in fact only shown in LevelTag; displayed level will
|
|
# still be empty like for any pattern) to pattern matches:
|
|
# if the effector detects the match at level 0, the match is given level 0,
|
|
# LevelTag "(0)", else LevelTag "(-1)"
|
|
# NOTE: effector hits are calculated here!
|
|
'EVALUATED_BY' => sub {
|
|
my ( $target, $effector, $mp_ac ) = @_;
|
|
# note $effector is empty as matches will be calculated here!
|
|
return unless $target and $mp_ac;
|
|
|
|
my $sequence = $target->[0]->[4] or return;
|
|
|
|
my $mini_tmpfile;
|
|
my $mini_exist;
|
|
my $miniprofile_dat = $opt_miniprofiles;
|
|
if ( exists $files_miniprofiles->{ $mp_ac } ) {
|
|
$mini_tmpfile = $files_miniprofiles->{ $mp_ac };
|
|
$mini_exist = "true";
|
|
}
|
|
else {# each new minifile is put in $files_miniprofiles
|
|
$mini_tmpfile = tmpnam();
|
|
open ( IN, $miniprofile_dat ) or
|
|
die "can't open evaluator.dat file [$miniprofile_dat]";
|
|
$/ = "\n//";# entry separator
|
|
# note full sep. is \n//\n but in DOS format is \r\n//\r\n, so
|
|
# only \n// matches both
|
|
while ( <IN> ) {
|
|
( my $entry = $_ ) =~ s/^\s+//;
|
|
if ( $entry =~ /$mp_ac;/ ) {
|
|
$mini_exist = 1;
|
|
open (OUT,">$mini_tmpfile");
|
|
print OUT $entry."\n";
|
|
close OUT;
|
|
last;
|
|
}
|
|
}
|
|
close IN;
|
|
|
|
$files_miniprofiles->{ $mp_ac } = $mini_tmpfile if -e "$mini_tmpfile";
|
|
}
|
|
$/ = "\n";
|
|
# if no miniprofile correspond to pattern then no evaluation
|
|
return unless $mini_exist;
|
|
|
|
my $struct = do_profile_scan( $mini_tmpfile, undef, $sequence );
|
|
# returned: hash ref where key = mini profile AC or AC|ID!
|
|
# take first (and only) value
|
|
$effector = ( values %$struct )[0];
|
|
foreach my $target_hit ( @$target ) {
|
|
my $level = -1;
|
|
foreach my $effector_hit ( @$effector ) {
|
|
last if $effector_hit->[1] > $target_hit->[2];
|
|
$level = 0, last if $effector_hit->[1] <= $target_hit->[1]
|
|
and $effector_hit->[2] >= $target_hit->[2];
|
|
}
|
|
# $target_hit->[8] = $level; do not change level, just tag!
|
|
$target_hit->[9] = "($level)";
|
|
}
|
|
},
|
|
|
|
# downgrade good match (substract 1 to level)
|
|
# if effector has at least one hit at level>=$minlevel or 0
|
|
'DEMOTED_BY' => sub {
|
|
my ( $target, $effector, $minlevel ) = @_;
|
|
return unless $target and $effector;
|
|
|
|
$minlevel ||= 0;
|
|
map { $_->[8]--, $_->[9] = 'DEMOTED' if $_->[8] >= 0 } @$target
|
|
if grep { defined( $_->[8] ) && $_->[8] >= $minlevel } @$effector;
|
|
},
|
|
|
|
# downgrade target/effector overlapping matches if their score is
|
|
# smaller than competitor score
|
|
'COMPETES_HIT_WITH' => sub {
|
|
my ( $hit_set_a,$hit_set_b,$overlap ) = @_;
|
|
return unless $hit_set_a and $hit_set_b;
|
|
|
|
my @set_a = sort { $a->[1] <=> $b->[1] || $a->[2] <=> $b->[2] }
|
|
#grep {$_->[1]||=0;$_->[2]||=0;$_->[8]>=0}
|
|
@$hit_set_a;
|
|
# create array of () hits sorted by
|
|
# position in set a
|
|
my @set_b = sort { $a->[1] <=> $b->[1] || $a->[2] <=> $b->[2] }
|
|
#grep {$_->[1]||=0;$_->[2]||=0;$_->[8]>=0}
|
|
@$hit_set_b;
|
|
# create array of () hits sorted by
|
|
# position in set b
|
|
return unless @set_a and @set_b;
|
|
( my $ac_a = $hit_set_a->[0]->[3] || '' ) =~ s/\|\w+$//;
|
|
( my $ac_b = $hit_set_b->[0]->[3] || '' ) =~ s/\|\w+$//;
|
|
$overlap=0 if !$overlap || $overlap!~m/^\d+$/;
|
|
# allowed overlap size (only consider overlap matches if
|
|
# overlap size > this value)
|
|
foreach my $hit_a ( @set_a ) {
|
|
# for each (sorted) level>=0 hits in set a
|
|
my $a_start = $hit_a->[1] + $overlap;
|
|
my $a_stop = $hit_a->[2] - $overlap;
|
|
next if $a_start > $a_stop;
|
|
foreach my $hit_b ( @set_b ) {
|
|
# for each (sorted) level>=0 hits in set b
|
|
next if $a_start > $hit_b->[2];
|
|
# hit b is before examined hit a: next
|
|
last if $a_stop < $hit_b->[1];
|
|
# hit b is after examined hit a: last
|
|
|
|
# overlap: demote hit with lowest score
|
|
if ( $hit_b->[7] < $hit_a->[7] ) {
|
|
$hit_b->[8]--;
|
|
$hit_b->[9] = "OUTCOMPETED_HIT_BY_$ac_a";
|
|
}
|
|
elsif ( $hit_a->[7] < $hit_b->[7]) {
|
|
$hit_a->[8]--;
|
|
$hit_a->[9] = "OUTCOMPETED_HIT_BY_$ac_b";
|
|
}
|
|
}
|
|
}
|
|
},
|
|
|
|
# downgrade target/effector matches from set with the lowest best score
|
|
'COMPETES_SEQ_WITH' => sub {
|
|
my ( $hit_set_a,$hit_set_b ) = @_;
|
|
return unless $hit_set_a and $hit_set_b;
|
|
( my $ac_a = $hit_set_a->[0]->[3] || '' ) =~ s/\|\w+$//;
|
|
( my $ac_b = $hit_set_b->[0]->[3] || '' ) =~ s/\|\w+$//;
|
|
|
|
my $max_a = 0.0; map { $max_a=$_->[7] if $_->[7]>$max_a } @$hit_set_a;
|
|
my $max_b = 0.0; map { $max_b=$_->[7] if $_->[7]>$max_b } @$hit_set_b;
|
|
if ( $max_a>$max_b ) {
|
|
map {
|
|
$_->[8]--;
|
|
$_->[9] = "OUTCOMPETED_SEQ_BY_$ac_a"
|
|
} @$hit_set_b
|
|
}
|
|
else {
|
|
map {
|
|
$_->[8]--;
|
|
$_->[9] = "OUTCOMPETED_SEQ_BY_$ac_b"
|
|
} @$hit_set_a
|
|
}
|
|
},
|
|
|
|
# promote weak match (add 1 to level)
|
|
# if effector has at least one hit with a (n)score >=$minscore
|
|
# and there are more than one target match
|
|
# Note: the target is usually the effector. Method used for repeats
|
|
'PROMOTE_REPEAT' => sub {
|
|
my ( $target, $effector, $cutoff ) = @_;
|
|
return unless $target and $effector;
|
|
|
|
$cutoff ||= 0;
|
|
#(my $ac = $effector->[0]->[3] || '') =~ s/\|\w+$//;
|
|
map { $_->[8]++, $_->[9] = "PROMOTED_REPEAT" if $_->[8] < 0 } @$target
|
|
if ( @$target >1 and
|
|
grep {
|
|
defined( $_->[7] ) && $_->[7] >= $cutoff
|
|
} @$effector
|
|
);
|
|
},
|
|
|
|
# promote weak match (add 1 to level), if there is at least one hit at
|
|
# level>=0 (RDM1)
|
|
'PROMOTER_RDM1' => sub {
|
|
my $target = shift or return;
|
|
# ('R?' in profile MA/CUT_OFF line)
|
|
map { $_->[8]++, $_->[9] = "R" if $_->[8] < 0 } @$target
|
|
if grep {defined($_->[8]) && $_->[8] >= 0} @$target;
|
|
},
|
|
|
|
# promote weak match (add 1 to level), if the sum of their (n)score
|
|
# is >= threshold OR one hit at level>=0 ('RR' in profile MA/CUT_OFF line)
|
|
'PROMOTER_RDM1_OR_RDM2' => sub {
|
|
my ( $target, $effector, $threshold ) = @_;
|
|
return unless $target and $threshold;
|
|
|
|
# RDM1 (>=1 level 0 match -> promote level<0)
|
|
if ( grep { defined( $_->[8] ) && $_->[8] >= 0 } @$target ) {
|
|
map { $_->[8]++, $_->[9] = "R" if $_->[8] < 0 } @$target;
|
|
}
|
|
# RDM2 (when no hit with level>=0) promote level<0 if hit score
|
|
# sum >= treshhold
|
|
else {
|
|
my $sum = 0.0; map { $sum+=$_->[7] || 0; } @$target;
|
|
map { $_->[8]++, $_->[9] = "r" if $_->[8] < 0 } @$target
|
|
if $sum >= $threshold;
|
|
}
|
|
},
|
|
|
|
# same as PROMOTE_RDM1_OR_RDM2 (but on repeats from '!' profile,
|
|
# only with -a option)
|
|
'PROMOTER_RDM1_OR_RDM2_4!REP' => sub {
|
|
my ( $target, $effector, $threshold ) = @_;
|
|
return unless $target and $threshold;
|
|
|
|
# RDM1 (>=1 level 0 match -> promote level<0)
|
|
if ( grep { defined($_->[8]) && $_->[8] >= 0 } @$target ) {
|
|
map { $_->[8]++, $_->[9] = "?R" if $_->[8] < 0 } @$target;
|
|
}
|
|
# RDM2 (when no hit with level >= 0) promote level < 0 if hit score
|
|
# sum >= treshhold
|
|
else {
|
|
my $sum = 0.0; map { $sum += $_->[7] || 0; } @$target;
|
|
map { $_->[8]++, $_->[9] = "?r" if $_->[8] < 0 } @$target
|
|
if $sum >= $threshold;
|
|
}
|
|
},
|
|
|
|
'FTREP' => sub {
|
|
my ( $target, $effector, $paramstr ) = @_;
|
|
return unless $target;
|
|
|
|
my ( $max_score, $ac_cutoff, $size, $gf, $filter ) =
|
|
split( '\|', $paramstr );
|
|
|
|
my $motif_param = {
|
|
'motif_id' => 'fake',
|
|
'motif_acc' => 'fake',
|
|
'motif_len' => $size,
|
|
'maxscore' => $max_score,
|
|
'filter' => undef,
|
|
'normscore1' => undef,
|
|
};
|
|
|
|
my $matchlist = [ map {
|
|
my $estart = $_->[1] - ($_->[4]||1) +1;
|
|
my $estop = $_->[2] - ($_->[5]||-1) -1;
|
|
my $size = 1 + $estop - $estart;
|
|
{
|
|
'match_start' => $_->[1],
|
|
'match_end' => $_->[2],
|
|
'match_nscore' => $_->[7],
|
|
'match_rscore' => $_->[6],
|
|
'match_pos' => $size/2.0 + $estart,
|
|
'match_size' => $size
|
|
}
|
|
} @$target ];
|
|
|
|
# Group matches
|
|
my $groups = group_matches( $matchlist, $gf * $size );
|
|
return unless @$groups;
|
|
|
|
# Fourier
|
|
my ( $ft, $fa ) = FT( $matchlist, $groups, $motif_param, 1 );
|
|
|
|
# Autocorrelation
|
|
my $ac = AC( $matchlist, $groups, $fa );
|
|
|
|
# Get evalues and select matches
|
|
select_match( $matchlist, $groups, $ac, $ac_cutoff, $filter );
|
|
|
|
for ( my $i = 0 ; $i < @$target ; $i++ ) {
|
|
my $match = $target->[ $i ];
|
|
my $ftrepmatch = $matchlist->[ $i ];
|
|
if ( $match->[ 8 ] >= 0 ) {
|
|
# good match: leave it as-is
|
|
}
|
|
elsif ( $ftrepmatch->{ selected } ) {# promote (to level 0)
|
|
$match->[ 8 ] = 0;
|
|
$match->[ 9 ] = 'PROMOTED_FTREP'
|
|
}
|
|
else {# demote (-1 to -2)
|
|
$match->[ 8 ] = -2;
|
|
$match->[ 9 ] = 'DEMOTED_FTREP'
|
|
}
|
|
}
|
|
}
|
|
|
|
};
|
|
|
|
my $allowBidirectionalPP = {
|
|
'PROMOTED_BY' => 1,
|
|
# note: COMPETES_... are 'bidirectional' inside their
|
|
# $postProcessDispatchTable sub but we don't wan't them to be called
|
|
# bidirectionally (in pp_scan)...
|
|
};
|
|
|
|
|
|
main();
|
|
exit $errcode;
|
|
|
|
|
|
##############################################################################
|
|
# Methods
|
|
|
|
sub tmpnam {
|
|
my $tmp;
|
|
do {
|
|
$tmp = $TMPDIR.$SLASH."ps${TMP_BASE}-".$TMP_COUNTER++.".tmp";
|
|
} while (-e $tmp);
|
|
return $tmp;
|
|
}
|
|
|
|
# -------------------------- output methods --------------------------
|
|
|
|
# format a field with a certain width
|
|
sub pf {
|
|
return $_[0] . ( " " x ( $_[1] - length $_[0] ) )
|
|
}
|
|
|
|
my $HIT_COUNT = 0;
|
|
sub dispHits {# display hits (for 1 sequence, 1 motif)
|
|
my ( $header, $sq, $hits, $seqid, $de, $aclist, $psac, $psid, $psde ) = @_;
|
|
return if !$hits || !$psac || $hideMotifByPSAC->{ $psac };
|
|
|
|
$sq ||= ''; $seqid ||= ''; $de ||= ''; $psid ||= ''; $psde ||= '';
|
|
( $de = $hits->[0]->[10] || '' ) =~ s/\.\s*$// if !$de && $opt_pfsearch;
|
|
$de =~ s/[\n\r\t]/ /g;
|
|
|
|
# Hit array:[matchseq, from, to, ac, pffrom, pfto, rawscore, nscore, level, leveltag, de, [repstuff!?]]]
|
|
# 0 1 2 3 4 5 6 7 8 9 10 11
|
|
# Note: might ~vary depending if pattern/profile pfscan, pfsearch!
|
|
|
|
# remove unwanted matches with level < opt_level
|
|
my $visible_hits;
|
|
@$visible_hits =
|
|
grep {
|
|
# if nscore [7] is undef = match is from a pattern : show match
|
|
# (even if has level < 0 as now pattern might (~canceled) have
|
|
# level if used with -b miniprofile option!) or if match
|
|
# level [8] >= user requested minimal level
|
|
!defined( $_->[7] ) || $_->[8] >= $opt_level
|
|
} @$hits;
|
|
@$visible_hits or return;
|
|
|
|
unroll_hits( $visible_hits );
|
|
if ( defined $opt_maxhits ) {
|
|
splice( @$visible_hits, $opt_maxhits - $HIT_COUNT );
|
|
exit 0 if $HIT_COUNT >= $opt_maxhits ;
|
|
}
|
|
my $hit_count = @$visible_hits;
|
|
$HIT_COUNT += $hit_count;
|
|
return unless $hit_count;
|
|
return if defined ( $opt_minhits ) and $opt_minhits > $hit_count;
|
|
|
|
$visible_hits->[0]->[4] = undef unless defined( $visible_hits->[0]->[7] );
|
|
# clear sequence put in [4] in 1st pattern hit!
|
|
|
|
print $header if defined $header;
|
|
|
|
if ( $opt_format eq "ipro" ) {
|
|
# fake & hidden output: gff with additional post-processing for patterns!
|
|
# (By pattern) if one pattern hit (on this entry) has been promoted by
|
|
# a corresponding mini profile (EVALUATED_BY PP) to "(0)", promote all
|
|
# other matches (of the same pattern, in the same entry)
|
|
# (n.b. only patterns have '(<level>)' leveltags)
|
|
# FIXME: why not make this a standard (& public) PP behavior
|
|
# (e.g. in EVALUATED_BY PP)
|
|
my $ok_by_ac = {};
|
|
map{ $ok_by_ac->{ $_->[3] } = 1 if $_->[9] eq '(0)' } @$visible_hits;
|
|
@$visible_hits =
|
|
map { $_->[9] = '(0)' if $ok_by_ac->{ $_->[3] };
|
|
$_ } @$visible_hits;
|
|
# n.b. pattern leveltags are only (-1) or (0) (if not, this
|
|
# will give unexpected results!)
|
|
$opt_format = 'gff';
|
|
}
|
|
|
|
if ( $opt_format eq "fasta" or $opt_psa_or_msa ) {
|
|
for my $hit ( @$visible_hits ) {
|
|
my ( $subseq, $from, $to, $_psac, $pffrom, $pfto, $rawscore,
|
|
$nscore, $leveln, $levelt, $seqde ) = @$hit;
|
|
my $print_level = defined( $levelt ) ? "L=$levelt " :
|
|
defined($leveln) ? "L=$leveln " : "";
|
|
if ( $opt_pfsearch ) {
|
|
$opt_psa_or_msa = 1;
|
|
print ">$seqid/$from-$to : $de : $psac $print_level\n";
|
|
}
|
|
else {
|
|
print ">$seqid/$from-$to : $psac $psid $print_level\n";
|
|
}
|
|
if ( $subseq and $opt_psa_or_msa ) {# pfscan output
|
|
$subseq =~ s/\n?$/\n/; # add \n to scanPattern output
|
|
print $subseq;
|
|
}
|
|
else {
|
|
$subseq = substr( $sq, $from-1, $to-$from+1 );
|
|
while ( $subseq =~ /(.{1,60})/g ) {
|
|
print "$1\n";
|
|
}
|
|
print "\n" if $subseq eq "";
|
|
}
|
|
}
|
|
}
|
|
elsif ( $opt_format eq "pff" || $opt_format eq "epff" ) {
|
|
for my $hit ( @$visible_hits ) {
|
|
my @pff = @$hit;
|
|
$pff[8] = $pff[9] if defined $pff[9];
|
|
# move LevelTag into Level if defined
|
|
pop @pff while @pff>9; # remove fields beyond numeric level
|
|
my $subseq = shift @pff; # remove subseq
|
|
$pff[2] = $psac || $pff[2];
|
|
if ($opt_format eq "epff") {
|
|
push @pff, "" while @pff < 8;
|
|
$subseq =~ s/\s//g;
|
|
push @pff, $subseq;
|
|
}
|
|
print $seqid, "\t", join( "\t", @pff ), "\n";
|
|
}
|
|
}
|
|
elsif ( $opt_format eq "gff" ) {
|
|
for my $hit ( @$visible_hits ) {
|
|
my ( $subseq, $from, $to, $_psac, $pffrom, $pfto, $rawscore,
|
|
$nscore, $leveln, $levelt, $seqde ) = @$hit;
|
|
print join( "\t", $seqid, "ps_scan|v$VERSION", $psac,
|
|
$from, $to, $nscore || ".", ".", "." );
|
|
my @attr;
|
|
if ( $psid ) {
|
|
$psid =~ s/.*\|//; # FIXME!: why!?
|
|
push @attr, "Name \"$psid\"";
|
|
}
|
|
push @attr, "AccessionNumbers " .
|
|
join( " ", map {"\"$_\""} @$aclist )
|
|
if defined $aclist and @$aclist;
|
|
push @attr, "Level $leveln" if defined $leveln;
|
|
push @attr, "LevelTag \"$levelt\"" if defined $levelt;
|
|
push @attr, "RawScore $rawscore" if defined $rawscore;
|
|
push @attr, "FeatureFrom $pffrom"
|
|
if defined $pffrom and $pffrom =~ /\d+/;
|
|
push @attr, "FeatureTo $pfto" if defined $pfto;
|
|
$subseq =~ s/\s//g;
|
|
push @attr, "Sequence \"$subseq\"" if defined $subseq;
|
|
push @attr, "SequenceDescription \"$de\"" if $de;
|
|
push @attr, "SkipFlag 1" if $SkipFlag{$psac};
|
|
push @attr, "KnownFalsePos $KnownFalsePos{$psac}"
|
|
if exists $KnownFalsePos{ $psac };
|
|
print "\t", join " ; ", @attr if @attr;
|
|
print "\n";
|
|
}
|
|
}
|
|
elsif ( $opt_format eq "sequence" ) {
|
|
print ">$seqid : $de : $psac\n", map { "$_\n" }
|
|
$sq =~ /(.{1,60})/g if @$visible_hits;
|
|
}
|
|
else {# default output format
|
|
# fasta like header
|
|
if ( $opt_pfsearch ) {
|
|
print ">$seqid : $de : $psac\n";
|
|
}
|
|
else {
|
|
print ">$seqid : $psac $psid $psde\n";
|
|
}
|
|
# hits
|
|
for my $hit ( @$visible_hits ) {
|
|
my ( $subseq, $from, $to, $_psac, $pffrom, $pfto, $rawscore,
|
|
$nscore, $leveln, $levelt, $seqde ) = @$hit;
|
|
my $print_level = defined( $levelt ) ? " L=$levelt" :
|
|
defined( $leveln ) ? " L=$leveln" : "";
|
|
my $fromto = "$from - $to";
|
|
print " " x ( 13-length $fromto ), $fromto;
|
|
if ( $subseq ) {# pfscan output
|
|
$subseq =~ s/\n?$/\n/;# add \n to scanPattern output
|
|
$subseq =~ s/^(?<!\A)(.*)/ $1/mg;
|
|
$subseq =~
|
|
s/^(.*)/$1 . (" " x (60-length $1)) . $print_level/e
|
|
if $print_level;
|
|
print " ", $subseq;
|
|
}
|
|
else {
|
|
$subseq = substr( $sq, $from - 1, $to-$from + 1 );
|
|
my $notfirst;
|
|
while ( $subseq =~ /(.{1,60})/g ) {
|
|
print " " x 13 if $notfirst++;
|
|
print " $1\n";
|
|
}
|
|
print "\n" if $subseq eq "";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
sub showMatchList {# show (protein sequence) match list of (selected) motifs
|
|
for my $pattern ( @MotifInfo ) {
|
|
my ( $psac, $psid, $type, $psde, $pat, $skip ) = @$pattern;
|
|
next if $type eq "MATRIX";# not available for profiles
|
|
my $n_hits = 0;
|
|
my $n_match = 0;
|
|
print "$psid\n";
|
|
for my $sprotfile ( @ARGV ) {
|
|
open SPROTFILE, $sprotfile or die "Cannot open $sprotfile : $!";
|
|
my $entry = "";
|
|
# Note: only works on Swiss-Prot flat file format not on fasta ...
|
|
while ( <SPROTFILE> ) {
|
|
$entry .= $_;
|
|
if ( /^\/\// ) {
|
|
if ( $entry =~
|
|
/^\s*ID\s+(\w+)/m ) {
|
|
my $id = $1;
|
|
if ( $entry =~ /^\s*AC\s+(\w+)/m ) {
|
|
my $ac = $1;
|
|
my @de = /^\s*DE\s+(.+)/mg;
|
|
my $de;
|
|
my $add_space = 0;
|
|
for ( @de ) {
|
|
$de .= " " if $add_space;
|
|
$de .= $_;
|
|
$add_space = !/-$/;
|
|
}
|
|
if ( $entry =~
|
|
/^\s*SQ\s+SEQUENCE\b.*\n((.+\n)+)/m ) {
|
|
my $sq = $1;
|
|
$sq =~ tr/A-Z//cd;
|
|
$de = "$id $de" if $id and $de;
|
|
my $hits;
|
|
$hits = scanPatternWrapper([$pat, $sq,
|
|
$scan_behavior, $opt_max_x,
|
|
$opt_miniprofiles], $psid);
|
|
if ($hits and @$hits) {
|
|
print pf($ac, 6), ", ", pf($id, 10),
|
|
", T; ", scalar(@$hits), "\n";
|
|
$n_hits += @$hits;
|
|
$n_match++;
|
|
}
|
|
}
|
|
else {
|
|
warn "No sequence found in entry $id";
|
|
}
|
|
}
|
|
}
|
|
$entry = '';
|
|
}
|
|
}
|
|
}
|
|
print " ", pf( $n_match, 13 ), " ", $n_hits, "\n";
|
|
}
|
|
}
|
|
|
|
|
|
sub process_external_gff {
|
|
# process existing external gff ps_scan results, post-process (unless --nopp)
|
|
# and/or redisplay them in any available output format
|
|
# (usefull when scanning profile by profile on a cluster = no correct pp
|
|
# could be performed, then reuse ps_scan with --gff on collected results
|
|
# where correct pp (as all the distinct motif matches are combined) could be
|
|
# now performed. Caution: memory usage!!...
|
|
# Nov25 2010
|
|
|
|
my $psinfo_bypsac = {};
|
|
map{# store motif info by motif ac
|
|
$_->[0] and $psinfo_bypsac->{ $_->[0] } = $_;
|
|
} @MotifInfo;
|
|
|
|
# parse input gff
|
|
my $hits_by_psac_by_seq = {};
|
|
my $info_by_seqid = {};
|
|
foreach my $gff_file ( @external_gff_files ) {
|
|
next unless -e $gff_file;
|
|
open( FH, $gff_file ) or die "can't open [$gff_file] file";
|
|
while( <FH> ) {# read gff line by line
|
|
s/\r?\n$//;
|
|
# parse core elements + extra (attributes)
|
|
my ( $id, $nil1, $psac, $from, $to, $nscore,
|
|
$nil2, $nil3, $extra ) = split /\t/;
|
|
$extra ||= '';
|
|
|
|
# parse gff attributes
|
|
#my $psid = $extra =~ /Name \"(\w+)\"/ ? $1 : undef;
|
|
my $aclist = [];
|
|
if ( $extra =~ /AccessionNumbers (.+?) ;/g ) {
|
|
@$aclist = map { s/['"]//g; $_ } split / /, $1
|
|
}
|
|
my $leveln = $extra =~ /Level ([-\d]+)/ ? $1 : undef;
|
|
my $levelt = $extra =~ /LevelTag \"(.+?)\"/ ? $1 : undef;
|
|
my $rawscore = $extra =~ /RawScore (\d+)/ ? $1 : undef;
|
|
my $pffrom = $extra =~ /FeatureFrom ([-\d]+)/ ? $1 : undef;
|
|
my $pfto = $extra =~ /FeatureTo ([-\d]+)/ ? $1 : undef;
|
|
my $subseq = $extra =~ /Sequence \"(.+?)\"/ ? $1 : undef;
|
|
my $seqde = $extra =~ /SequenceDescription \"(.+?)\"/ ? $1 : '';
|
|
|
|
push @{ $hits_by_psac_by_seq->{ $id }->{ $psac } },
|
|
[ $subseq, $from, $to, $psac, $pffrom, $pfto, $rawscore,
|
|
$nscore, $leveln, $levelt, '' ];
|
|
$info_by_seqid->{ $id }->{ de } = $seqde;
|
|
$info_by_seqid->{ $id }->{ acs } = $aclist;
|
|
}
|
|
}
|
|
# for each seq do pp
|
|
foreach my $seqid ( keys %$hits_by_psac_by_seq ) {
|
|
postProcess( $hits_by_psac_by_seq->{ $seqid } )
|
|
unless $opt_no_postprocessing;
|
|
# if no pp: usefull only for format conversion
|
|
foreach my $psac ( sort { $a cmp $b }
|
|
keys %{ $hits_by_psac_by_seq->{ $seqid } } ) {
|
|
# display hits (foreach motif ac)
|
|
my $hits = $hits_by_psac_by_seq->{ $seqid }->{ $psac };
|
|
my $seqde = $info_by_seqid->{ $seqid }->{ de } || '';
|
|
my $aclist = $info_by_seqid->{ $seqid }->{ acs } || [];
|
|
my $psid = $psinfo_bypsac->{ $psac }->[ 1 ] || $psac;
|
|
my $psde = $psinfo_bypsac->{ $psac }->[ 3 ] || '';
|
|
# n.b. psde not in gff source, but in motif file
|
|
dispHits( undef, '', $hits, $seqid, $seqde, $aclist, $psac,
|
|
$psid, $psde );
|
|
}
|
|
$hits_by_psac_by_seq->{ $seqid } = undef;
|
|
$info_by_seqid->{ $seqid } = undef;
|
|
}
|
|
}
|
|
|
|
|
|
# -------------------------- scan methods --------------------------
|
|
|
|
# scan entries in specified sequence file 'collection' against motif
|
|
# collection and display results through called subs)
|
|
sub scanSeqFile {
|
|
my $seqfile = shift or return;
|
|
|
|
my $entry = "";
|
|
my $opt_fasta;
|
|
my @motifs_4_normal_scan;
|
|
my @motifs_4_pfsearch;
|
|
my $psinfo_bypsac = {};
|
|
map {
|
|
# store motif info by motif ac
|
|
$_->[0] and $psinfo_bypsac->{$_->[0]} = $_;
|
|
$opt_pfsearch && $_->[2] && $_->[2] eq 'MATRIX' ?
|
|
push @motifs_4_pfsearch, $_ : push @motifs_4_normal_scan, $_
|
|
} @MotifInfo;
|
|
my $all_hits_bypsac_byseqid = {};
|
|
# -W PFSEARCH SCAN (every profile against sequence collection)
|
|
if (grep {$_->[7]} @motifs_4_pfsearch) {
|
|
# with pfsearch (faster to scan each profile
|
|
# against a sequence collection):
|
|
# use whole sequence collection files, do not extract each entry!
|
|
foreach my $prf_info (@motifs_4_pfsearch) {# for each Motifinfo
|
|
next if $prf_info->[5]; # next if skip...
|
|
my $psac = $prf_info->[0] or next;
|
|
# needs an associated tmp profile file
|
|
my $profile_tmp_file = $prf_info->[7] or next;
|
|
my $psid = $prf_info->[1] || '';
|
|
# scan seq collection file against 1 profile with pfsearch
|
|
#! with pfsearch do_profile_scan
|
|
# will return hits in a hash (ref) by seq ac|id
|
|
# (not by psac like with pfscan)!
|
|
my $pfsearchhits =
|
|
do_profile_scan($profile_tmp_file, $seqfile, undef, $prf_info);
|
|
# group all diff hits by psac and by seq :
|
|
foreach my $seq_id (keys %$pfsearchhits) {
|
|
# for each hit
|
|
foreach my $hit (@{$pfsearchhits->{$seq_id}}) {
|
|
# set Name field (in -o gff) to psid (instead of seqid)
|
|
# = same as with normal scan
|
|
# (my $seqid=$seq_id)=~s/^\w+\|//;
|
|
# FIXME: should be done at the do_profile_scan level
|
|
# fix seqid (do_profile_scan sometimes
|
|
# returns ac|id instead of id)
|
|
|
|
# EXPGE-313
|
|
next if $opt_filterheader && $hit->[10] !~ /$opt_filterheader/o;
|
|
$hit->[3] = $psid;
|
|
push @{$all_hits_bypsac_byseqid->{$seq_id}->{$psac}},
|
|
$hit;
|
|
}
|
|
}
|
|
}
|
|
# post process and show results (all profile hits grouped by protein)
|
|
# note: $all_hits_bypsac_byseqid only used for -w pfsearch profile.
|
|
# for normal scan post-processed and displayed is done in scanSeq!
|
|
# side effect: with -w pfsearch,
|
|
# post processing between profiles and patterns
|
|
# won't work
|
|
|
|
# for each matched sequence ac
|
|
foreach my $seq_id ( keys %$all_hits_bypsac_byseqid ) {
|
|
my $pfhits_bypsac = $all_hits_bypsac_byseqid->{ $seq_id } or next;
|
|
# do postprocessing (on all sequence hits)
|
|
postProcess( $pfhits_bypsac ) unless $opt_no_postprocessing;
|
|
foreach my $psac ( keys %$pfhits_bypsac ) {# display hits
|
|
my $psid = $psinfo_bypsac->{$psac}->[1] || $psac;
|
|
# won't be displayed (in pfsearch mode)!
|
|
# but in case this will be changed
|
|
my $psde = $psinfo_bypsac->{$psac}->[3] || '';
|
|
dispHits( undef, '',
|
|
$all_hits_bypsac_byseqid->{$seq_id}->{$psac},
|
|
$seq_id, '', undef, $psac, $psid, $psde );
|
|
}
|
|
}
|
|
}
|
|
# 'NORMAL' SCAN (every sequence against profile collection)
|
|
if ( @motifs_4_normal_scan ) {# scan entries in sequence file collection
|
|
my $seqfile_h = new IO::File $seqfile
|
|
or die "Cannot open $seqfile: $!";
|
|
while ( <$seqfile_h> ) {# read sequence file
|
|
my $line = $_;
|
|
if ( /^>(.*)/ ) {# entry is in fasta
|
|
#(my $seqid = $1) =~ s/\s+\.*$//;
|
|
# $all_hits_bypsac_byseqid->{$seqid}
|
|
# = scanFromFastaEntry($entry);
|
|
scanFromFastaEntry( $entry );
|
|
# note: do not store hits (by ac) into
|
|
# $all_hits_bypsac_byseqid
|
|
# might use too much memory!; side effect: with -w pfsearch,
|
|
# post processing between profiles and patterns won't work
|
|
$opt_fasta = 1;
|
|
$entry = '';
|
|
}
|
|
$entry .= $line;
|
|
# entry separator (in Swiss-Prot format) => an entry has been
|
|
# fully read
|
|
if ( /^\/\// ) {
|
|
$opt_fasta = 0;
|
|
my @id = $entry =~ /^\s*ID\s+(\w+)/mg;
|
|
my $ac_lines;
|
|
$ac_lines .= $1 while $entry =~ /^\s*AC\s+(.*)/mg;
|
|
my @ac;
|
|
while ( $ac_lines =~ /(\w+)/g ) { push @ac, $1; }
|
|
if ( @id ) {
|
|
if ( not (@entries) or
|
|
grep { my $ent = $_;
|
|
grep{ $_ eq $ent } @id, @ac } @entries) {
|
|
my $id = $id[ 0 ];
|
|
# get entry DE:
|
|
my $de = #? 'new' format
|
|
$entry =~ /^DE (?:Rec|Sub)Name: Full=(.+);/m ?
|
|
$1 : '';
|
|
unless ( $de ) {# old DE format
|
|
my @de = $entry =~ /^\s*DE\s+(.+)/mg;
|
|
my $add_space = 0;
|
|
for ( @de ) {
|
|
$de .= " " if $add_space;
|
|
$de .= $_;
|
|
$add_space = !/-$/;
|
|
}
|
|
}
|
|
|
|
if ($entry =~ /^\s*SQ\s+SEQUENCE\b.*\n((.+\n)+)/m) {
|
|
my $sq = $1;
|
|
$sq =~ tr/A-Z//cd;
|
|
# $all_hits_bypsac_byseqid->{$id} =
|
|
# scanSeq($id, \@ac, $de || $id, $sq);
|
|
# note: do not store hits (by ac) into
|
|
# $all_hits_bypsac_byseqid would use too much
|
|
# memory; side effect: with -w pfsearch,
|
|
# post processing between a profile and and
|
|
# pattern won't work
|
|
scanSeq($id, \@ac, $de || $id, $sq);
|
|
}
|
|
else {
|
|
warn "No sequence found in entry $id";
|
|
}
|
|
}
|
|
}
|
|
# ignore entries which have "id" in lowercase
|
|
elsif ($entry =~ /^\s*id /m) {
|
|
}
|
|
elsif ($entry =~ /(.*\S.*)/) {
|
|
warn "Bad sequence found in file, first line: $1\n";
|
|
$errcode = 1;
|
|
}
|
|
$entry = "";
|
|
}
|
|
}
|
|
close $seqfile_h;
|
|
if ($entry =~ /^>/) {# process last fasta entry
|
|
scanFromFastaEntry($entry);
|
|
}
|
|
elsif ($entry =~ /(.*\S.*)/) {
|
|
warn "Bad sequence found in file, first line : $1\n";
|
|
$errcode = 1
|
|
}
|
|
}
|
|
}
|
|
|
|
# scan from 1 fasta entry (from a sequence entry collection)
|
|
sub scanFromFastaEntry {
|
|
my $entry = shift or return;
|
|
return unless $entry =~ s/^>((\S*).*)\n//;
|
|
my ($fasta_header, $primary_id) = ($1, $2);
|
|
return if defined($opt_filterheader)
|
|
and $fasta_header !~ /$opt_filterheader/o;
|
|
if (not (@entries) or grep {$_ eq $primary_id} @entries) {
|
|
$entry =~ tr/A-Z//cd;
|
|
return scanSeq($primary_id, [], $fasta_header, $entry);
|
|
}
|
|
}
|
|
|
|
# scans one sequence against prosite motifs (only normal scan; if
|
|
# opt_pfsearch, all scans are done before in scanSeqFile)
|
|
sub scanSeq {
|
|
my ($id, $aclist, $de, $sq) = @_;
|
|
if ($opt_reverse) {# reverses sequence (if opt on)
|
|
$sq = reverse $sq;
|
|
}
|
|
if ($opt_shuffle) {# shuffles sequence (if opt on)
|
|
srand 0;
|
|
my @seq = grep {$_ ne "\n"} split(//,$sq);
|
|
$sq = "";
|
|
for (my $start_win = 0; $start_win < @seq;
|
|
$start_win += $opt_shuffle) {
|
|
my $stop_win = $start_win + $opt_shuffle - 1;
|
|
$stop_win = @seq - 1 if $stop_win >= @seq;
|
|
my @residues = @seq[$start_win..$stop_win];
|
|
while (@residues) {
|
|
$sq .= splice(@residues, int(rand(scalar @residues)) ,1);
|
|
}
|
|
}
|
|
}
|
|
# scan sequence against profile first
|
|
my $all_pfhits = {};
|
|
# if there are some profile to be scanned...
|
|
if ($scan_profiles && !$opt_pfsearch) {
|
|
# in normal scan (not -w pfsearch) only 1 temp file is used (or
|
|
# user specified prosite files)
|
|
my $file_source = ($last_profile_tmp_filename ?
|
|
[$last_profile_tmp_filename] : \@prosite_files);
|
|
foreach my $profile_file (@$file_source) {#
|
|
my $pfhits = do_profile_scan($profile_file, undef, $sq);
|
|
# store results (do_profile_scan result is a hash ref
|
|
# to [[hitdata],[],[]...] by profile AC)
|
|
foreach my $ps_acid (keys %$pfhits) {
|
|
(my $psac = $ps_acid) =~ s/\|.+//;
|
|
$all_pfhits->{$psac} = $pfhits->{$ps_acid};
|
|
}
|
|
}
|
|
}
|
|
# loop through each motif, scan sequence against MATRIX, builds hit
|
|
# structure
|
|
my $hits_by_motif_ac = {};
|
|
foreach (@MotifInfo) {# for each motif info 'data'
|
|
my ($psac, $psid, $type, $psde, $pat, $skip) = @$_;
|
|
next if $skip;
|
|
my $hits;
|
|
# MATRIX (profile) (profiles were run just before...)
|
|
if ($type eq "MATRIX") {
|
|
# if opt_pfsearch: scan and results display already done
|
|
# (in scanSeqFile sub)
|
|
next if $opt_pfsearch;
|
|
$hits = $hits_by_motif_ac->{$psac} = $all_pfhits->{$psac} or next;
|
|
}
|
|
else {# PATTERN
|
|
warn("Empty pattern for $psac\n"), next unless $pat;
|
|
$hits = scanPatternWrapper([$pat, $sq, $scan_behavior,
|
|
$opt_max_x, $opt_miniprofiles], $psid);
|
|
$hits_by_motif_ac->{$psac} = $hits if @$hits;
|
|
}
|
|
dispHits(undef, $sq, $hits, $id, $de, $aclist, $psac, $psid, $psde)
|
|
if $opt_no_postprocessing;
|
|
# if no pp matches can be displayed right now (as they go)
|
|
# otherwise, first have to collect all matches (this loop), then pp
|
|
}
|
|
unless ($opt_no_postprocessing) {
|
|
# if pp: only can be done after having collected all matches
|
|
# do post-processing (between all hits within one entry)
|
|
# note: messy: pp.+disp. performed here, but also before in scanSeqFile
|
|
# for -w pfsearch profiles!
|
|
postProcess($hits_by_motif_ac);
|
|
# display results (after post processing) for each motif
|
|
# (keep order of motifs in prosite motif file or order of
|
|
# specified (-p) motifs)
|
|
foreach (@MotifInfo) {
|
|
my ($psac, $psid, $type, $psde, $pat, $skip) = @$_;
|
|
my $hits = $hits_by_motif_ac->{$psac} or next;
|
|
dispHits(undef, $sq, $hits, $id, $de, $aclist, $psac, $psid, $psde);
|
|
}
|
|
}
|
|
return $hits_by_motif_ac;
|
|
}
|
|
|
|
|
|
# run pfscan on 1 sequence against profile collection or pfsearch (if -w...):
|
|
# 1 profile against sequence collection sequence can be specified either
|
|
# as a filename or a string
|
|
sub do_profile_scan {
|
|
my ($PROSITE, $seqfile_to_scan, $sequence, $prf_info) = @_;
|
|
#return {} unless grep {$_->[2] eq "MATRIX"} @MotifInfo;
|
|
my $PFSCAN_TMP = tmpnam();
|
|
my $level_arg = defined($opt_level) ? "$opt_level" : "0";
|
|
my $level_min = 0;
|
|
my(@pre_command, @post_command);
|
|
my $use_pfscan_v3 = $PFSCAN =~ /V3/ ? 1 : 0; # fugly! but needed as new pfscanV3 has different cmd line!
|
|
if ($opt_pfsearch) {# pfsearch (1 prf against a seq collection)
|
|
# retrive cut off value from profile to set pfsearch -C option
|
|
# (not needed with pfsearch 2.3 (can use -l <level>)...
|
|
# but to be compatible with 2.2...)
|
|
my $pfsearch_cutoff;
|
|
# get cut-off corresponding to level -1
|
|
# (needed for the detection of the repeats)
|
|
map {
|
|
$pfsearch_cutoff = ($opt_raw ? $_->{SCORE} : $_->{N_SCORE})
|
|
if ($_->{LEVEL} eq '-1');
|
|
$level_min = $_->{LEVEL} if ($_->{LEVEL} < $level_min);
|
|
} @{$prf_info->[6]}
|
|
if ($prf_info && $prf_info->[6]);
|
|
# if option pfsearch is selected, get the user specified C=?? parameter
|
|
# if not defined: pfsearch is run as default at L=-1
|
|
# (pfsearch aut. detects if C is SCORE (integer) or N_SCORE(float))
|
|
my $cutoff = defined($opt_cutoff) ? "$opt_cutoff" :
|
|
"$pfsearch_cutoff";$cutoff ||= '0';
|
|
# detect format
|
|
my $fasta = $sequence ? "-f" : "";
|
|
if ( open( DETECT, $seqfile_to_scan ) ) {
|
|
while ( <DETECT> ) {
|
|
next unless /\S/;
|
|
$fasta = /^\s*>/ ? "-f" : "";
|
|
last;
|
|
}
|
|
close DETECT;
|
|
}
|
|
@pre_command = ( $use_pfsearchV3 ?
|
|
"$opt_pfsearch -o1 -L $level_min $PROSITE" : "$opt_pfsearch $fasta -lxz $PROSITE" );
|
|
# p.s. if input is not fasta, pfsearchV3 will fail!
|
|
@post_command = ( $use_pfsearchV3 ? "" : "C=$cutoff" );
|
|
}
|
|
else {# 'normal' pfscan (1seq against a prf collection)
|
|
if ($level_arg eq 0) { # yes default scan l=0 is internally run at -1 (so that postprocessing could consider weak matches)
|
|
$level_arg = -1;
|
|
}
|
|
@pre_command = $use_pfscan_v3 ? "$PFSCAN --matrix-only -o4 -L$level_arg $PROSITE" : "$PFSCAN -flxz -v -C$level_arg";
|
|
@post_command = $use_pfscan_v3 ? "" : $PROSITE;
|
|
}
|
|
my $out;
|
|
if ( $use_pfsearchV3 || $NO_DIRECT_PIPE || defined $seqfile_to_scan || $use_pfscan_v3 ) {
|
|
my $seqfile;
|
|
unless (defined $seqfile_to_scan) {
|
|
$seqfile = tmpnam();
|
|
open SEQ_TMP, ">$seqfile" or die "Cannot create $seqfile : $!";
|
|
print SEQ_TMP ">seq for pfscan\n";
|
|
print SEQ_TMP "$1\n" while $sequence =~ /(.{1,60})/g;
|
|
close SEQ_TMP;
|
|
}
|
|
else {
|
|
$seqfile = $seqfile_to_scan;
|
|
}
|
|
my $cmd = "@pre_command $seqfile @post_command > $PFSCAN_TMP";
|
|
# launch pftool scan command
|
|
system $cmd and die "Could not execute $cmd: $!";
|
|
unlink $seqfile unless defined $seqfile_to_scan;
|
|
my $pfscan_fh = new IO::File($PFSCAN_TMP)
|
|
or die "Cannot open $PFSCAN_TMP: $!";
|
|
$out = scanProfiles($pfscan_fh, $level_min-1);# parse scan output
|
|
close $pfscan_fh or die "Error $? with $PFSCAN_TMP";
|
|
}
|
|
else {
|
|
# directly feed data to pfscan via pipe, do not use any temporary
|
|
# files
|
|
require IPC::Open2;
|
|
my ($reader, $writer);
|
|
my $cmd = "@pre_command - @post_command";
|
|
my $pid =
|
|
eval {
|
|
IPC::Open2::open2($reader, $writer, $cmd)
|
|
or die "Could not fork pipe to $cmd: $!";
|
|
};
|
|
if ($@) {
|
|
die "$@\n" . '>' x 62 . "\nERROR: 'pfscan' execution failed.
|
|
Check pfscan is in your PATH\n" . '>' x 62 . "\n";
|
|
}
|
|
local $/ = \32767;# buffer size
|
|
print $writer ">seq for pfscan\n";
|
|
print $writer "$1\n" while $sequence =~ /(.{1,60})/g;
|
|
close $writer;
|
|
$out = scanProfiles($reader, $level_min-1);
|
|
close $reader or die "Error $? with $cmd";
|
|
waitpid $pid, 0; #avoid defunct kid processes
|
|
die "Could not execute $cmd; returned: $?" if $?;
|
|
}
|
|
unlink $PFSCAN_TMP;
|
|
if ($opt_format eq "msa") {# if output format is MSA run psa2msa
|
|
for my $ac (keys %$out) {
|
|
# for each prosite profile ac write to tmp file (reuse $PFSCAN_TMP)
|
|
open PFSCANTMP, ">$PFSCAN_TMP"
|
|
or die "Cannot create $PFSCAN_TMP : $!";
|
|
# for each hit (by ps ac)
|
|
for (my $i = 0; $i < @{$out->{$ac}}; $i++) {
|
|
my $hit = $out->{$ac}->[$i];
|
|
print PFSCANTMP ">$i\n$hit->[0]";
|
|
}
|
|
close PFSCANTMP;
|
|
my $PSA2MSA_TMP = tmpnam();
|
|
my $cmd = "$PSA2MSA $PFSCAN_TMP > $PSA2MSA_TMP";# run psa2msa
|
|
system $cmd and die "Cannot execute $cmd";
|
|
open MSATMP, $PSA2MSA_TMP or die "Cannot read $PSA2MSA_TMP : $!";
|
|
my %msa;
|
|
my $cur_pos;
|
|
while(defined (local $_ = <MSATMP>)) {# read psa2msa result output
|
|
if (/^>(\d+)/) {
|
|
$cur_pos = $1;# hit number
|
|
}
|
|
elsif (defined $cur_pos) {
|
|
$msa{$cur_pos} .= $_;# sequence 'alignment'
|
|
}
|
|
}
|
|
close MSATMP;
|
|
unlink $PSA2MSA_TMP;
|
|
while (my ($number, $seq) = each %msa) {
|
|
# change match sequence (?) for hit (for prosite $ac)
|
|
# number $number
|
|
$out->{$ac}->[$number]->[0] = $seq;
|
|
}
|
|
}
|
|
}
|
|
return $out;
|
|
}
|
|
|
|
# replace hits on a repeat region by individual repeat elements
|
|
# (with pfsearch/pfscan 2.3)
|
|
sub unroll_hits {
|
|
my ($hits) = @_;
|
|
return unless(grep {$_->[11] && @{$_->[11]}} @$hits);
|
|
for (my $i = 0; $i < @$hits; $i++) {
|
|
my ($subseq, $from, $to, $pfid, $pffrom, $pfto, $rawscore,
|
|
$nscore, $leveln, $levelt, $seqde, $subhits) = @{$hits->[$i]};
|
|
next unless($subhits && @$subhits);
|
|
map {$_->[8] = $leveln; $_->[9] = $levelt; $_->[10] = $seqde}
|
|
@$subhits;
|
|
splice @$hits, $i--, 1, @$subhits;
|
|
}
|
|
}
|
|
|
|
sub scanPatternWrapper {
|
|
my ($args, $id) = @_;
|
|
my $out = scanPattern(@$args);
|
|
# in PSA format, remove the '.' character from inserts.
|
|
# these can be reintroduced with the 'psa2msa' program.
|
|
if ($opt_format eq "psa") {
|
|
$_->[0] =~ s/\.//g for @$out;
|
|
}
|
|
if ($id) {
|
|
$_->[3] = $id for @$out;
|
|
}
|
|
return $out;
|
|
}
|
|
|
|
sub prositeToRegexpWrapper {
|
|
my $out = prositeToRegexp(@_);
|
|
unless (defined $out) {
|
|
print STDERR "ps_scan.pl: Syntax error in pattern".
|
|
" at position $Prosite::errpos\n";
|
|
print STDERR "$_[0]\n";
|
|
print STDERR " " x $Prosite::errpos, "^--- $Prosite::errstr\n";
|
|
exit 1;
|
|
}
|
|
return $out;
|
|
}
|
|
|
|
# -------------------------- post processing --------------------------
|
|
|
|
sub pp_scan {
|
|
my ($hits, $all_psac_in_seq, $intra_not_inter) = @_;
|
|
return unless($hits && $all_psac_in_seq);
|
|
my $seen = {};
|
|
foreach my $target_psac (@$all_psac_in_seq) {
|
|
# loops through all motifs in result structure
|
|
my $hit_target_set = $hits->{$target_psac} or next;
|
|
@$hit_target_set or next;
|
|
foreach my $potential_pp (grep {($intra_not_inter ? $_->{effector} eq
|
|
$target_psac : $_->{effector} ne $target_psac) }
|
|
@{$postProcessingByPSAC->{$target_psac}}) {
|
|
# loops through all motif associated pp
|
|
# (order = order of PP lines...)
|
|
# (+filter for intra/inter pp...)
|
|
#print STDERR "##PP ($intra_not_inter) target [$target_psac]",
|
|
#"-> pp effector[$potential_pp->{effector}] ",
|
|
#"type[$potential_pp->{type}] value[$potential_pp->{value}]\n";
|
|
my $effector_psac = $potential_pp->{effector} or next;
|
|
my $pp_type = $potential_pp->{type} or next;
|
|
next if $seen->{$pp_type}->{$target_psac}->{$effector_psac};
|
|
#next if !$allowBidirectionalPP->{$pp_type} &&
|
|
# $seen->{$target_psac}->{$effector_psac};
|
|
# next if 'reverse' (target<->effector) pp was already seen!
|
|
# might happen with COMPETE_ PPs where PP is specified in
|
|
# both effector & target...
|
|
$seen->{$pp_type}->{$effector_psac}->{$target_psac} = 1
|
|
unless $allowBidirectionalPP->{$pp_type};
|
|
|
|
my $hit_effector_set = $hits->{$effector_psac};
|
|
next unless $hit_effector_set or
|
|
$potential_pp->{allow_no_effector_matches};
|
|
next if $hit_effector_set && !@$hit_effector_set;
|
|
|
|
my $pp_value = $potential_pp->{value};
|
|
# perform post-processing on target:
|
|
if ($postProcessDispatchTable->{$pp_type}) {
|
|
&{$postProcessDispatchTable->{$pp_type}}
|
|
($hit_target_set,$hit_effector_set,$pp_value);
|
|
# call pp subs (from dispatch table)
|
|
} else {
|
|
print STDERR "unknown post-processing key [$pp_type]: ignored";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
# post process results (inside same entry)
|
|
sub postProcess {
|
|
my $hits = shift or return;# input: ref to hash of hits by motif ac...
|
|
my @all_psac_in_seq = sort {
|
|
($motifRank4PostProcessing->{$b}||0) <=>
|
|
($motifRank4PostProcessing->{$a}||0)
|
|
} keys %$hits;
|
|
# perform intra motif pp first (repeats...)
|
|
pp_scan($hits, \@all_psac_in_seq, 1);
|
|
# perform inter motif pp
|
|
pp_scan($hits, \@all_psac_in_seq, 0);
|
|
}
|
|
|
|
|
|
#############################################################################
|
|
# CORE
|
|
|
|
my $isPPLinkedTo = {};
|
|
my $isInSpecifiedPPGroup = {};
|
|
|
|
sub addPPLinkedMotifs2FetchStruct {
|
|
foreach my $psfile (@prosite_files) {
|
|
# loop through all specified prosite files
|
|
open PSFILE, $psfile or die "Cannot open $psfile : $!";
|
|
my $ac = '';
|
|
while (<PSFILE>) {# read prosite file
|
|
$ac = $1 if m/^AC (\w+);/;
|
|
$ac = '' if m/^\/\//;
|
|
map {
|
|
s/\(.+\)$//;
|
|
$isPPLinkedTo->{$ac}->{$_} = 1, $isPPLinkedTo->{$_}->{$ac} = 1
|
|
if ($_ && $ac ne $_)
|
|
} split(';\s*',$1)
|
|
if $ac && /^PP \/?\w+:\s*([^\r\n]+)/;
|
|
}
|
|
}
|
|
close PSFILE;
|
|
# collect all profile 'linked by pp' to specified psac
|
|
my @followpp_stack = @followpp;
|
|
while (my $psac = pop @followpp_stack) {# for each specified psac
|
|
$specifiedPrositeMotifByAc->{$psac} = 1;
|
|
# add motif ac to 'motif to fetch/use' struct
|
|
foreach my $linked_motifac (keys(%{$isPPLinkedTo->{$psac}})) {
|
|
# add linked profile ac to profile to be hidden at display:
|
|
$hideMotifByPSAC->{$linked_motifac} = 1
|
|
unless grep {$linked_motifac eq $_} @followpp,
|
|
@motifAC_or_userpattern;
|
|
# add linked profile ac to list of profile to be 'pp-followed'
|
|
# (unless already in specifiedPrositeMotifByAc struct):
|
|
push @followpp_stack, $linked_motifac
|
|
unless $specifiedPrositeMotifByAc->{$linked_motifac};
|
|
}
|
|
}
|
|
}
|
|
|
|
# processes PROSITE motif data/info (into global @MotifInfo)
|
|
sub processMotif {
|
|
my $ps_entry = shift or return;# prosite entry (string)
|
|
# parse entry
|
|
my ($ac, $id, $type, $de, $pa, $rule, $pdoc, $skipflag,
|
|
$tax, $rep, $sites, $cutoffs, $pps) = parseProsite($ps_entry);
|
|
$ac or $id or die "ERROR: can't parse prosite entry";
|
|
# store post processing data (if any)
|
|
map {
|
|
my $type = $_->{type};
|
|
map {
|
|
my ($effector,$value) = (m/^(\w+\d+)(?:\((.+)\))*/);
|
|
push @{$postProcessingByPSAC->{$ac}},
|
|
{'type'=>$type,'effector'=>$effector,'value'=>$value};
|
|
# print "PP: ac[$ac], type[$type], effector[$effector],
|
|
# value[$value]\n";
|
|
}
|
|
@{$_->{discriminators}}
|
|
} @$pps
|
|
if ($pps && @$pps);
|
|
my $mp_ac;
|
|
if ($type eq "PATTERN" and $opt_miniprofiles
|
|
and $ac=~/PS(\d{5})/ and !$skipflag
|
|
and $opt_format ne "sequence" and $opt_format ne "matchlist") {
|
|
$mp_ac = "MP$1";
|
|
unshift @{$postProcessingByPSAC->{$ac}},
|
|
{'type'=>'EVALUATED_BY', 'effector'=>$mp_ac,
|
|
'value'=>$mp_ac, 'allow_no_effector_matches'=>1 };
|
|
}
|
|
# store implicit 'repeat' post-processing data
|
|
# (RR, R? in Text element of MA CUT_OFF profile lines)
|
|
my $sum=0.0;
|
|
foreach my $cutoff (@$cutoffs) {
|
|
$sum += $cutoff->{N_SCORE} || 0.0 if ($cutoff->{LEVEL} >= -1);
|
|
if ($cutoff->{TEXT} eq "'R?'") {
|
|
# RDM1 pp: promote weak match (add 1 to level)
|
|
# if there is at least one hit at level>=0
|
|
unshift @{$postProcessingByPSAC->{$ac}},
|
|
{'type'=>'PROMOTER_RDM1','effector'=>$ac,'value'=>undef};
|
|
}
|
|
elsif ($cutoff->{TEXT} eq "'RR'") {
|
|
# RDM1 or RDM2 pp: RDM2=
|
|
unshift @{$postProcessingByPSAC->{$ac}},
|
|
{'type'=>'PROMOTER_RDM1_OR_RDM2',
|
|
'effector'=>$ac,'value'=>$sum};
|
|
}
|
|
elsif ($opt_rep_pp_4allprofiles && $cutoff->{TEXT} eq "'!'") {
|
|
# if -a option (opt_rep_pp_4allprofiles):
|
|
# apply RDM1 RDM2 on any hits
|
|
# from a level0 '!' profile (within a seq)!
|
|
unshift @{$postProcessingByPSAC->{$ac}},
|
|
{'type'=>'PROMOTER_RDM1_OR_RDM2_4!REP',
|
|
'effector'=>$ac,'value'=>$sum};
|
|
}
|
|
}
|
|
# sets skipflag, knowfalsepos flags (hashes by ac)
|
|
$SkipFlag{$ac} = 1 if $skipflag;
|
|
my $skip = $skipflag && $opt_skip;
|
|
$skip = !$skipflag if $opt_skiponly;
|
|
if ($ps_entry =~ /^DR.*, T;/m) {
|
|
my $nbfp = 0;
|
|
for my $line ($ps_entry =~ /^DR(.*)/mg) {
|
|
$nbfp += $line =~ s/, F;//g;
|
|
}
|
|
$KnownFalsePos{$ac} = $nbfp;
|
|
}
|
|
# build motif info 'arrays'
|
|
if ($type eq "PATTERN") {# for PATTERN
|
|
$scan_pattern = 1;
|
|
#die "ERROR: can't scan against PATTERN motifs
|
|
#when -w (pfsearch) option is used"
|
|
#if ($opt_pfsearch);
|
|
$scan_pattern = 1;
|
|
push @MotifInfo,
|
|
[$ac, $id, $type, $de,
|
|
prositeToRegexpWrapper($pa, $opt_nongreedy,
|
|
defined $opt_max_x && !$opt_max_x ? 1 : 0),
|
|
$skip];
|
|
}
|
|
elsif ($type eq "MATRIX") {# for profiles
|
|
$scan_profiles = 1;
|
|
if ($specifiedPrositeMotifByAc || $opt_pfsearch) {
|
|
# if not all motifs are used (when some were specified) or
|
|
# if $opt_pfsearch: save data to tmp profile
|
|
if (!$last_profile_tmp_filename || $opt_pfsearch) {
|
|
# open new profile temp
|
|
close PROFILE_TMP if ($last_profile_tmp_filename);
|
|
# close previous one
|
|
$last_profile_tmp_filename=tmpnam();
|
|
# get new temp file name
|
|
open PROFILE_TMP, ">$last_profile_tmp_filename"
|
|
or die "Cannot open $last_profile_tmp_filename: $!";
|
|
}
|
|
# save profile to temp file
|
|
print PROFILE_TMP $ps_entry
|
|
or die "can't print to $last_profile_tmp_filename: $!\n";
|
|
}
|
|
push @MotifInfo,
|
|
[$ac, $id, $type, $de, undef, $skip, $cutoffs,
|
|
$last_profile_tmp_filename];
|
|
}
|
|
else {
|
|
warn "Unknown prosite entry type $type";
|
|
}
|
|
}
|
|
|
|
sub main {
|
|
|
|
######################## Motif data parsing/processing #####################
|
|
addPPLinkedMotifs2FetchStruct() if !$opt_no_postprocessing && @followpp;
|
|
my $rank = 0;
|
|
for my $psfile ( @prosite_files ) {
|
|
# loop through all specified prosite files
|
|
open PSFILE, $psfile or die "Cannot open $psfile : $!";
|
|
my $ps_entry = ""; my $ac = ""; my $id = ""; my $type = "";
|
|
my $pos = 0;
|
|
|
|
PROSITE: while ( <PSFILE> ) {# read prosite file
|
|
$ps_entry .= $_;
|
|
$ac=$1 if !$ac && m/^AC (\w+);/;
|
|
$id=$1, $type=$2
|
|
if !$id && !$type && m/^ID (\w+);\s+(\w+)\./;
|
|
|
|
if (/^\/\//) {# end of an entry: PROCESS MOTIF DATA
|
|
my $use_motif = 1;
|
|
$use_motif = 0 if (!$id || !$ac);# skip 'bad' entries
|
|
if ($use_motif && $opt_noprofiles && $type eq "MATRIX") {
|
|
# skip profiles if opt_noprofiles
|
|
$use_motif = 0;
|
|
die "ERROR: profile(s) specified".
|
|
" with -p or -f option, but -r".
|
|
" (do not scan profile) option used"
|
|
if $specifiedPrositeMotifByAc;
|
|
}
|
|
# skip non profiles if opt_onlyprofiles
|
|
if ($use_motif && $opt_onlyprofiles && $type ne "MATRIX") {
|
|
$use_motif = 0;
|
|
die "ERROR: pattern(s) specified with -p or -f option,".
|
|
" but -m (only scan profile) option used"
|
|
if $specifiedPrositeMotifByAc;
|
|
}
|
|
$use_motif = 0 if $use_motif && $specifiedPrositeMotifByAc &&
|
|
!$specifiedPrositeMotifByAc->{$ac};
|
|
# if AC were specified (PSxxxx format), skip entry that do
|
|
# not match
|
|
delete ($specifiedPrositeMotifByAc->{$ac})
|
|
if $specifiedPrositeMotifByAc;
|
|
# delete specified AC so that we can check later if they
|
|
# were all found...
|
|
$motifRank4PostProcessing->{$ac} =++ $rank
|
|
if $use_motif && !$opt_no_postprocessing;
|
|
# ... precedence order for PostProcessing = inverse of motif
|
|
# position in file (lower will be pp first)
|
|
processMotif( $ps_entry ) if $use_motif;
|
|
# process motif data into @MotifInfo
|
|
$ps_entry = $ac = $id = $type = "";
|
|
$pos = tell PSFILE;
|
|
}
|
|
}
|
|
close PSFILE;
|
|
}
|
|
close PROFILE_TMP if $last_profile_tmp_filename;
|
|
if ( $specifiedPrositeMotifByAc ) {
|
|
# look if there are unfound specified motif AC
|
|
my @notfound;
|
|
foreach my $ac_id_not_found ( keys( %$specifiedPrositeMotifByAc ) ) {
|
|
push @notfound, $ac_id_not_found;
|
|
}
|
|
die "Prosite entry [@notfound] not found in".
|
|
" specified prosite file(s)\n"
|
|
if @notfound;
|
|
}
|
|
# process -p option user defined patterns
|
|
my $user_ctr = 1;
|
|
for ( @userpat ) {# transform userpat into correct array representation
|
|
my $i = "0" x ( 3 - length( $user_ctr ) ) . $user_ctr++;
|
|
push @MotifInfo, [ "USER$i", undef, "USER", undef,
|
|
prositeToRegexpWrapper( $_, $opt_nongreedy, 1 ), 0 ];
|
|
}
|
|
|
|
###################### Perform scan (and show results) #####################
|
|
|
|
if ( @external_gff_files ) {
|
|
# do not scan: just (post-process and/or convert) gff results
|
|
process_external_gff();
|
|
}
|
|
elsif ( $opt_format ne "matchlist" ) {# SCAN (non matchlist output format)
|
|
# get & show motif hits
|
|
unshift( @ARGV, '-' ) unless @ARGV;
|
|
while ( my $seqfile = shift @ARGV ) {
|
|
# for each specified sequence file: scan
|
|
scanSeqFile($seqfile);
|
|
}
|
|
}
|
|
else {# SCAN (matchlist output format) (fugly!?)
|
|
# FIXME: used!? + seen in comments: only works with sp flat input!?
|
|
showMatchList();
|
|
}
|
|
|
|
############################ exit ############################
|
|
|
|
my @tmp_mp_files = values %$files_miniprofiles;
|
|
# delete each temporary miniprofile file
|
|
foreach my $tmp_mp( @tmp_mp_files ) {
|
|
unlink $tmp_mp if $tmp_mp;
|
|
}
|
|
# delete profile temp file(s)
|
|
foreach my $tmp_prf ( map { $_->[7] } @MotifInfo ) {
|
|
unlink $tmp_prf if $tmp_prf;
|
|
}
|
|
}
|