Changeset 83


Ignore:
Timestamp:
Jul 16, 2012, 3:35:06 PM (4 years ago)
Author:
victor.de.jager@…
Message:

Fixed bug in unique data recovery
Fixed bug in data formatting (getting rid of slashes/backslashes)
Fixed bug in sam2rovar data

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/rovar.pl

    r81 r83  
    4343                                       
    4444$ini{'-q2'}{val}             = '';
    45 #$ini{'-q2'}{checkisexisting} = 1;
     45$ini{'-q2'}{checkisexisting} = 0;
    4646$ini{'-q2'}{checkisstring}   = 1;
    47 #$ini{'-q2'}{mandatory}       = 0;
     47$ini{'-q2'}{mandatory}       = 0;
    4848
    4949$ini{'-desc'}{desc} = 'Description of the sample; this will be the directory name under which the data will be saved.';
     
    167167my $prog_bowtie_build = $pathprog . 'aligners/bowtie/bowtie-build';
    168168
     169my $prog_bowtie2          = $pathprog . 'aligners/bowtie2/bowtie2';
     170my $prog_bowtie2_build= $pathprog . 'aligners/bowtie2/bowtie2-build';
     171
     172
    169173my $prog_soap             = $pathprog . 'aligners/soap/soap';
    170174my $prog_soap_build   = $pathprog . 'aligners/soap/2bwt-builder';
     
    173177my $prog_smalt            = $pathprog . 'aligners/smalt/smalt'; #create a symlink to the smalt binary for your platform
    174178
    175 my $prog_blast         = $pathprog . 'aligners/blast/bin/blastn';
    176 my $prog_makeblastdb   = $pathprog . 'aligners/blast/bin/makeblastdb';
     179my $prog_blast         = $pathprog . 'aligners/blast+/bin/blastn';
     180my $prog_makeblastdb   = $pathprog . 'aligners/blast+/bin/makeblastdb';
    177181
    178182#other auxillary programs
     
    428432                        $format = 'fasta';
    429433                }
    430                 if ($line1 =~ /[\\\/]/){
    431                         #replace all slashes and backslashes
    432                         print dotime() . " ".$format." replacing slashes and backslashes in query id's\n";
    433                         system("sed -i  '/[>]/ s/[/\//]/_/' $fname");
    434                 }
     434                #replace all slashes and backslashes
     435                print dotime() . " ".$format." replacing slashes and backslashes in query id's of $fname\n";
     436                system("sed -i  '/[>]/ s%[/\]%_%g' $fname");
    435437               
    436438        }
     
    439441                        $format = 'fastq';
    440442                }
    441                 if ($line1 =~ /[\\\/]/){
    442                         #replace all slashes and backslashes
    443                         print dotime() . " ".$format." replacing slashes and backslashes in query id's\n";
    444                         system("sed -i  '/[@+]/ s/[/\//]/_/' $fname");
    445                 }
     443                #replace all slashes and backslashes
     444                print dotime() . " ".$format." replacing slashes and backslashes in query id's of $fname\n";
     445                system("sed -i  '/[@+]/ s%[/\]%_%g' $fname");
     446               
    446447        }
    447448        print dotime() . " ".$format." format detected\n";
     
    665666        }
    666667        return $fname;
    667 }    # clean_fasta
     668}   # clean_fasta
    668669
    669670sub repmask {
     
    710711        # cleaning reference fasta file
    711712        $reference = clean_fasta($reference);
     713       
     714       
    712715        if ( lc( $ini{'-repmask'}{val} ) eq 'yes' or $ini{'-repmask'}{val} eq '1' )
    713716        {
     
    805808    my $t3    = shift();
    806809    my $q1    = $ini{'-q1'}{val};
    807     my $q2    = $ini{'-q2'}{val};
     810    my $q2    = $ini{'-q2'}{val} if exists($ini{'-q2'}{val});
    808811
    809812   
     
    830833    # choose which aligner to use.
    831834    if ($alnprog eq 'blat'){
     835        if ($ini{'-q2'}{val}){
     836                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     837        }
     838       
    832839        &run_blat_aligner($dir,$q1, $type);     
    833840    }
    834841    if ($alnprog eq 'blast'){
    835         &run_blast_aligner($dir,$q1, $type);   
     842        if ($ini{'-q2'}{val}){
     843                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     844        }
     845        &run_blast_aligner($dir,$q1, $type);   
    836846    }
    837847    if ($alnprog eq 'bwa' || $alnprog eq 'bwasw'){
     848        if ($ini{'-q2'}{val}){
     849                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     850        }
     851       
    838852        &run_bwa_aligner($dir,$q1, $type);   
    839853    }
    840854   
    841855    if ($alnprog eq 'ssaha2'){
     856        if ($ini{'-q2'}{val}){
     857                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     858        }
     859       
    842860        &run_ssaha2_aligner($dir,$q1, $type);   
    843861    }
    844862    if ($alnprog eq 'bowtie'){
    845         &run_bowtie_aligner($dir,$q1, $type);   
     863        if ($ini{'-q2'}{val}){
     864                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     865        }
     866        &run_bowtie_aligner($dir,$q1, $type);   
    846867    }
    847868   
    848869    if ($alnprog eq 'soap'){
     870        if ($ini{'-q2'}{val}){
     871                $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     872        }
     873       
    849874        &run_soap_aligner($dir,$q1, $type);   
    850875    }
     
    858883           
    859884    }
    860    
    861885   
    862886    if ($alnprog eq 'pass'){
     
    871895                print dotime() . " Concatenate Pass results $pass1 and $pass2 \nCommandline:\n";
    872896                system($cmd);
    873                
     897                       
     898                        $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
     899                                       
    874900                &blat2rovar($dir,$q1,$pass); #
    875901               
     
    11041130    my $q1    = shift();
    11051131    my $type = shift();
    1106     my $bwa_sai = $dir . 'tmp/bwa.sai';
     1132   
     1133    my $bwa1_sai = $dir . 'tmp/bwa1.sai';
     1134   
    11071135    my $bwa_sam  = $dir . 'tmp/bwa.sam';
    11081136    my $bwa_bam  = $dir . 'tmp/bwa.bam';
     
    11201148                                                    'REF' => $reference,
    11211149                                                    'QUERY' => $q1,
    1122                                                     'OUTPUT' => $bwa_sai);
     1150                                                    'OUTPUT' => $bwa1_sai);
    11231151            print dotime() . " Performing bwa search of first query to reference sequence\nCommandline:".$cmd."\n";
    11241152            system( $cmd );
    1125        
    1126             general::error("no bwa output for $q1") if ( !-e $bwa_sai );
     1153           
     1154            general::error("no bwa output for $q1") if ( !-e $bwa1_sai );
     1155           
    11271156            print dotime() . " Converting bwa output from $q1 to bam format\n";
    1128                 system( $prog_bwa . " samse " . $reference . " " . $bwa_sai . " " .$q1 ." -f ". $bwa_sam );
     1157                system( $prog_bwa . " samse " . $reference . " " . $bwa1_sai . " " .$q1 ." -f ". $bwa_sam );
    11291158        }
    11301159       
     
    12491278        print dotime() . " Skipping smalt search\n";
    12501279    }
     1280    $q1 = &concatenate_input_files($q1,$q2, './concatenated_data')
    12511281    &bam2rovar($dir, $q1, $smalt_bam);
    12521282} # smalt_aligner
     
    14331463} # run_blast_aligner
    14341464
     1465sub concatenate_input_files {
     1466        my $q1 = shift;
     1467        my $q2 = shift;
     1468        my $q1q2 = shift;
     1469       
     1470        if (! -e $q1q2){
     1471                my $cmd = "cat $q1 $q2 > $q1q2";
     1472                print dotime()
     1473                . " Combining both datafiles\n";
     1474                system($cmd);
     1475        }
     1476        return $q1q2;
     1477}
    14351478
    14361479sub summarize_results {
    14371480        my $dir = 'results';
    14381481        my $cmd;
    1439 
     1482        my $ref;
    14401483        general::makedir($dir);
    14411484        $cmd = $prog_cat . ' ' . $text_explanation . ' > explanation.txt';
     
    14471490        system($cmd);
    14481491        print dotime() . " Determining unique reference sequences\n";
    1449         $cmd = $prog_uniqueseq . ' -i=' . $ini{'-s'}{val} . '_clean_repmask -dir=' . $dir;
     1492       
     1493        if (-e $ini{'-s'}{val} . '_clean_repmask'){
     1494                $ref = $ini{'-s'}{val} . '_clean_repmask';
     1495        }
     1496       
     1497        if (-e $ini{'-s'}{val} . '_repmask'){
     1498                $ref = $ini{'-s'}{val} . '_repmask';
     1499        }
     1500               
     1501        $cmd = $prog_uniqueseq . ' -i=' . $ref . ' -dir=' . $dir;
    14501502        system($cmd);
    14511503        print dotime()
  • trunk/sam2filter.pl

    r73 r83  
    11#!/usr/bin/perl
    2 #use Data::Dumper::Simple;
     2use Data::Dumper::Simple;
    33#
    44use strict ;
     
    3131$ini{'-i1'}{checkisstring} = 1 ;
    3232$ini{'-i1'}{mandatory} = 1 ;
     33
     34#$ini{'-i2'}{desc} = 'second query sequence file.' ;
     35#$ini{'-i2'}{val}  = '' ;
     36#$ini{'-i2'}{checkisexisting} = 1 ;
     37#$ini{'-i2'}{checkisstring} = 1 ;
     38#$ini{'-i2'}{mandatory} = 0 ;
    3339
    3440$ini{'-o'}{desc} = 'Filtered ROAST file, used for input for snp indel classification of ROAST' ;
     
    115121} # readfasta
    116122
    117 sub readquery2 {
     123sub readqueries {
     124       
     125        if ($ini{'-i1'}{val}) {
     126                readquery($ini{'-i1'}{val});
     127        }
     128        if ($ini{'-i2'}{val}) {
     129                readquery($ini{'-i2'}{val});
     130        }
     131}
     132
     133sub readquery {
    118134    my $line ;
    119135    my $name ;
     
    121137    my $fasta ;
    122138    my $idtag ;
    123    
    124     print dotime()." Reading query FASTA/FASTQ file $ini{'-i1'}{val}\n" ;
    125     open(FILE,$ini{'-i1'}{val}) or general::error('could not open '.$ini{'-i1'}{val}) ;
     139    my $infile = shift;
     140   
     141    print dotime()." Reading query FASTA/FASTQ file $infile\n" ;
     142    open(FILE,$infile) or general::error('could not open '.$infile) ;
    126143   
    127144    #get first line
     
    134151       $fasta=1;
    135152    }
    136    
    137    
     153       
    138154    while (<FILE>) {
    139155                $line = $_ ;
     
    146162                        }
    147163                        $fasta = 1;
     164                   
    148165                    #print "name: ".$name."\n";
     166                   
    149167                    next();
    150168                } else {
     
    319337                $subject_end    = $a->end;
    320338               
    321                 print $line." start position subject larger than stop position\n" if ($subject_end < $subject_start) ;
     339                # print $line." start position subject larger than stop position\n" if ($subject_end < $subject_start) ;
    322340                # where does the alignment start in the query sequence
    323341                $query = $a->query->seq_id;
     
    334352                $query_end   = $a->query->end;
    335353               
    336                 #print "$subject $subject_start $subject_start $query $query_start $query_end  $queries{$query}{len} \n";
    337                        
    338                 if (abs($query_start - $query_end) + 1 > $queries{$query}{len} ) {
    339                                 $queries{$query}{numhits}++ ;
    340                 }
     354                    if (exists($queries{$query})){
     355                        if (abs($query_start - $query_end) + 1 > $queries{$query}{len} ) {
     356                                        $queries{$query}{numhits}++ ;
     357                        }       
     358                   
     359                    }else{
     360                        print "$query not found in sequence file\n";
     361                        print "$subject $subject_start $subject_start $query $query_start $query_end  $queries{$query}{len} \n\n";
     362                        #exit();
     363                    }
     364               
    341365                       
    342366                        #my @scores    = $a->qscore;     # per-base quality scores
     
    451475&parseparams() ;
    452476&readfasta() ;
    453 &readquery2() ;
     477&readqueries() ;
    454478
    455479#system ("ps aux | grep axt2filter | grep -v grep") ; print "wait read\n" ; my $name = <STDIN>;
Note: See TracChangeset for help on using the changeset viewer.