Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Parsing Blast xml using Perl (BioPerl)

    Hi people,
    I'd like to resolve the following problem(s):

    I'm trying to resolve the structure of reads, partial of retroviral origin, partial of host genomic origin. In addition, I'd like to detect primer binding sites for a certain primer pair.

    My approach was to blast my sequences against the virus, the genome and the primer. In theory I could construct one big database, but for the primers I used blastn, while megablast was much faster and appropriate for the other two databases.

    Eyeballing the results, I can "easily" determine virus,genome and primer. Now I tried to parse my blast results to filter out overlapping hit (the genome has multiple insertions of the virus which will match every viral structure and I only want to extract genome hits without any overlap to any other db ref).

    important would be, that a genomic hit overlapping with any other hits will be sorted out (marked with x)

    Right at the beginning, the following problem came up:
    Parsing blast xml with BioPerls SearchIO gives me only one query all the time! Here is the test code I used:

    Code:
    #! /usr/bin/perl -w
    use strict;
    use Bio::Perl;
    use Bio::SearchIO;
    
    my $hitcountMAN=0;
    my $test;
    
      # Get the report 
    my  $searchio = Bio::SearchIO  ->new (
                                            -format => 'blastxml',
                                            -file   => $ARGV[0]);
    
    while(my $result = $searchio->next_result){
            my  $algorithm_type =  $result->algorithm;
            my $algorithm_version = $result ->algorithm_version;
    
            while (my $hit = $result->next_hit) {
                    my  $sseqid = $hit->name ;
                    my $qseqid= $result->query_description;
    #               print "$qseqid\t$sseqid\n";
    
    #initialize test
    
                    if($hitcountMAN==0){
                            $test=$qseqid;
                    }
    
                    if($qseqid!~$test){
                            print "new qseqid: $test\n";
                            $test=$qseqid;
                    }
                    $hitcountMAN++;
            }
    # Lets do some statistics
            my $resultcount=$searchio->result_count();
            my $hitcount=$result->num_hits;
    
            print "SearchIO parsed result from $ARGV[0]: $algorithm_type\t $algorithm_version\t $resultcount result(s) with $hitcount hit(s)\n";
    }
    
    open (INFILE, $ARGV[0]) or die $!;
    my @blastxml=<INFILE>;
    close INFILE;
    I would be very glad if anyone could tell me where I went into the wrong direction. The script should work with any blastxml output.

    If you have other suggestions solving the general problem, please tell me.

    Thanks in advance!

  • #2
    Update

    Okay, I gave up on bioperl and wrapped my tabular output.
    I got a multidimensional hash structure now.
    How may I compare all qstart/qend values of %result{$consensus} and keep only the once which do not overlap or overlap more then a given tolerance?
    $result{$qseqid}{$count}{'qstart'}
    $result{$qseqid}{$count}{'qend'}

    Any suggestions?

    Code:
    #! /usr/bin/perl -w
    use strict;
    
    
    #input -outfmt "6 qseqid sseqid qlen qstart qend slen sstart send length pident evalue"
    open (INPUT, $ARGV[0]) or die $!;
    my @blastin=<INPUT>;
    close INPUT;
    
    my $lqseqid="";
    my %result;
    my $count=0;
    my @primres;
    my @primer=("LTR","pol");
    my @virus=("KoRV");
    my %subres;
    
    foreach my $line(@blastin){
            chomp $line;
            my @array=split "\t",$line;
            my $qseqid=$array[0];
    #       my $sseqid=$array[1];
    #       my $qlen=$array[2];
    #       my $qstart=$array[3];
    #       my $qend=$array[4];
    #       my $slen=$array[5];
    #       my $sstart=$array[6];
    #       my $send=$array[7];
    #       my $length=$array[8];
            if($qseqid=~$lqseqid){
                    $result{$qseqid}{$count}{'qseqid'}=$array[0];
                    $result{$qseqid}{$count}{'sseqid'}=$array[1];
                    $result{$qseqid}{$count}{'qlen'}=$array[2];
    #For later comparison, it's important, that the smaller position always is the start position
                    if($array[3]<$array[4]){
                            $result{$qseqid}{$count}{'qstart'}=$array[3];
                            $result{$qseqid}{$count}{'qend'}=$array[4];
                    }
                    else{
                            $result{$qseqid}{$count}{'qstart'}=$array[4];
                            $result{$qseqid}{$count}{'qend'}=$array[3];
                    }
                    $result{$qseqid}{$count}{'slen'}=$array[5];
                    $result{$qseqid}{$count}{'sstart'}=$array[6];
                    $result{$qseqid}{$count}{'send'}=$array[7];
                    $result{$qseqid}{$count}{'length'}=$array[8];
     $count++;
            }
            else{
                    $lqseqid=$qseqid;
                    $count=0;
            }
    } 
    #print "Key: $_ and Value: $result{$_}\n" foreach(keys%result);
    #access toplevel (qseqid)
    foreach my $consensus(keys%result){
                    my $hitnum=keys%{$result{$consensus}};
    #               print "$consensus\t$hitnum\n";          #print the number of hits for each query
            foreach my $hit(keys %{$result{$consensus}}){
    
    #check primers depending on lenth ratio (alignment length/subject length) and sseqid
                    my $lengthratio=0.9;
                    for (@primer){
                            if($result{$consensus}{$hit}{length}/$result{$consensus}{$hit}{slen}>=$lengthratio){
    
                                    if($result{$consensus}{$hit}{sseqid} =~ $_ ){           #compares sseqid to all elements of @primer match TRUE
    #                                       push @primres, "$result{$consensus}{$hit}{qseqid}\t$result{$consensus}{$hit}{sseqid}\t$result{$consensus}{$hit}{qstart}\t$result{$consensus}{$hit}{qend}\t$result{$$
    #                                       print "$result{$consensus}{$hit}";
                                    }
                            }
    # other subjects (not in primer list)                                 
                            else{                    #compares sseqid to all elements of @primer match FALSE
                                    foreach(@virus){
                                            if($result{$consensus}{$hit}{sseqid} =~ $_){
                                                    # search non overlapping significant hits
    
                                            }
                                            else{
                                                    # search uncovered (non viral) regions
                                            }
                                    }
                            }
                    }
    
            }
    }

    Comment

    Latest Articles

    Collapse

    • seqadmin
      Strategies for Sequencing Challenging Samples
      by seqadmin


      Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
      03-22-2024, 06:39 AM
    • seqadmin
      Techniques and Challenges in Conservation Genomics
      by seqadmin



      The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

      Avian Conservation
      Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
      03-08-2024, 10:41 AM

    ad_right_rmr

    Collapse

    News

    Collapse

    Topics Statistics Last Post
    Started by seqadmin, Yesterday, 06:37 PM
    0 responses
    8 views
    0 likes
    Last Post seqadmin  
    Started by seqadmin, Yesterday, 06:07 PM
    0 responses
    8 views
    0 likes
    Last Post seqadmin  
    Started by seqadmin, 03-22-2024, 10:03 AM
    0 responses
    49 views
    0 likes
    Last Post seqadmin  
    Started by seqadmin, 03-21-2024, 07:32 AM
    0 responses
    66 views
    0 likes
    Last Post seqadmin  
    Working...
    X