restructured directory #noissue

arachne/generate_bogus_qual_for_tobacco.pl

@@ -0,0 +1,50 @@
+#!/usr/bin/perl
+
+use Carp;
+use strict;
+use warnings;
+use Bio::DB::Fasta;
+use Cwd qw{ abs_path };
+use File::Temp;
+
+my $fname = shift @ARGV;
+
+$fname = abs_path($fname);
+
+my $out_fname = $fname.".qual";
+
+#open (OUT, ">$out_fname") || confess "Failed to open '$out_fname' for writing";
+
+my $db = new Bio::DB::Fasta($fname);
+my @ids = $db->get_all_ids();
+
+foreach my $id(@ids) {
+
+    my $seq_len = $db->length($id);
+    my $defline = $db->header($id);
+
+    my ($gb_id, $trace_name, $lib_id);
+
+    if ($defline =~ /^(\S+)\s+(\S+)\s+(\S+)/) {
+        ($gb_id, $trace_name, $lib_id) = ($1, $2, $3);
+    } else {
+        confess "Failed to match defline '$defline'";
+    }
+
+    $trace_name =~ /^((.*)(\d{3})x([a-z]\d+))([a-z])\d+\.ab1/ || confess "Failed to match on trace name '$trace_name'";
+    my $template_id = $1;
+#    my $plate_id = $2;
+    my $plate_id = $2.$3;
+    my $well_id = $4;
+    $well_id =~ tr/a-z/A-Z/;
+    my $direction = $5;
+    $direction =~ tr/a-z/A-Z/;
+
+    print ">$defline\n";
+
+    my @quals = ('15') x $seq_len;
+
+    while (my @q_line = splice(@quals, 0, 19)) {
+        print join(" ", @q_line)."\n";
+    }
+}

arachne/generate_bogus_xml_for_tobacco.pl

@@ -0,0 +1,64 @@
+#!/usr/bin/perl
+
+use Carp;
+use strict;
+use warnings;
+use Bio::DB::Fasta;
+use Cwd qw{ abs_path };
+
+my $fname = shift @ARGV;
+
+$fname = abs_path($fname);
+
+my $out_fname = $fname.".xml";
+
+open (OUT, ">$out_fname") || confess "Failed to open '$out_fname' for writing";
+
+my $db = new Bio::DB::Fasta($fname);
+my @ids = $db->get_all_ids();
+
+print OUT <<XML;
+<?xml version="1.0"?>
+<trace_volume>
+XML
+
+foreach my $id(@ids) {
+
+    my $seq_len = $db->length($id);
+    my $defline = $db->header($id);
+
+    my ($gb_id, $trace_name, $lib_id);
+
+    if ($defline =~ /^(\S+)\s+(\S+)\s+(\S+)/) {
+        ($gb_id, $trace_name, $lib_id) = ($1, $2, $3);
+    } else {
+        confess "Failed to match defline '$defline'";
+    }
+
+    $trace_name =~ /^((.*)(\d{3})x([a-z]\d+))([a-z])\d+\.ab1/ || confess "Failed to match on trace name '$trace_name'";
+    my $template_id = $1;
+#    my $plate_id = $2;
+    my $plate_id = $2.$3;
+    my $well_id = $4;
+    $well_id =~ tr/a-z/A-Z/;
+    my $direction = $5;
+    $direction =~ tr/a-z/A-Z/;
+
+print OUT <<XML;
+    <trace>
+       <trace_name>$gb_id</trace_name>
+       <type>paired_production</type>
+       <library_id>$lib_id</library_id>
+       <plate_id>$plate_id</plate_id>
+       <well_id>$well_id</well_id>
+       <template_id>$template_id</template_id>
+       <trace_end>$direction</trace_end>
+       <clip_vector_left>1</clip_vector_left>
+       <clip_vector_right>$seq_len</clip_vector_right>
+       <insert_size>2500</insert_size>
+       <insert_stdev>500</insert_stdev>
+    </trace>
+XML
+
+}
+print OUT "</trace_volume>\n";

blast/blast_output_to_simple_table.pl

@@ -0,0 +1,51 @@
+#!/usr/bin/perl
+
+## quick script to do tblastx with cutoffs for percent identity and percent coverage
+## of the query sequence for all HSPs against a given subject sequence
+
+use warnings;
+use strict;
+
+use lib "/home/whitty/SVN/lib";
+
+use MyIO;
+
+use Data::Dumper;
+use Bio::SearchIO;
+use File::Basename;
+
+my $result_file = shift @ARGV;
+my $evalue_cutoff = shift @ARGV;
+
+unless ($evalue_cutoff) { $evalue_cutoff = '1e-5';}
+
+    my $infh = get_infh($result_file);
+
+    my $blast_report = new Bio::SearchIO(
+                                            -format =>  'blast',
+                                            -fh   =>    $infh,
+                                        );
+
+    while( my $result = $blast_report->next_result ) {
+
+        my $query_name = $result->query_name;
+        my $hit_count = 0;
+
+        while( my $hit = $result->next_hit ) {
+
+            my $hit_name = $hit->name;
+            my $rank = $hit->rank;
+            my $score = $hit->score;
+            my $sig = $hit->significance;
+
+                print join("\t", (
+                                    $query_name,
+                                    $hit_name,
+                                    $rank,
+                                    $score,
+                                    $sig,
+                                 ))."\n";
+            }
+
+    }
+

blast/blast_output_to_table.pl

@@ -0,0 +1,83 @@
+#!/opt/rocks/bin/perl
+
+## quick script to do tblastx with cutoffs for percent identity and percent coverage
+## of the query sequence for all HSPs against a given subject sequence
+
+use warnings;
+use strict;
+
+use Data::Dumper;
+use Bio::SearchIO;
+use File::Basename;
+
+my $result_file = shift @ARGV;
+my $identity_cutoff = shift @ARGV;
+my $coverage_cutoff = shift @ARGV;
+my $evalue_cutoff = shift @ARGV;
+
+unless ($identity_cutoff) { $identity_cutoff = 1;}
+unless ($coverage_cutoff) { $coverage_cutoff = 1;}
+unless ($evalue_cutoff) { $evalue_cutoff = '10';}
+
+
+    my $blast_report = new Bio::SearchIO(
+                                            -format =>  'blast',
+                                            -file   =>  $result_file,
+                                        );
+
+    while( my $result = $blast_report->next_result ) {
+
+        my $query_name = $result->query_name;
+        my $hit_count = 0;
+
+        while( my $hit = $result->next_hit ) {
+
+            my $hit_name = $hit->name;
+            my $rank = $hit->rank;
+            my $score = $hit->score;
+            my $sig = $hit->significance;
+            my $coverage = sprintf("%.2f", $hit->frac_aligned_query());
+            my $percent_coverage = $coverage * 100;
+
+            my $total_length = 0;
+            my $sum = 0;
+            my $simsum = 0;
+            while( my $hsp = $hit->next_hsp ) {
+                ## calculate average % identity across all hsps
+                $total_length += $hsp->length('query');
+                $sum += $hsp->frac_identical * $hsp->length('query');
+                $simsum += $hsp->frac_conserved * $hsp->length('query');
+            }
+            my $identity = sprintf("%.2f", $sum/$total_length);
+            my $similarity = sprintf("%.2f", $simsum/$total_length);
+            my $percent_identity = $identity * 100;
+            my $percent_similarity = $similarity * 100;
+
+            if ($percent_coverage >= $coverage_cutoff && $percent_identity >= $identity_cutoff) {
+
+                $hit_count++;
+
+                print join("\t", (
+                                    $query_name,
+                                    $hit_name,
+                                    $rank,
+                                    $score,
+                                    $sig,
+                                    $identity,
+                                    $similarity,
+                                    $coverage,
+                                    #$percent_identity.'%',
+                                    #$percent_coverage.'%',
+                                 ))."\n";
+            }
+
+        }
+
+        if ($hit_count) {
+            #print STDOUT $query_name."\n"; 
+        } else {
+            #print STDERR $query_name."\n";
+            print $query_name."\t-\n";
+        }
+
+    }