--- /dev/null
+#!/usr/bin/perl
+
+# hack to pull the alignments and profiles from psiblast
+
+# -m6 switch is needed to get the alignments
+# in the BLAST report, also generates the
+# psiblast frequency profiles, as they
+# are no longer present in the output
+
+use Getopt::Std;
+getopt(hfgud);
+
+if ($opt_f eq"req"){$freq=1}
+elsif ($opt_g eq "ap"){$gap=1}
+elsif ($opt_u eq "ngap"){$ungap=1}
+elsif ($opt_d eq "egap"){$degap=1}
+elsif ($opt_h eq "elp"){
+ print "
+# Parse_PSI
+
+This little script allows parsing of psiblast report files
+to extract the embedded sequence alignment, and write out
+in different formats.
+
+parse_psi -help : help
+parse_psi -degap : remove all gaps in alignment and print as
+ a fasta format file
+parse_psi -ungap : remove gaps in query sequence, and those
+ below the gap in the query sequence, print
+ as fasta format
+parse_psi -gap : print alignment with gaps in fasta format
+parse_psi -freq : print frequency profile needed for Jnet
+
+Usage:
+
+parse_psi -freq file.blast > file.freq
+\n";exit(0);
+}
+
+else {print "NO options given, try -help\n";exit}
+if (!-e "$ARGV[0]"){
+ print ("No such file as $ARGV[0] exists in this path...\n");
+ exit (0);
+}
+
+if ($ARGV[0] eq ""){
+ print "No filename given... \n";
+ exit;
+}
+
+open IN, "$ARGV[0]";
+$num=$thisone=$c=0;
+while(<IN>){
+ if (/^Sequences producing significant alignments/){
+ $num++;
+ }
+}
+print STDERR "$num alignments were found in PSIBlast file $ARGV[0]\n";
+$num--;
+
+close IN;
+open IN, "$ARGV[0]";
+while(<IN>){
+ if (/^Sequences producing significant alignments/){
+ if ($thisone == $num){
+ while (($ln2 = <IN>)){
+ if ($ln2 =~ /^QUERY/){
+ $alignout[$c]=$ln2;
+ $c++;
+ while (($align = <IN>) !~ /^ Database/){
+ $alignout[$c]=$align;
+ $c++;
+ }
+ last;
+ }
+ }
+ }
+ $thisone++;
+ }
+}
+
+$c=$j=$z=$blocks=0;
+
+# get number of blocks and number of sequences in a block
+
+foreach (@alignout){
+ if ($_ =~ /^\n/){
+ $blocks++;
+ }
+ if ($_ !~ /^\n/){
+ $store[$blocks][$z]=$_;
+ $z++;
+ }
+}
+$npb=($z/($blocks+1));
+
+for ($i=0; $i < ($blocks+1); $i++){
+ for ($j=0; $j < $npb ; $j++){
+ # get the alignment itself
+ #print "$store[$i][$q]";
+ (@tokens)=split(" ",$store[$i][$q]);
+ #print "$tokens[0] $tokens[1] $tokens[2] $tokens[3]\n";
+ if($#tokens == 3){
+ $id=$tokens[0];
+ $interest=$tokens[2];
+ }
+ if($#tokens == 1){
+ $id=$tokens[0];
+ $interest=$tokens[1];
+ }
+ if($#tokens == 2){
+ $id=$tokens[0];
+ $interest=$tokens[1];
+ if ($interest =~ /\d/g){
+ # get the next token, number
+ # must be at this end
+ $interest=$tokens[2];
+ }
+ }
+ # remove any digits that may creep in with the 3 token,
+ # this would mean the coord is on the end not in the
+ # middle.
+ $interest =~ s/\d//g;
+ undef @tokens;
+
+ $namekeep[$j]=$id;
+ $seqkeep[$j]=$seqkeep[$j].$interest;
+ $q++;
+ }
+}
+
+
+if ($gap==1 || $degap==1){
+ $c=0;
+ foreach (@seqkeep){
+
+ if ($degap==1){
+ $seqkeep[$c] =~ s/-//g;
+ }
+ print ">$namekeep[$c]\n$seqkeep[$c]\n";
+ $c++;
+ }
+}
+
+if ($ungap==1 || $freq == 1){
+ $c=0;
+
+ @sequence=split("",$seqkeep[0]);
+
+ foreach (@sequence){
+ if($sequence[$count] eq "\-"){
+ $printer[$count]=1;
+ }
+ else{$printer[$count]=0;}
+
+ $count++;
+ }
+ $d=0;
+ foreach $erm (@seqkeep){
+ undef @sequence;
+ @sequence=split("",$erm);
+ $count=0;
+ foreach (@sequence){
+ if ($printer[$count]==0){
+ $ung[$d]=$ung[$d].$sequence[$count];
+ }
+ $count++;
+ }
+ $d++;
+ }
+ if ($ungap == 1){
+ $c=0;
+ foreach $re (@ung){
+ print ">$namekeep[$c]\n$re\n";
+ $c++;
+ }
+ }
+ if ($freq == 1){
+ # do the counting....
+ $seqlen=length($ung[0]);
+
+ #get a double array
+ for ($i = 0; $i < $#ung; $i++){
+ undef @sequence;
+ @sequence=split("",$ung[$i]);
+
+ # convert sequence to an integer value
+ @ints = seq2int(@sequence);
+ $count=0;
+ $j=0;
+ foreach (@sequence) {
+ $ar[$i][$j]=$ints[$j];
+ $j++;
+ }
+ }
+
+ # zero countup
+ for ($r=0; $r <20; $r++){
+ $countup[$r]=0;
+ }
+
+ for ($i = 0; $i < $seqlen; $i++) {
+ # skip the query sequence
+ for ($j = 1; $j < $#ung ; $j++){
+ # look at array
+ $countup[$ar[$j][$i]]++;
+ }
+
+ # get total counts for this column
+ $tot = 0;
+ for ($r = 0; $r < 20; $r++){
+ $tot += $countup[$r];
+ }
+ for ($r=0; $r <20; $r++){
+ # get percent and round up to 10
+ # this seems to be slightly different to how
+ # psiblast works god knowns how....
+ # works for now. :-(
+ $outnum= ($countup[$r] / $tot ) * 10;
+ printf ("%2.0f ",$outnum);
+ $countup[$r]=0;
+ }
+ print "\n";
+ }
+ }
+}
+
+sub seq2int{
+ my(@s)=@_;
+ $g=0;
+ for ($g=0; $g <= $#s; $g++){
+ if ($s[$g] eq 'A'){ $intseq[$g] = "0"};
+ if ($s[$g] eq 'R'){ $intseq[$g] = "1"};
+ if ($s[$g] eq 'N'){ $intseq[$g] = "2"};
+ if ($s[$g] eq 'D'){ $intseq[$g] = "3"};
+ if ($s[$g] eq 'C'){ $intseq[$g] = "4"};
+ if ($s[$g] eq 'Q'){ $intseq[$g] = "5"};
+ if ($s[$g] eq 'E'){ $intseq[$g] = "6"};
+ if ($s[$g] eq 'G'){ $intseq[$g] = "7"};
+ if ($s[$g] eq 'H'){ $intseq[$g] = "8"};
+ if ($s[$g] eq 'I'){ $intseq[$g] = "9"};
+ if ($s[$g] eq 'L'){ $intseq[$g] = "10"};
+ if ($s[$g] eq 'K'){ $intseq[$g] = "11"};
+ if ($s[$g] eq 'M'){ $intseq[$g] = "12"};
+ if ($s[$g] eq 'F'){ $intseq[$g] = "13"};
+ if ($s[$g] eq 'P'){ $intseq[$g] = "14"};
+ if ($s[$g] eq 'S'){ $intseq[$g] = "15"};
+ if ($s[$g] eq 'T'){ $intseq[$g] = "16"};
+ if ($s[$g] eq 'W'){ $intseq[$g] = "17"};
+ if ($s[$g] eq 'Y'){ $intseq[$g] = "18"};
+ if ($s[$g] eq 'V'){ $intseq[$g] = "19"};
+ if ($s[$g] eq '-'){ $intseq[$g] = "25"};
+ }
+ return (@intseq);
+}
git://source.jalview.org / jpred.git / blobdiff