Unified Matlab and Python visualization

* Matlab will use the python visualization engine. This commit prepares
  future developments for spladder_viz. Functions loading data in a
  matlab or python specific manner have been moved into an identity
  module.

Author: Andre Kahles

.gitignore

@@ -2,3 +2,4 @@
 *.bam
 *.bai
 *.mexa64
+.history

matlab/modules/classes/__init__.py

@@ -0,0 +1,69 @@
+import scipy as sp
+
+class Event:
+
+    def __init__(self, event_type, chr=None, strand=None):
+        
+        self.event_type = event_type
+        self.chr = chr
+        self.strand = strand
+        self.strain = ''
+        self.exons1 = sp.zeros((0, 2), dtype = 'int')
+        self.exons2 = sp.zeros((0, 2), dtype = 'int')
+        self.exons1_col = sp.zeros((2, 0), dtype = 'int')
+        self.exons2_col = sp.zeros((2, 0), dtype = 'int')
+        self.p_values = None
+        self.gene_name = None
+        self.transcript_type = None
+        self.num_detected = None
+        self.id = None
+        self.detected = None
+        self.verified = []
+        self.num_verfied = 0
+        self.confirmed = []
+
+    def get_len(self, trafo=False):
+
+        if trafo:
+            return max(self.exons1_col.max(), self.exons2_col.max()) - min(self.exons1_col.min(), self.exons2_col.min())
+        else:
+            return max(self.exons1.max(), self.exons2.max()) - min(self.exons1.min(), self.exons2.min())
+
+    def get_inner_coords(self, trafo=False):
+        
+        if self.event_type == 'mult_exon_skip':
+            if trafo:
+                return sp.sort(sp.unique(sp.r_[sp.sort(self.exons2_col.ravel())[1:4], sp.sort(self.exons2_col.ravel())[-4:-1]]))
+                #return sp.unique(self.exons2_col.ravel())[1:-1]
+            else:
+                return sp.sort(sp.unique(sp.r_[sp.sort(self.exons2.ravel())[1:4], sp.sort(self.exons2.ravel())[-4:-1]]))
+                #return sp.unique(self.exons2.ravel())[1:-1]
+        elif self.event_type == 'mutex_exons':
+            if trafo:
+                return sp.sort(sp.r_[self.exons1_col.ravel()[1:4], self.exons2_col[1, :], self.exons1_col[2, 0]])
+            else:
+                return sp.sort(sp.r_[self.exons1.ravel()[1:4], self.exons2[1, :], self.exons1[2, 0]])
+        else:
+            if trafo:
+                return sp.sort(sp.unique(sp.r_[sp.sort(self.exons1_col.ravel())[1:-1], sp.sort(self.exons2_col.ravel())[1:-1]]))
+            else:
+                return sp.sort(sp.unique(sp.r_[sp.sort(self.exons1.ravel())[1:-1], sp.sort(self.exons2.ravel())[1:-1]]))
+            
+        
+
+    def get_coords(self, trafo=False):
+        
+        if self.event_type != 'mult_exon_skip':
+            if trafo:
+                #return sp.sort(sp.unique(sp.c_[self.exons1_col.ravel(), self.exons2_col.ravel()]))
+                return sp.sort(sp.r_[self.exons1_col.ravel(), self.exons2_col.ravel()])
+            else:
+                #return sp.sort(sp.unique(sp.c_[self.exons1.ravel(), self.exons2.ravel()]))
+                return sp.sort(sp.r_[self.exons1.ravel(), self.exons2.ravel()])
+        else:
+            if trafo:
+                return sp.sort(sp.r_[self.exons1_col.ravel()[:4], self.exons2_col.ravel()[-4:]])
+            else:
+                return sp.sort(sp.r_[self.exons1.ravel()[:4], self.exons2.ravel()[-4:]])
+            
+        

matlab/modules/classes/event.py

@@ -1,6 +1,8 @@
 import scipy as sp
 
-from .utils import *
+from ..utils import *
+from .segmentgraph import Segmentgraph
+from .splicegraph import Splicegraph
 
 class Gene:
 
@@ -13,9 +15,9 @@ def __init__(self, name=None, start=None, stop=None, chr=None, strand=None, sour
         self.strand = strand
         self.transcripts = []
         self.source = source
-        self.splicegraph = None
-        self.segmentgraph = None
-        self.gene_type=gene_type
+        self.splicegraph = Splicegraph()
+        self.segmentgraph = Segmentgraph()
+        self.gene_type = gene_type
         self.is_alt = None
         self.is_alt_spliced = None
 

matlab/modules/gene.py renamed to matlab/modules/classes/gene.py

@@ -0,0 +1,44 @@
+import scipy as sp
+
+class Segmentgraph:
+
+    def __init__(self, gene = None):
+    
+        self.segments = sp.zeros((2, 0), dtype='int')
+        self.seg_match = sp.zeros((0, 0), dtype='bool')
+        self.seg_edges = sp.zeros((0, 0), dtype='bool')
+
+        if gene is not None:
+            self.from_gene(gene)
+
+    def from_gene(self, gene): 
+
+        sg = gene.splicegraph.vertices
+        breakpoints = sp.unique(sg.ravel())
+        self.segments = sp.zeros((2, 0), dtype='int')
+        for j in range(1, breakpoints.shape[0]):
+            s = sp.sum(sg[0, :] < breakpoints[j])
+            e = sp.sum(sg[1, :] < breakpoints[j])
+            if s > e:
+                self.segments = sp.c_[self.segments, [breakpoints[j-1], breakpoints[j]]]
+
+        ### match nodes to segments
+        self.seg_match = sp.zeros((0, sg.shape[1]), dtype='bool')
+        for j in range(sg.shape[1]):
+            tmp = ((sg[0, j] <= self.segments[0, :]) & (sg[1, j] >= self.segments[1, :]))
+            if self.seg_match.shape[0] == 0:
+                self.seg_match = tmp.copy().reshape((1, tmp.shape[0]))
+            else:
+                self.seg_match = sp.r_[self.seg_match, tmp.reshape((1, tmp.shape[0]))]
+
+        ### create edge graph between segments
+        self.seg_edges = sp.zeros((self.segments.shape[1], self.segments.shape[1]), dtype='bool')
+        k, l = sp.where(sp.triu(gene.splicegraph.edges))
+        for m in range(k.shape[0]):
+            ### donor segment
+            d = sp.where(self.seg_match[k[m], :])[0][-1]
+            ### acceptor segment
+            a = sp.where(self.seg_match[l[m], :])[0][0]
+            self.seg_edges[d, a] = True
+
+

matlab/modules/classes/segmentgraph.py

@@ -0,0 +1,230 @@
+import scipy as sp
+
+from ..utils import *
+
+import pdb
+
+class Splicegraph:
+
+    def __init__(self, gene = None):
+        
+        self.vertices = sp.zeros((2, 0), dtype='int')
+        self.edges = sp.zeros((0, 0), dtype='int')
+        self.terminals = sp.zeros((2, 0), dtype='int')
+
+        if gene:
+            self.from_gene(gene)
+
+    def get_len(self):
+        
+        return self.vertices.shape[1]
+
+    def new_edge(self):
+
+        self.edges = sp.c_[self.edges, sp.zeros((self.edges.shape[0], 1), dtype='int')]
+        self.edges = sp.r_[self.edges, sp.zeros((1, self.edges.shape[1]), dtype='int')]
+    
+    def subset(self, keep_idx):
+        
+        self.vertices = self.vertices[:, keep_idx]
+        self.edges = self.edges[keep_idx, :][:, keep_idx]
+        self.terminals = self.terminals[:, keep_idx]
+
+    def update_terminals(self):
+        
+        self.terminals = sp.zeros(self.vertices.shape, dtype='int')
+        self.terminals[0, sp.where(sp.sum(sp.tril(self.edges), axis=1) == 0)[0]] = 1
+        self.terminals[1, sp.where(sp.sum(sp.triu(self.edges), axis=1) == 0)[0]] = 1
+
+    def reorder(self, idx):
+        
+        self.vertices = self.vertices[:, idx]
+        self.edges = self.edges[idx, :][:, idx]
+        self.terminals = self.terminals[:, idx]
+
+    def sort(self):
+        
+        s_idx = sp.lexsort([self.vertices[1, :], self.vertices[0, :]])
+        self.reorder(s_idx)
+
+    def from_gene(self, gene):
+        
+        for transcript_idx in range(len(gene.transcripts)):
+            exon_start_end = gene.exons[transcript_idx]
+            
+            ### only one exon in the transcript
+            if exon_start_end.shape[0] == 1:
+                exon1_start = exon_start_end[0, 0]
+                exon1_end = exon_start_end[0, 1]
+
+                if self.vertices.shape[1] == 0:
+                    self.vertices = sp.array([[exon1_start], [exon1_end]], dtype='int')
+                    self.edges = sp.array([[0]], dtype='int')
+                else:
+                    self.vertices = sp.c_[self.vertices, [exon1_start, exon1_end]]
+                    self.new_edge()
+            ### more than one exon in the transcript
+            else:
+                for exon_idx in xrange(exon_start_end.shape[0] - 1):
+                    exon1_start = exon_start_end[exon_idx , 0]
+                    exon1_end = exon_start_end[exon_idx, 1]
+                    exon2_start = exon_start_end[exon_idx + 1, 0]
+                    exon2_end = exon_start_end[exon_idx + 1, 1]
+          
+                    if self.vertices.shape[1] == 0:
+                        self.vertices = sp.array([[exon1_start, exon2_start], [exon1_end, exon2_end]], dtype='int')
+                        self.edges = sp.array([[0, 1], [1, 0]], dtype='int')
+                    else:
+                        exon1_idx = -1
+                        exon2_idx = -1
+                        ### check if current exon already occurred
+                        for idx in range(self.vertices.shape[1]):
+                            if ((self.vertices[0, idx] == exon1_start) and (self.vertices[1, idx] == exon1_end)):
+                                 exon1_idx = idx
+                            if ((self.vertices[0, idx] == exon2_start) and (self.vertices[1, idx] == exon2_end)):
+                                 exon2_idx = idx
+
+                        ### both exons already occured -> only add an edge
+                        if (exon1_idx != -1) and (exon2_idx != -1):
+                            self.edges[exon1_idx, exon2_idx] = 1
+                            self.edges[exon2_idx, exon1_idx] = 1
+                        else:
+                            ### 2nd exon occured
+                            if ((exon1_idx == -1) and (exon2_idx != -1)):
+                                self.vertices = sp.c_[self.vertices, [exon1_start, exon1_end]]
+                                self.new_edge()
+                                self.edges[exon2_idx, -1] = 1
+                                self.edges[-1, exon2_idx] = 1
+                            ### 1st exon occured
+                            elif ((exon2_idx == -1) and (exon1_idx != -1)):
+                                self.vertices = sp.c_[self.vertices, [exon2_start, exon2_end]]
+                                self.new_edge()
+                                self.edges[exon1_idx, -1] = 1
+                                self.edges[-1, exon1_idx] = 1
+                            ### no exon occured
+                            else:
+                                assert((exon1_idx == -1) and (exon2_idx == -1))
+                                self.vertices = sp.c_[self.vertices, [exon1_start, exon1_end]]
+                                self.vertices = sp.c_[self.vertices, [exon2_start, exon2_end]]
+                                self.new_edge()
+                                self.new_edge()
+                                self.edges[-2, -1] = 1
+                                self.edges[-1, -2] = 1
+
+        ### take care of the sorting by exon start
+        s_idx = sp.argsort(self.vertices[0, :])
+        self.vertices = self.vertices[:, s_idx]
+        self.edges = self.edges[s_idx, :][:, s_idx]
+        self.terminals = sp.zeros(self.vertices.shape, dtype='int')
+        self.terminals[0, sp.where(sp.tril(self.edges).sum(axis=1) == 0)[0]] = 1
+        self.terminals[1, sp.where(sp.triu(self.edges).sum(axis=1) == 0)[0]] = 1
+
+
+    def from_matfile(self, mat_struct):
+        """generates a splicing graph structure from a matfile structure"""
+
+        self.vertices = mat_struct['splicegraph'][0, 0].astype('int')
+        self.edges = mat_struct['splicegraph'][0, 1].astype('int')
+        self.terminals = mat_struct['splicegraph'][0, 2].astype('int')
+        
+
+    def add_intron(self, idx1, flag1, idx2, flag2):
+        """adds new introns into splicegraph between idx1 and idx2"""
+
+        ### if flag1, all end terminal exons in idx1 are preserved
+        ### if flag2, all start terminal exons in idx2 are preserved
+
+        if idx2.shape[0] > 0:
+            adj_mat = sp.triu(self.edges)
+
+            if flag1:
+                for i1 in idx1:
+
+                    ### if exon is end-terminal
+                    if sp.all(adj_mat[i1, :] == 0):
+
+                        self.vertices = sp.c_[self.vertices, self.vertices[:, i1]]
+
+                        self.new_edge()
+                        self.edges[:, -1] = self.edges[:, i1]
+                        self.edges[-1, :] = self.edges[i1, :]
+
+                        self.terminals = sp.c_[self.terminals, self.terminals[:, i1]]
+            if flag2:
+                for i2 in idx2:
+                    ### if exon is start-terminal
+                    if sp.all(adj_mat[:, i2] == 0):
+                        self.vertices = sp.c_[self.vertices, self.vertices[:, i2]]
+
+                        self.new_edge()
+                        self.edges[:, -1] = self.edges[:, i2]
+                        self.edges[-1, :] = self.edges[i2, :]
+
+                        self.terminals = sp.c_[self.terminals, self.terminals[:, i2]]
+
+        for i1 in idx1:
+            for i2 in idx2:
+                self.edges[i1, i2] = 1
+                self.edges[i2, i1] = 1
+        
+        self.uniquify()
+
+    def add_cassette_exon(self, new_exon, exons_pre, exons_aft):
+        ### exon_pre contains the indices of preceding exons
+        ### exon_aft contains the indices of successing exons
+        
+        self.vertices = sp.c_[self.vertices, new_exon]
+
+        self.new_edge()
+
+        self.edges[exons_pre, -1] = 1
+        self.edges[exons_aft, -1] = 1
+        self.edges[-1, :] = self.edges[:, -1].T
+
+        self.terminals = sp.c_[self.terminals, sp.zeros((2,), dtype='int')]
+
+
+    def add_intron_retention(self, idx1, idx2):
+        
+        adj_mat = sp.triu(self.edges)
+
+        self.vertices = sp.c_[self.vertices, sp.array([self.vertices[0, idx1], self.vertices[1, idx2]], dtype='int')]
+
+        self.new_edge()
+
+        adj_mat = sp.r_[adj_mat, sp.zeros((1, adj_mat.shape[1]), dtype='int')]
+        adj_mat = sp.c_[adj_mat, sp.zeros((adj_mat.shape[0], 1), dtype='int')]
+
+        ### check if adjacency matrix is symmetric
+        ### otherwise or is not justyfied
+        assert(sp.all(sp.all(adj_mat - (self.edges - adj_mat).T == 0)))
+
+        ### AK: under the assumption that our splice graph representation is symmetric
+        ### I preserve symmetry by using OR over the adj_mat column and row
+        
+        self.edges[:, -1] = adj_mat[:, idx1] | adj_mat[idx2, :].T
+        self.edges[-1, :] = adj_mat[:, idx1].T | adj_mat[idx2, :]
+
+        self.terminals = sp.c_[self.terminals, sp.array([self.terminals[0, idx1], self.terminals[1, idx2]], dtype='int')]
+
+    def uniquify(self):
+        # OUTPUT: splice graph that has been made unique on exons for each gene
+
+        self.sort()
+        (s_tmp, s_idx) = sort_rows(self.vertices.T, index=True)
+        self.vertices = s_tmp.T
+        self.edges = self.edges[s_idx, :][:, s_idx]
+        self.terminals = self.terminals[:, s_idx]
+
+        rm_idx = []
+        for j in range(1, self.vertices.shape[1]):
+            if sp.all(self.vertices[:, j-1] == self.vertices[:, j]):
+                self.edges[:, j] = self.edges[:, j-1] | self.edges[:, j]
+                self.edges[j, :] = self.edges[j-1, :] | self.edges[j, :]
+                rm_idx.append(j - 1)
+
+        keep_idx = sp.where(~sp.in1d(sp.array(range(self.vertices.shape[1])), rm_idx))[0]
+        self.vertices = self.vertices[:, keep_idx]
+        self.edges = self.edges[keep_idx, :][:, keep_idx]
+        self.terminals = self.terminals[:, keep_idx]
+