Changeset 83

Show
Ignore:
Timestamp:
16-07-12 15:35:06 (22 months 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 modified

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>;