updates

scripts/lpls/lpls.py

@@ -403,7 +403,5 @@ def jk_lpls(X, Y, Z, a_max, nsets=None, alpha=.5):
         WWx[i,:,:] = W
         WWz[i,:,:] = L
         WWy[i,:,:] = Q
-        print "Q"
-        print Q
 
     return WWx, WWz, WWy

scripts/lpls/plots_lpls.py

@@ -31,8 +31,8 @@ def plot_corrloads(R, pc1=0,pc2=1,s=20, c='b', zorder=5,expvar=None,ax=None,draw
         pylab.xlabel(xstring)
         ystring = "Comp: %d   expl.var.:  %.1f " %(pc2+1, expvar[pc2])
         pylab.ylabel(ystring)
-    if labels:
+    if labels!=None:
         assert(len(labels)==R.shape[0])
         for name, r in zip(labels, R):
-            ax.text(r[pc1], r[pc2], "  " + name)
+            pylab.text(r[pc1], r[pc2], "  " + name)
     #pylab.show()

scripts/lpls/rpy_go.py

@@ -45,7 +45,7 @@ def goterms_from_gene(genelist, ontology=['BP'], garbage = ['IEA', 'ISS', 'ND'])
     _ONTOLOGIES = ['BP', 'CC', 'MF']
     assert(scipy.all([(code in _CODES) for code in garbage]))
     assert(scipy.all([(ont in _ONTOLOGIES) for ont in ontology]))
-    
+    have_these = rpy.r('as.list(GOTERM)').keys()
     goterms = {}
     for gene in genelist:
         goterms[gene] = []
@@ -53,8 +53,12 @@ def goterms_from_gene(genelist, ontology=['BP'], garbage = ['IEA', 'ISS', 'ND'])
         #print info
         if info:
             for term, desc in info.items():
+                if term not in have_these:
+                    print "GO miss:"
+                    print term
                 if desc['Ontology'] in ontology and desc['Evidence'] not in garbage:
                         goterms[gene].append(term)
+
     return goterms
 
 def genego_matrix(goterms, tmat, gene_ids, term_ids, func=min):
@@ -97,7 +101,12 @@ def goterm2desc(gotermlist):
     return term2desc
 
 def parents_dag(go_terms, ontology=['BP']):
-    """ Returns a list of lists representation of a GO DAG parents of goterms."""
+    """ Returns a list of lists representation of a GO DAG parents of goterms.
+
+    make the networkx graph by:
+        G = networkx.Digraph()
+        G = networkx.from_dict_of_lists(edge_dict, G)
+    """
     try:
         rpy.r.library("GOstats")
     except:
@@ -105,6 +114,11 @@ def parents_dag(go_terms, ontology=['BP']):
     assert(go_terms[0][:3]=='GO:')
 
     # go valid namespace
-    go_env = {'BP':rpy.r.BPPARENTS, 'MF':rpy.r.MFPARENTS, 'CC': rpy.r.CCPARENTS}
+    go_env = {'BP':rpy.r.GOBPPARENTS, 'MF':rpy.r.GOMFPARENTS, 'CC': rpy.r.GOCCPARENTS}
-    
+    graph = rpy.r.GOGraph(go_terms, go_env[ontology[0]])
-    
+    edges = rpy.r.edges(graph)
+    edges.pop('all')
+    edge_dict = {}
+    for head, nei in edges.items():
+        edge_dict[head] = nei.values()
+    return edge_dict

scripts/lpls/run_smoker.py

@@ -32,7 +32,7 @@ print "SAM done"
 qq = rpy.r('qobj<-qvalue(sam.out@p.value)')
 qvals = asarray(qq['qvalues'])
 # cut off
-co = 0.001
+co = 0.1
 index = where(qvals<0.01)[0]
 
 # Subset data
@@ -51,7 +51,7 @@ terms = list(terms)
 rpy.r.library("GOSim")
 # Go-term similarity matrix
 methods = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
-meth = methods[0]
+meth = methods[2]
 print "Term-term similarity matrix (method = %s)" %meth
 if meth=="CoutoEnriched":
     rpy.r('setEnrichmentFactors(alpha=0.1,beta=0.5)')
@@ -100,7 +100,6 @@ T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(Xr,Y,Z, a_max, alpha)
 dx,Rx,ssx= correlation_loadings(Xr, T, P)
 dx,Ry,ssx= correlation_loadings(Y, T, Q)
 cadx,Rz,ssx= correlation_loadings(Z.T, K, L)
-
 # Prediction error
 rmsep , yhat, class_error = cv_lpls(Xr, Y, Z, a_max, alpha=alpha)
 
@@ -112,11 +111,19 @@ tsqz = cx_stats.hotelling(Wz,L[:,:aopt])
 
 ## plots ##
 figure(1) #rmsep
-#bar()
+bar_w = .2
+bar_col = 'rgb'*5
+m = Y.shape[1]
+for a in range(m):
+    bar(arange(a_max)+a*bar_w+.1, rmsep[:,a], width=bar_w, color=bar_col[a])
+ylim([rmsep.min()-.05, rmsep.max()+.05])
+
 figure(2) # Hypoid correlations
+title('RMSEP')
 plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz/10.0, c='b', zorder=5, expvar=evz, ax=None)
 ax = gca()
-plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', zorder=5, expvar=evy, ax=ax)
+ylabels = DY.get_identifiers('_cat', sorted=True)
+plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', zorder=5, expvar=evy, ax=ax,labels=ylabels)
 
 figure(3)
 subplot(221)