Merge branch 'master' of https://github.com/PASApipeline/pasapipeline

PerlLib/Exons_to_geneobj.pm

@@ -27,9 +27,13 @@ sub create_gene_obj {
     ## exons_ref should be end5's keyed to end3's for all exons.
     my ($gene_struct_mod, $cdna_seq)  = &get_cdna_seq ($exons_href, $sequence_ref);
     
+    
+    
     my $cdna_seq_length = length $cdna_seq;
     my $long_orf_obj = new Longest_orf();
     
+    #print STDERR "CDNA_SEQ: [$cdna_seq], length: $cdna_seq_length\n";
+
     # establish long orf finding parameters.
     $long_orf_obj->forward_strand_only();
     if ($partial_info_href->{"5prime"}) {
@@ -43,24 +47,28 @@ sub create_gene_obj {
     $long_orf_obj->get_longest_orf($cdna_seq);
     my ($end5, $end3) = $long_orf_obj->get_end5_end3(); 
     
-    print "CDS: $end5, $end3\n" if $SEE;
+    #print STDERR "***   CDS: $end5, $end3\n";# if $SEE;
     my $gene_obj = &create_gene ($gene_struct_mod, $end5, $end3);
     
+    #print STDERR $gene_obj->toString();
+    
     $gene_obj->create_all_sequence_types($sequence_ref);
     my $protein = $gene_obj->get_protein_sequence();
-    
+    my $recons_cds = $gene_obj->get_CDS_sequence();
+    #print STDERR "reconsCDS: $recons_cds\n";
+
     ## check partiality
     if ($protein) { # it is possible that we won't have any cds structure
         if ($protein !~ /^M/) {
             # this would require that we allowed for 5prime partials
             unless ($partial_info_href->{"5prime"}) {
-                confess "Error, have 5' partial protein when 5prime partials weren't allowed!\n$protein\n";
+                confess "Error, have 5' partial protein when 5prime partials weren't allowed!\n$protein\n$cdna_seq\n";
             }
         }
         if ($protein !~ /\*$/) {
             # this would require that we allowed for 3prime partials
             unless ($partial_info_href->{"3prime"}) {
-                confess "Error, have 3' partial protein when 3prime partials weren't allowed!\n$protein\n";
+                confess "Error, have 3' partial protein when 3prime partials weren't allowed!\n$protein\n$cdna_seq\n";
             }
         }
    
@@ -80,6 +88,9 @@ sub create_gene_obj {
 ####
 sub get_cdna_seq {
     my ($gene_struct, $assembly_seq_ref) = @_;
+    
+    my $seq_length = length($$assembly_seq_ref);
+
     my (@end5s) = sort {$a<=>$b} keys %$gene_struct;
     my $strand = "?";
     foreach my $end5 (@end5s) {
@@ -90,7 +101,7 @@ sub get_cdna_seq {
     }
     if ($strand eq "?") {
         print Dumper ($gene_struct);
-        die "ERROR: I can't determine what orientation the cDNA is in!\n";
+        confess "ERROR: I can't determine what orientation the cDNA is in!\n";
     }
     print NOTES "strand: $strand\n";
     my $cdna_seq;
@@ -99,6 +110,11 @@ sub get_cdna_seq {
     foreach my $end5 (@end5s) {
         #print $end5;
         my $end3 = $gene_struct->{$end5};
+        
+        if ($end5 > $seq_length || $end3 > $seq_length) {
+            confess "Error, coords are out of bounds of sequence length: $seq_length:\n" . Dumper(\$gene_struct);
+        }
+        
         my ($coord1, $coord2) = sort {$a<=>$b} ($end5, $end3);
         my $exon_seq = substr ($$assembly_seq_ref, $coord1 - 1, ($coord2 - $coord1 + 1));
         $cdna_seq .= $exon_seq;
@@ -114,21 +130,26 @@ sub get_cdna_seq {
 ####
 sub create_gene {
     my ($gene_struct_mod, $cds_pointer_lend, $cds_pointer_rend) = @_;
+    
+    #use Data::Dumper;
+    #print STDERR Dumper($gene_struct_mod) . "CDS: $cds_pointer_lend, $cds_pointer_rend\n";
+        
     my $strand = $gene_struct_mod->{strand};
-    my @exons = @{$gene_struct_mod->{exons}};
+    my @exons = sort {$a->[0]<=>$b->[0]} @{$gene_struct_mod->{exons}};
     if ($strand eq '-') {
         @exons = reverse (@exons);
     }
     my $mRNA_pointer_lend = 1;
     my $mRNA_pointer_rend = 0;
     my $gene_obj = new Gene_obj();
     foreach my $coordset_ref (@exons) {
-        my ($coord1, $coord2) = @$coordset_ref;
+        my ($coord1, $coord2) = sort {$a<=>$b} @$coordset_ref;
         my ($end5, $end3) = ($strand eq '+') ? ($coord1, $coord2) : ($coord2, $coord1);
         my $exon_obj = new mRNA_exon_obj($end5, $end3);
         my $exon_length = ($coord2 - $coord1 + 1);
         $mRNA_pointer_rend = $mRNA_pointer_lend + $exon_length - 1;
         ## see if cds is within current cDNA range.
+        #print STDERR "mRNA coords: $mRNA_pointer_lend-$mRNA_pointer_rend\n";
         if ( $cds_pointer_rend >= $mRNA_pointer_lend && $cds_pointer_lend <= $mRNA_pointer_rend) { #overlap
             my $diff = $cds_pointer_lend - $mRNA_pointer_lend;
             my $delta_lend = ($diff >0) ? $diff : 0;

PerlLib/Fastq_reader.pm

@@ -3,8 +3,6 @@ package Fastq_reader;
 use strict;
 use warnings;
 
-use PerlIO::gzip;
-
 sub new {
     my ($packagename, $fastqFile) = @_;
 
@@ -24,7 +22,11 @@ sub new {
 	}
 	else {
 		if ( $fastqFile =~ /\.gz$/ ) {
-		    open ($filehandle, "<:gzip", $fastqFile) or die "Error: Couldn't open compressed $fastqFile\n";
+		    open ($filehandle, "gunzip -c $fastqFile | ") or die "Error: Couldn't open compressed $fastqFile\n";
+        }
+        elsif ($fastqFile =~ /\.bz2$/) {
+            open ($filehandle, "bunzip2 -c $fastqFile | ") or die "Error, couldn't open compressed $fastqFile $!";
+            
         } else {
             open ($filehandle, $fastqFile) or die "Error: Couldn't open $fastqFile\n";
         }

PerlLib/GFF3_alignment_utils.pm

@@ -1,14 +1,20 @@
+#!/usr/bin/env perl
+
 package GFF3_alignment_utils;
 
 use strict;
 use warnings;
 use Carp;
 
 use Gene_obj;
+use Gene_obj_indexer;
 use CDNA::Alignment_segment;
 use CDNA::CDNA_alignment;
+use File::Basename;
+
+__run_test() unless caller;
 
-sub index_GFF3_alignment_objs {
+sub index_alignment_objs {
     my ($gff3_alignment_file, $genome_alignment_indexer_href) = @_;
     
     unless ($gff3_alignment_file && -s $gff3_alignment_file) {
@@ -17,10 +23,11 @@ sub index_GFF3_alignment_objs {
     unless (ref $genome_alignment_indexer_href) {
         confess "Error, need genome indexer href as param ";
     }
-
     
-    	
+       	
 	my %genome_trans_to_alignment_segments;
+    my %trans_to_gene_id;
+
 	
 	open (my $fh, $gff3_alignment_file) or die "Error, cannot open file $gff3_alignment_file";
 	while (<$fh>) {
@@ -69,17 +76,19 @@ sub index_GFF3_alignment_objs {
 		
 		my ($end5, $end3) = ($orient eq '+') ? ($lend, $rend) : ($rend, $lend);
         
-        $info =~ /Target=\S+ (\d+) (\d+) ([\+\-])/ or die "Error, cannot extract match coordinates from info: $info";
+        $info =~ /Target=\S+ (\d+) (\d+)/ or die "Error, cannot extract match coordinates from info: $info";
         my $cdna_seg_lend = $1;
         my $cdna_seg_rend = $2;
-        my $cdna_orient = $3; # always set to + in pasa
+        
+        ($cdna_seg_lend, $cdna_seg_rend) = sort {$a<=>$b} ($cdna_seg_lend, $cdna_seg_rend); # always + orient for transcript coords.
+        
         
         my $alignment_segment = new CDNA::Alignment_segment($end5, $end3, $cdna_seg_lend, $cdna_seg_rend, $per_id);
         
 
-        push (@{$genome_trans_to_alignment_segments{$scaff}->{$gene_id}}, $alignment_segment);
+        push (@{$genome_trans_to_alignment_segments{$scaff}->{$trans_id}}, $alignment_segment);
         
-		
+		$trans_to_gene_id{$trans_id} = $gene_id;
         
 	}
     
@@ -109,14 +118,69 @@ sub index_GFF3_alignment_objs {
             $cdna_alignment_obj->set_acc($alignment_acc);
             $cdna_alignment_obj->{genome_acc} = $scaff;
             
-            $genome_alignment_indexer_href->{$alignment_acc} = $cdna_alignment_obj;
+            my $gene_id = $trans_to_gene_id{$alignment_acc} or confess "Error no gene_id for acc: $alignment_acc";
+
+            $cdna_alignment_obj->{gene_id} = $gene_id;
+            
+
+            $cdna_alignment_obj->{source} = basename($gff3_alignment_file);
             
+            if (ref $genome_alignment_indexer_href eq "Gene_obj_indexer") {
+                $genome_alignment_indexer_href->store_gene($alignment_acc, $cdna_alignment_obj);
+                
+            }
+            else {
+                                
+                $genome_alignment_indexer_href->{$alignment_acc} = $cdna_alignment_obj;
+            }
             push (@{$scaff_to_align_list{$scaff}}, $alignment_acc);
         }
     }
 
     return(%scaff_to_align_list);
 }
 
+
+
+#################
+## Testing
+#################
+
+
+sub __run_test {
+    
+    my $usage = "usage: $0 file.alignment.gff3\n\n";
+
+    my $gff3_file = $ARGV[0] or die $usage;
+
+    my $indexer = {};
+    my %scaff_to_alignments = &index_alignment_objs($gff3_file, $indexer);
+    
+    
+    foreach my $scaffold (keys %scaff_to_alignments) {
+
+        my @align_ids = @{$scaff_to_alignments{$scaffold}};
+
+        foreach my $align_id (@align_ids) {
+            my $cdna_obj = $indexer->{$align_id};
+
+            print $cdna_obj->toString();
+        }
+    }
+    
+    
+    
+    exit(0);
+    
+
+
+}
+
+
+
+
+
+
+
 1; #EOM
 

PerlLib/GFF3_utils.pm

@@ -66,7 +66,7 @@ sub index_GFF3_gene_objs {
 
         unless ($feat_type) { die "Error, $_, no feat_type: line\[$_\]"; }
         
-        unless ($feat_type =~ /^(gene|mRNA|CDS|exon)$/) { next;} 
+        unless ($feat_type =~ /^(gene|mRNA|transcript|CDS|exon)$/) { next;} 
 
         $gene_info = uri_unescape($gene_info);
         
@@ -107,7 +107,7 @@ sub index_GFF3_gene_objs {
                 
         # print "id: $id, parent: $parent\n";
         
-        if ($feat_type eq 'mRNA') {
+        if ($feat_type eq 'mRNA' || $feat_type eq 'transcript') {
             ## just get the identifier info
             $transcript_to_gene{$id} = $parent;
             next;