Home > Archive > PERL CGI Beginners > December 2005 > am i right?
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
|
|
|
| I have Bioinformatics Related question and I am thinking on it from 4-5
days but I couldn't figure it out . I post this Problam in some Groups
too but i didn't get satisfactory answer may be there no one is
related to BioInformatics .
Is someone in this groups or somebuddy know that which group i should
use so please inform me.
Thanks
| |
| bahhab@hotmail.com 2005-11-23, 3:55 am |
| I guess the place to post depends on the question and with
bioinformatics there are lots of questions! I'm in a bioinformatics
lab with four others so you could try us with the question. If we
haven't got an answer we may know where you should ask.
| |
|
|
> I'm in a bioinformatics
> lab with four others so you could try us with the question.
Thanks,
I am making a programm to find a Gap Sequences for every sequence of
file after blasting . I made program for 2 sequence Which has 1 Gap
and now I want to generlaize my program But in this case may be some
sequence has no gap some has 2 gap some has 3 so i couldn't figureout
that how i can do this .i am posting my code for 2 sequence.
#!/usr/bin/perl
# A Program to take 2 sequences ,Blast it and find the gap into output
#sequence
use strict;
use warnings;
use bio::SeqIO;
use Bio::Seq;
use Bio::searchIo;
#Get First Sequence :
my$seq1=" ACTGGGGGGGATATTAACTAGCAAAAAACTCCTTCATTGA
AGCTTAATCAGAACTATCAAAGTTTGAGTTCATCAATCAA
ATAAAAAGAGGTGATACATTTCGCTACTATTCTAACACAA
CAAAATCGGATTTTTCACACGTCCGCGTGGCTTCTTCTTC
TTGTCTTGCTTCGCCACTTTGGGAGTTTGTCCTCTCACTT
TCCCCCCCCGCA";
#write in temprary file:
my $seqio = Bio::SeqIO ->new(-file =>'>temp1',-format =>'fasta');
$seq1= Bio::Seq -> new (-seq =>$seq1, -display_id =>'test1');
$seqio ->write_seq($seq1);
#Get Second Sequence :
my$seq2=" AACACACTCTTGTGTAGAAGAAACTCTTACACAATTGCCA
GTGAGAAAGAAAAATGTGGTTTGATATTAACTAGCAAAAA
ACTCCTTCATTGAAGCTTAATCAGAACTATCAAAGTTTGA
GTTCATCAATCAAATAAAAAGAGGTGATACATTTCGCAAA
CTAAGAAAAACCATCATCTTTTCTTGTCTTGCTTCGCCAC
TTTGGGAGTTTGTCCTCTCACTTTCCCCGCACGTGCGAGA
GATCCG
TGAACCTTACCCATTGTGATGAATTACTGCTTCTCTCACT
CTCTCTCCCTCTGTGTTTCCTGGTGCCGAATTC";
#write in temprary file:
$seqio = Bio::SeqIO ->new(-file =>'>temp2',-format =>'fasta');
$seq2= Bio::Seq -> new (-seq =>$seq2, -display_id =>'test2');
$seqio ->write_seq($seq2);
#Blast to both sequences and store output in file name temp.blast:
system ("C:/downloads/bin/bl2seq -i temp1 -j temp2 -p blastn -o
temp.blast");
my $in = Bio::SearchIO ->new(-format => 'blast',
-file => 'temp.blast');
#Print Starting and Ending Point of 1 Hsps-
my $result=$in->next_result;
my $hit =$result->next_hit;
my $hsp1 =$hit->next_hsp;
my$start1=$hsp1->query->start;
my $end1=$hsp1->query->end;
print "start First Hsp ",$start1,"\n";
print "End First Hsp ",$end1,"\n\n";
#Print Starting and Ending Point of 2 Hsps-
my $hsp2 =$hit->next_hsp;
my$start2=$hsp2->query->start;
my $end2=$hsp2->query->end;
print "start Second Hsp ",$start2,"\n";
print "End Second Hsp ",$end2,"\n\n";
#Print Gap Sequence
print "Gap Sequence- ",$seq1->subseq($end1+1,$start2-1);
Do you have any Idea about it? People think in all groups that i am
rude but don't mind if you feel that I write something rudely bcoz i
dont know too much about English so Please...
| |
| bahhab@hotmail.com 2005-11-24, 3:55 am |
| At least you are trying with another language, I never got further than
english.
If I understand correctly you want to put in a loop that reads the
start and end of each hsp then reports the sequence between the start
of the current hsp and the end of the pervious one.
I haven't tried running this so there could be a typo in it. Give it a
go and come back if you've still got issues.
my $in = Bio::SearchIO ->new(-format => 'blast',
-file => 'temp.blast');
my $result=$in->next_result;
my $hit =$result->next_hit;
#define start and end of hsps 1 & 2 as zero
my $start1= my $end1= 0;
my $start2= my $end2= 0;
#if you want to count the hsps...
my $hsp_count=0;
#go throught each hsp in turn
while (my $hsp = $hit->next_hsp) {
$hsp_count++; #add 1 to the count
#set start and end of 1 to the values of the previous hsp
# that will be 0 if it is first hsp
$start1=$start2;
$end1=$end2;
#set start and end of 2 to the values current hsp
$start2=$hsp->query->start;
$end2=$hsp->query->end;
#unless this is the first hsp
# (you need to be on the 2nd hsp to have a gap)
unless ($start1==0) {
#Print Starting and Ending Point of 1 Hsps-
print "start previous Hsp ",$start1,"\n";
print "End previous Hsp ",$end1,"\n\n";
#Print Starting and Ending Point of 2 Hsps-
print "start $hsp_count Hsp ",$start2,"\n";
print "End $hsp_count Hsp ",$end2,"\n\n";
#Print Gap Sequence
print "Gap Sequence- ",$seq1->subseq($end1+1,$start2-1);
}
}
| |
|
|
> I haven't tried running this so there could be a typo in it. Give it a
> go and come back if you've still got issues.
> #Print Starting and Ending Point of 2 Hsps-
> print "start $hsp_count Hsp ",$start2,"\n";
> print "End $hsp_count Hsp ",$end2,"\n\n";
>
> #Print Gap Sequence
> print "Gap Sequence- ",$seq1->subseq($end1+1,$start2-1);
Thanks ,Your logic is right and my programm is also working but i think
i am doing something wrong because i am using subseq function to get
Gap Sequence.
which is depend on hsp value .
and in blast programm Hsp value is depend on identities or something
but thats why it gives me valus like
start previous Hsp 461
End previous Hsp 703
start 2 Hsp 146
End 2 Hsp 243
start previous Hsp 146
End previous Hsp 243
start 3 Hsp 332
End 3 Hsp 363
start previous Hsp 332
End previous Hsp 363
start 4 Hsp 14
End 4 Hsp 47
start previous Hsp 14
End previous Hsp 47
start 5 Hsp 613
End 5 Hsp 624
start previous Hsp 613
End previous Hsp 624
start 6 Hsp 565
End 6 Hsp 575
so it can not give me Gap because grater value is coming first .what
you think about it ?
Is i have to use some another function or I can sort out this values ?
| |
|
|
> start previous Hsp 461
> End previous Hsp 703
>
> start 2 Hsp 146
> End 2 Hsp 243
>
> start previous Hsp 146
> End previous Hsp 243
>
> start 3 Hsp 332
> End 3 Hsp 363
>
> start previous Hsp 332
> End previous Hsp 363
>
> start 4 Hsp 14
> End 4 Hsp 47
>
> start previous Hsp 14
> End previous Hsp 47
>
> start 5 Hsp 613
> End 5 Hsp 624
>
> start previous Hsp 613
> End previous Hsp 624
>
> start 6 Hsp 565
> End 6 Hsp 575
Okay i decided to Sort it .and i am doing
my $in = Bio::SearchIO ->new(-format => 'blast',
-file => 'temp.blast');
#define start and end of hsps 1 & 2 as zero
my $start1= my $end1= 0;
my $start2= my $end2= 0;
my ( @start, @end );
# count the hsps...
my $hsp_count=0;
while( my $result = $in->next_result ) {
while( my $hit = $result->next_hit ){
#go throught each hsp in turn
while (my $hsp = $hit->next_hsp) {
if( $hsp->length('total') > 15) {
if ( $hsp->percent_identity >= 75 ) {
$hsp_count++; #add 1 to the count
#set start and end of 1 to the values of the previous hsp
# that will be 0 if it is first hsp
$start1=$start2;
$end1=$end2;
#set start and end of 2 to the values current hsp
$start2=$hsp->hit->start;
$end2=$hsp->hit->end;
# it will store all the Starting and Ending values of Hsp in Array:
push @start,$start2;
push @end,$start2;
#Sort that values
my @start3= sort{$a <=> $b} @start;
my @end3= sort{$a <=> $b} @end;
#unless this is the first hsp
# (we need to be on the 2nd hsp to have a gap)
unless ($start1==0) {
my $counter=0;
#Print Sorted Starting and Ending Point of 1 Hsps-
print "start previous Hsp ",$start1,"\n";
print "End previous Hsp ",$end1,"\n\n";
$start1=$start3[$counter];
$end1=$end3[$counter];
#Print Sorted Starting and Ending Point of 2 Hsps-
print "start $hsp_count Hsp ",$start3[$counter],"\n";
print "End $hsp_count Hsp ",$end3[$counter],"\n\n";
$counter++;
}
}
}
}
}
}
But it is not working. Can you Help me Please.
Thanks
| |
| bahhab@hotmail.com 2005-12-01, 6:55 pm |
| The only obvious typo I can see in your script is that you are pushing
$start2 into both the start and end array.
However I think your theory is a bit wrong. You need to collect all
the hsps for a hit and store start/end for each in a 2 dimentional
array. Then sort the array on the start values. Then loop through the
array and print the starts/ends in order. Hopefully this is what the
following script does.
When you get to the point of calculating the gaps you need to remember
that there could be hsps that overlap or even one hsp that is inside a
longer one. So you'll need an if statment along the lines of 'if
(current_hsp_start < previous_hsp_end ){ last;} #go straight to next
pass through the loop without calculating gap
Hope that helps.
while( my $result = $in->next_result ) {
while( my $hit = $result->next_hit ){
#go throught each hsp in turn
while (my $hsp = $hit->next_hsp) {
if( $hsp->length('total') > 15) {
if ( $hsp->percent_identity >= 75 ) {
$hsp_count++; #add 1 to the count
#set start and end of 1 to the values of the previous
hsp
# that will be 0 if it is first hsp
$start1=$start2;
#end1=$end2;
#set start and end of 2 to the values current hsp
$start2=$hsp->hit->start;
$end2=$hsp->hit->end;
# it will store all the Starting and Ending values of
Hsp in Array:
$array[0][$hsp_count]=$start2;
$array[1][$hsp_count]=$end2;
} #do that for all hsps
}
}
@sorted = sort {$a->[0] <=> $b->[0]} @array; #sort array
for ($array_element=0; $array_element<scalar(@array);
$array_element++){
#(we need to be on the 2nd hsp to have a gap)
unless ($array_element==0) {
#Print Sorted Starting and Ending Point of 1 Hsps-
print "start previous Hsp ",$array[0][$array_element-1],"\n";
print "End previous Hsp ",$array[1][$array_element-1],"\n\n";
#Print Sorted Starting and Ending Point of 2 Hsps-
print "start $array_element Hsp
",$array[0][$array_element],"\n";
print "End $array_element Hsp
",$array[1][$array_element],"\n\n";
}
}
}
}
Rita wrote:
>
> Okay i decided to Sort it .and i am doing
> my $in = Bio::SearchIO ->new(-format => 'blast',
> -file => 'temp.blast');
> #define start and end of hsps 1 & 2 as zero
> my $start1= my $end1= 0;
> my $start2= my $end2= 0;
> my ( @start, @end );
> # count the hsps...
> my $hsp_count=0;
> while( my $result = $in->next_result ) {
> while( my $hit = $result->next_hit ){
>
> #go throught each hsp in turn
> while (my $hsp = $hit->next_hsp) {
> if( $hsp->length('total') > 15) {
> if ( $hsp->percent_identity >= 75 ) {
> $hsp_count++; #add 1 to the count
> #set start and end of 1 to the values of the previous hsp
> # that will be 0 if it is first hsp
> $start1=$start2;
> $end1=$end2;
> #set start and end of 2 to the values current hsp
> $start2=$hsp->hit->start;
> $end2=$hsp->hit->end;
> # it will store all the Starting and Ending values of Hsp in Array:
> push @start,$start2;
> push @end,$start2;
> #Sort that values
> my @start3= sort{$a <=> $b} @start;
> my @end3= sort{$a <=> $b} @end;
> #unless this is the first hsp
> # (we need to be on the 2nd hsp to have a gap)
> unless ($start1==0) {
> my $counter=0;
> #Print Sorted Starting and Ending Point of 1 Hsps-
> print "start previous Hsp ",$start1,"\n";
> print "End previous Hsp ",$end1,"\n\n";
> $start1=$start3[$counter];
> $end1=$end3[$counter];
>
> #Print Sorted Starting and Ending Point of 2 Hsps-
> print "start $hsp_count Hsp ",$start3[$counter],"\n";
> print "End $hsp_count Hsp ",$end3[$counter],"\n\n";
> $counter++;
> }
> }
> }
> }
> }
> }
> But it is not working. Can you Help me Please.
> Thanks
| |
|
|
> while( my $result = $in->next_result ) {
> while( my $hit = $result->next_hit ){
>
> #go throught each hsp in turn
> while (my $hsp = $hit->next_hsp) {
> if( $hsp->length('total') > 15) {
> if ( $hsp->percent_identity >= 75 ) {
> $hsp_count++; #add 1 to the count
> #set start and end of 1 to the values of the previous
> hsp
> # that will be 0 if it is first hsp
> $start1=$start2;
> #end1=$end2;
> #set start and end of 2 to the values current hsp
> $start2=$hsp->hit->start;
> $end2=$hsp->hit->end;
> # it will store all the Starting and Ending values of
> Hsp in Array:
> $array[0][$hsp_count]=$start2;
> $array[1][$hsp_count]=$end2;
This line is giving me syntax error near "$array[".
|
|
|
|
|