Whatnot

scripts/lpls/plots_lpls.py

@@ -11,24 +11,30 @@ def plot_corrloads(R, pc1=0,pc2=1,s=20, c='b', zorder=5,expvar=None,ax=None,draw
     if ax==None or drawback==True:
         radius = 1
         center = (0,0)
-        c100 = matplotlib.patches.Circle(center,radius=radius,
+        c100 = matplotlib.patches.Circle(center,
-                                         facecolor='gray',
+                                         radius=radius,
-                                         alpha=.1,
+                                         facecolor=(0.97, .97, .97),
-                                         zorder=1)
+                                         zorder=1,
-        c50 = matplotlib.patches.Circle(center, radius=radius/2.0,
+                                         linewidth=1,
-                                        facecolor='gray',
+                                         edgecolor=(0,0,0))
-                                        alpha=.1,
+        
-                                        zorder=2)
+        c50 = matplotlib.patches.Circle(center,
+                                        radius=radius/2.0,
+                                        facecolor=(.85,.85,.85),
+                                        zorder=1,
+                                        linewidth=1,
+                                        edgecolor=(0,0,0))
+        
         ax = pylab.gca()
         ax.add_patch(c100)
         ax.add_patch(c50)
-        ax.axhline(lw=1.5,color='k')
+        ax.axhline(lw=1.5,color='k', zorder=4)
-        ax.axvline(lw=1.5,color='k')
+        ax.axvline(lw=1.5,color='k', zorder=4)
 
     # corrloads
     ax.scatter(R[:,pc1], R[:,pc2], s=s, c=c,zorder=zorder, **kwds)
-    ax.set_xlim([-1,1])
+    ax.set_xlim([-1.1,1.1])
-    ax.set_ylim([-1,1])
+    ax.set_ylim([-1.1,1.1])
     if expvar!=None:
         xstring = "Comp: %d   expl.var:  %.1f " %(pc1+1, expvar[pc1])
         pylab.xlabel(xstring)
@@ -51,7 +57,7 @@ def dag(terms, ontology):
     #setattr(dag, "_ontology", ontology)
     return dag
 
-def plot_dag(dag, node_color='b', node_size=30,with_labels=False,nodelist=None,pos=None):
+def plot_dag(dag, node_color='b', node_size=30,with_labels=False,nodelist=None,pos=None,**kwd):
 
     rpy.r.library("GOstats")
     
@@ -76,7 +82,7 @@ def plot_dag(dag, node_color='b', node_size=30,with_labels=False,nodelist=None,p
     if len(node_color)>1:
         assert(len(node_color)==len(nodelist))
         
-    nx.draw_networkx(G,pos, with_labels=with_labels, node_size=node_size, node_color=node_color, nodelist=nodelist)
+    nx.draw_networkx(G,pos, with_labels=with_labels, node_size=node_size, node_color=node_color, nodelist=nodelist, **kwd)
     return pos
 
 

scripts/lpls/rpy_go.py

@@ -4,7 +4,7 @@ import rpy
 silent_eval = rpy.with_mode(rpy.NO_CONVERSION, rpy.r)
 import collections
 
-def goterms_from_gene(genelist, ontology='BP', garbage=None, ic_cutoff=2.0):
+def goterms_from_gene(genelist, ontology='BP', garbage=['IEA'], ic_cutoff=2.0, verbose=False):
     """ Returns the go-terms from a specified genelist (Entrez id).
 
     Recalculates the information content if needed based on selected evidence codes.
@@ -30,10 +30,17 @@ def goterms_from_gene(genelist, ontology='BP', garbage=None, ic_cutoff=2.0):
     ddef = False
     if ontology=='BP' and garbage!=None:
         # This is for ont=BP and garbage =['IEA', 'ISS', 'ND']
-        rpy.r.load("ICsBPIMP_IGI_IPI_ISS_IDA_IEP_TAS_NAS_IC.rda")
+        rpy.r.load("ICsBP_small.rda") # Excludes IEA
         ic = rpy.r.assign("IC",rpy.r.IC, envir=rpy.r.GOSimEnv)
-        print len(ic)
+        max_val = 0
+        for key, val in ic.items():
+            if val != scipy.inf:
+                if val>max_val:
+                    max_val = val
+        for key, val in ic.items():
+            ic[key] = val/max_val
     else:
+        # NB! this IC is just for BP
         ic = rpy.r('get("IC", envir=GOSimEnv)')
     print "loading GO definitions environment"
     
@@ -42,25 +49,57 @@ def goterms_from_gene(genelist, ontology='BP', garbage=None, ic_cutoff=2.0):
     dd = 0
     ii = 0
     jj = 0
-    all = rpy.r.mget(gene_ids, rpy.r.GOENTREZID2GO,ifnotfound="NA")
+    kk = 0
+    all = rpy.r.mget(genelist, rpy.r.GOENTREZID2GO,ifnotfound="NA")
+    n_ic = len(ic)
+    print "Number of terms with IC: %d" %n_ic
+    stopp = False
     for gene, terms in all.items():
+        if verbose:
+            print "\n\n ======ITEM========\n"
+            print "Gene: " + str(gene)
+            print "Number of terms: %d " %len(terms)
+            print terms
+            print "---\n"
+            if stopp:
+                1/0
         if terms!="NA":
-            for term,desc in terms.items():
+            for term, desc in terms.items():
-                if desc['Ontology'].lower() == ontology and term in ic:
+                if verbose:
-                    if ic[term]>.88:
+                    print "\nChecking term: " + str(term)
+                    print "With description: " + str(desc)
+                if desc['Ontology'].lower() == ontology.lower() and term in ic:
+                    if ic[term]>ic_cutoff:
+                        #print ic[term]
                         jj+=1
+                        if verbose:
+                            print "too high" + str((gene, term))
+                            stopp = True
                         continue
-                    cc+=1
+                    cc += 1
+                    if verbose:
+                        print "accepted" + str((gene, term))
                     gene2terms[gene].append(term)
                 else:
+                    if verbose:
+                        print "Not accepted: " + str((gene, term))
+                    if term not in ic:
+                        if verbose:
+                            print "Not in IC: "  + str((gene, term))
+                        kk+=1
+                    if desc['Ontology'].lower() != ontology:
+                        if verbose:
+                            print "Not in Ontology" + str((gene, term))
                         dd+=1
         else:
             ii+=1
                     
-        
+    print "Number of genes total: %d" %len(all)
-    print "\nNumber of genes without annotation: %d" %ii 
+    print "\nNumber of genes without annotation: (%d (NA))" %ii
-    print "\nNumber of genes not in %s : %d " %(ontology, dd)
+    print "\nNumber of terms with annoation but no IC: %d" %kk
-    print "\nNumber of genes with too high IC : %d " %jj
+    print "\nNumber of terms not in %s : %d " %(ontology, dd)
+    print "\nNumber of terms with too high IC : %d " %jj
+    print "\n Number of accepted terms: %d" %cc
     
     return gene2terms
 
@@ -156,7 +195,7 @@ def parents_dag(go_terms, ontology=['BP']):
 def gene_GO_hypergeo_test(genelist,universe="entrezUniverse",ontology="BP",chip = "hgu133a",pval_cutoff=0.01,cond=False,test_direction="over"):
     
     #assert(scipy.alltrue([True for i in genelist if i in universe]))
-    universeGeneIds=universe
+    universeGeneIds = universe
     params = rpy.r.new("GOHyperGParams",
                        geneIds=genelist,
                        annotation="hgu133a",
@@ -222,5 +261,3 @@ def R_PLS(x,y,ncomp=3, validation='"LOO"'):
     result = rpy.r(callstr)
     return result
 
-
-def gogene()

scripts/lpls/run_smoker.py

@@ -1,4 +1,4 @@
-import sys
+import sys,time,cPickle
 import rpy
 from pylab import gca, figure, subplot,plot
 from scipy import *
@@ -9,70 +9,214 @@ import rpy_go
 sys.path.append("../../fluents") # home of dataset
 sys.path.append("../../fluents/lib") # home of cx_stats
 sys.path.append("/home/flatberg/fluents/scripts/lpls")
+sys.path.append("/home/flatberg/pyblm/")
 import dataset
 import cx_stats
-from engines import nipals_lpls
+import pyblm
-from validation import lpls_val, lpls_jk
+from pyblm.engines import nipals_lpls, pls
+from pyblm.crossvalidation import lpls_val, lpls_jk
+from pyblm.statistics import pls_qvals
 from plots_lpls import plot_corrloads, plot_dag
 import plots_lpls
 
+
+def iqr(X, axis=0):
+    """Interquartile range filtering."""
+    def _iqr(c):
+        return stats.scoreatpercentile(c, 75) - stats.scoreatpercentile(c, 25)
+    return apply_along_axis(_iqr, axis, X)
+
+
+
 # Possible outliers
 # http://www.pubmedcentral.nih.gov/articlerender.fcgi?tool=pubmed&pubmedid=16817967
 sample_outliers = ['OV:NCI_ADR_RES', 'CNS:SF_295', 'CNS:SF_539', 'RE:SN12C', 'LC:NCI_H226', 'LC:NCI_H522', 'PR:PC_3', 'PR:DU_145']
 
+outlier = 'ME:LOXIMVI' # 19
 
 ####### OPTIONS ###########
 # data
+
 chip = "hgu133a"
 use_data = 'uma'
-subset = 'plsr'
+#use_data = 'scherf'
-small_test = False
+#use_data = 'uma'
-use_sbg_subset = True # the sandberg nci-Ygroups subset
-std_y = True
-std_z = False
-# go
-ontology = "bp"
-min_genes = 5
-similarities = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
-meth = similarities[2]
-go_term_sim = "OA"
-# lpls
-a_max = 5
-aopt = 2
-xz_alpha = .4
-w_alpha = .3
-mean_ctr = [2, 0, 2]
-nsets = None
-qval_cutoff = 0.01
-n_iter = 200
 
-alpha_check = False
+if use_data == 'scherf':
-calc_rmsep = False
+    data_cached = False
+    use_saved_plsr_result = False
+    subset = 'plsr'
+    small_test = False
+    use_sbg_subset = True # the sandberg nci-Ygroups subset
+    std_y = False
+    std_z = False
+    # go
+    ontology = "bp"
+    min_genes = 5
+    similarities = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
+    meth = similarities[2]
+    go_term_sim = "OA"
+    # lpls
+    a_max = 10
+    aopt = 4
+    aopt = 2 # doubling-time
+    xz_alpha = .5
+    w_alpha = .3
+    center_axis = [2, 0, 2]
+    zorth = True
+    nsets = None
+    qval_cutoff = 0.1
+    n_iter = 50
     
+    alpha_check = True
+    calc_rmsep = True
+    
+    bevel_check = False
+    
+    save_calc = True
+
+elif use_data == 'uma':
+    data_cached = False
+    use_saved_plsr_result = False
+    subset = 'iqr'
+    small_test = False
+    use_sbg_subset = True # the sandberg nci-Ygroups subset
+    std_y = False
+    std_z = False
+    # go
+    ontology = "bp"
+    min_genes = 5
+    similarities = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
+    meth = similarities[2]
+    go_term_sim = "OA"
+    # lpls
+    a_max = 10
+    aopt = 5
+    xz_alpha = .5
+    w_alpha = .3
+    center_axis = [2, 0, 2]
+    zorth = True
+    nsets = None
+    qval_cutoff = 0.01
+    n_iter = 50
+    
+    alpha_check = True
+    calc_rmsep = True
+    
+    bevel_check = False
+    
+    save_calc = True
+
+elif use_data == 'smoker':
+    data_cached = False
+    use_saved_plsr_result = False
+    #subset = 'plsr'
+    subset = 'plsr'
+    small_test = False
+    use_sbg_subset = False # the sandberg nci-Ygroups subset
+    std_y = False
+    std_z = False
+    # go
+    ontology = "bp"
+    min_genes = 5
+    similarities = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
+    meth = similarities[2]
+    go_term_sim = "OA"
+    # lpls
+    a_max = 5
+    aopt = 2
+    xz_alpha = .5
+    w_alpha = .3
+    center_axis = [2, 0, 2]
+    zorth = True
+    nsets = None
+    qval_cutoff = 0.01
+    n_iter = 50
+    
+    alpha_check = True
+    calc_rmsep = True
+    
+    bevel_check = False
+    
+    save_calc = True
+else:
+    raise ValueError
+
+print "Using options for : " + use_data
 
 ######## DATA ##########
 if use_data=='smoker':
     # full smoker data
-    DX = dataset.read_ftsv(open("../../data/smokers-full/Smokers.ftsv"))
+    DX = dataset.read_ftsv(open("/home/flatberg/datasets/smokers/full/Smokers.ftsv"))
-    DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv"))
+    DY = dataset.read_ftsv(open("/home/flatberg/datasets/smokers/full/Yg.ftsv"))
-    Y = DY.asarray().astype('d')
+    DYr = dataset.read_ftsv(open("/home/flatberg/datasets/smokers/full/Ypy.ftsv"))
+    Y = DYr.asarray().astype('d')
     gene_ids = DX.get_identifiers('gene_ids', sorted=True)
+    sample_ids = DY.get_identifiers('_patient', sorted=True)
+
 elif use_data=='scherf':
-    DX = dataset.read_ftsv(open("../../data/scherf/scherfX.ftsv"))
+    print "hepp"
-    DY = dataset.read_ftsv(open("../../data/scherf/scherfY.ftsv"))
+    #DX = dataset.read_ftsv(open("../../data/scherf/old_data/scherfX.ftsv"))
-    Yg = DY.asarray().astype('d')
+    #DY = dataset.read_ftsv(open("../../data/scherf/old_data/scherfY.ftsv"))
+    DX = dataset.read_ftsv(open("nci60/X5964.ftsv", "r"))
+    DYg = dataset.read_ftsv(open("../../data/uma/Yg133.ftsv"))
+    DYr = dataset.read_ftsv(open("../../data/uma/Yd.ftsv"))
+    Y = DYg.asarray().astype('d')
+    DY = DYg.copy()
+    Yg = Y
+    Yr = DYr.asarray().astype('d')
+    X = DX.asarray()
     gene_ids = DX.get_identifiers('gene_ids', sorted=True)
+    sample_ids = DY.get_identifiers('cline', sorted=True)
+
 elif use_data=='staunton':
     pass
-elif use_data=='uma':
-    DX = dataset.read_ftsv(open("../../data/uma/X133.ftsv"))
-    DYg = dataset.read_ftsv(open("../../data/uma/Yg133.ftsv"))
-    DY = dataset.read_ftsv(open("../../data/uma/Yd.ftsv"))
-    Y = DY.asarray().astype('d')
-    Yg = DYg.asarray().astype('d')
-    gene_ids = DX.get_identifiers('gene_ids', sorted=True)
 
-# use subset with defined GO-terms
+elif use_data=='uma':
+    DX = dataset.read_ftsv(open("/home/flatberg/datasets/uma/X133.ftsv"))
+    DYg = dataset.read_ftsv(open("/home/flatberg/datasets/uma/Yg133.ftsv"))
+    DYr = dataset.read_ftsv(open("/home/flatberg/datasets/uma/Yd.ftsv"))
+    X = DX.asarray()
+    Y = DYg.asarray().astype('d')
+    DY = DYg.copy()
+    Yg = Y
+    Yr = DYr.asarray().astype('d')
+    gene_ids = DX.get_identifiers('gene_ids', sorted=True)
+    sample_ids = DY.get_identifiers('cline', sorted=True)
+
+else:
+    print "use_data argument: (%s) not valid" %use_method
+
+if use_sbg_subset and use_data in ['uma', 'scherf', 'staunton']:
+    print "Using sbg subset of cancers"
+    Y = Yg
+    Y_old = Y.copy()
+    Yr_old = Yr.copy()
+    X_old = X.copy() 
+    keep_samples = ['CN', 'ME', 'LE', 'CO', 'RE']
+    #keep_samples = ['CN', 'ME', 'LE', 'CO', 'RE']
+    sample_ids_original = DY.get_identifiers('cline', sorted=True)
+    sample_ids= [i for i in sample_ids if i[:2] in keep_samples]
+    rows_ind = [i for i,name in enumerate(sample_ids_original) if name[:2] in keep_samples]
+    # take out rows in X,Y
+    X = X[rows_ind,:]
+    Y = Y[rows_ind,:]
+    Yr = Yr[rows_ind,:]
+
+    #  identify redundant columns in Y
+    cols_ind = where(Y.sum(0)>1)[0]
+    Y = Y[:, cols_ind]
+
+    # create new datasets with updated idents
+    cat_ids = [name for i,name in enumerate(DYg.get_identifiers('_cancer', sorted=True)) if i in cols_ind]
+    DX = dataset.Dataset(X, [['cline', sample_ids], ['gene_ids', gene_ids]], name='Dxr')
+    DYg = dataset.CategoryDataset(Y, [['cline', sample_ids], ['_cancer', cat_ids]], name='Dyr')
+    DYr = dataset.Dataset(Yr, [['cline', sample_ids], ['_time', ['doubling_time']]], name='Dyrr')
+    DY_old = DY.copy()
+    DY = DYg
+    print "Now there are %d samples in X" %X.shape[0]
+
+# use subset of genes with defined GO-terms
 ic_all = 2026006.0 # sum of all ic in BP
 max_ic = -log(1/ic_all)
 ic_cutoff = -log(min_genes/ic_all)/max_ic
@@ -99,27 +243,41 @@ print "\n\nFiltering genes by Go terms "
 
 # use subset based on SAM,PLSR or (IQR)
 
-if subset=='sam':
+if subset=='plsr':
-    # select subset genes by SAM
+    print "plsr filter on genes"
-    rpy.r.library("siggenes")
+    if use_saved_plsr_result:
-    rpy.r.library("qvalue")
+        index = cPickle.load(open('plsr_index.pkl'))
-    rpy.r.assign("data", X.T)
+         # Subset data
-    cl = dot(DYg.asarray(), diag(arange(Yg.shape[1])+1)).sum(1)
+        X = X[:,index]
-    rpy.r.assign("cl", cl)
+
-    rpy.r.assign("B", 20)
+        gene_ids = [gid for i, gid in enumerate(gene_ids) if i in index]
-    # Perform a SAM analysis.
+        print "\nNumber of genes: %s" %len(gene_ids)
-    print "Starting SAM"
+        print "\nWorking on subset with %s genes " %len(gene_ids)
-    sam = rpy.r('sam.out<-sam(data=data,cl=cl,B=B,rand=123)')
+
-    print "SAM done"
+        # update valid go-terms
-    # Compute the q-values of the genes.
+        gene2goterms = rpy_go.goterms_from_gene(gene_ids, ic_cutoff=ic_cutoff)
-    qq = rpy.r('qobj<-qvalue(sam.out@p.value)')
+        all_terms = set()
-    qvals = asarray(qq['qvalues'])
+        for t in gene2goterms.values():
+            all_terms.update(t)
+        terms = list(all_terms)
+        print "\nNumber of go-terms: %s" %len(terms)
+        # update genelist
+        gene_ids = gene2goterms.keys()
+    else:
+        
+        print "Initial plsr qvals"
+        
+        xcal_tsq_x, xpert_tsq_x = pyblm.pls_qvals(X, Y, aopt=aopt, n_iter=n_iter, center_axis=[0,0], nsets=None)
+        qvals = pyblm.statistics._fdr(xcal_tsq_x, xpert_tsq_x, median)
+        
         # cut off
-    cutoff = 0.001
+        #sort_index = qvals.argsort()
-    index = where(qvals<cutoff)[0]
+        #index = sort_index[:800]
+        #qval_cutoff = qvals[sort_index[500]]
+        print "Using cuf off: %.2f" %qval_cutoff
+        index = where(qvals<qval_cutoff)[0]
         if small_test:
             index = index[:20]
-    
         # Subset data
         X = X[:,index]
 
@@ -137,10 +295,27 @@ if subset=='sam':
         # update genelist
         gene_ids = gene2goterms.keys()
     print "\nNumber of genes: %s" %len(gene_ids)
-elif subset=='plsr':
+
-    cx_stats.pls_qvals(X, Y)
+elif subset == 'iqr':
+    iqr_vals = iqr(X)
+    index = where(iqr_vals>1)[0]
+    X = X[:,index]
+    gene_ids = [gid for i, gid in enumerate(gene_ids) if i in index]
+    print "\nNumber of genes: %s" %len(gene_ids)
+    print "\nWorking on subset with %s genes " %len(gene_ids)
+    
+    # update valid go-terms
+    gene2goterms = rpy_go.goterms_from_gene(gene_ids, ic_cutoff=ic_cutoff)
+    all_terms = set()
+    for t in gene2goterms.values():
+        all_terms.update(t)
+    terms = list(all_terms)
+    print "\nNumber of go-terms: %s" %len(terms)
+    # update genelist
+    gene_ids = gene2goterms.keys()
 else:
     # noimp (smoker data is prefiltered)
+    print "No prefiltering on data used"
     pass 
 
 
@@ -154,40 +329,91 @@ if meth=="CoutoEnriched":
     aa = 0
     ba = 0
     rpy.r.setEnrichmentFactors(alpha = aa, beta =ba)
-rpytmat = rpy.with_mode(rpy.NO_CONVERSION, rpy.r.getTermSim)(terms, method=meth,verbose=False)
+if not data_cached:
-tmat = rpy.r.assign("haha", rpytmat)
+    rpytmat = rpy.with_mode(rpy.NO_CONVERSION, rpy.r.getTermSim)(terms, method=meth,verbose=False)
-print "\n Calculating Z matrix"
+    tmat = rpy.r.assign("haha", rpytmat)
-Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim=go_term_sim,term_sim=meth)
+    print "\n Calculating Z matrix"
-
+    Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim=go_term_sim,term_sim=meth)
-# update data (X) matrix
-newind = DX.get_indices('gene_ids', gene_ids)
-Xr = DX.asarray()[:,newind]
-
-
-######## LPLSR ######## 
-print "LPLSR ..."
-Y = Yg
-
-if use_sbg_subset:
-    Y_old = Y.copy()
-    Xr_old = Xr.copy() 
-    keep_samples = ['CN', 'ME', 'LE', 'CO', 'RE']
-    sample_ids = DY.get_identifiers('cline', sorted=True)
-    keep_ind = [i for i,name in enumerate(sample_ids) if name[:2] in keep_samples]
-    Xr = Xr[keep_ind,:]
-    Y = Y[keep_ind,:]
-    Y = Y[:, where(Y.sum(0)>1)[0]]
 
+    DZ = dataset.Dataset(Z, [['go-terms', terms], ['gene_ids', gene_ids]], name='Dz_'+str(meth))
+    # update data (X) matrix
+    newind = DX.get_indices('gene_ids', gene_ids)
+    Xr = DX.asarray()[:,newind]
+    DXr = dataset.Dataset(Xr, [['cline', sample_ids], ['gene_ids', gene_ids]], name='Dxr')
+else:
+    #DXr = dataset.read_ftsv(open('Xr.ftsv', 'r'))
+    newind = DX.get_indices('gene_ids', gene_ids)
+    Xr = DX.asarray()[:,newind]
+    DXr = dataset.Dataset(Xr, [['cline', sample_ids], ['gene_ids', gene_ids]], name='Dxr')
+    DY = dataset.read_ftsv(open('Y.ftsv', 'r'))
+    DZ = dataset.read_ftsv(open('Z.ftsv', 'r')) 
+    Xr = DXr.asarray()
+    Y = DY.asarray()
+    Z = DZ.asarray()
+    sample_ids = DX.get_identifiers('cline', sorted=True)
 
 # standardize Z?
 sdtz = False
 if sdtz:
-    Z = Z/Z.std(0)
+    DZ._array = DZ._array/Dz._array.std(0)
-
+sdty = False
-sdty = True
 if sdty:
-    Y = Y/Y.std(0)
+    DY._array = DY._array/DY._array.std(0)
-lpls_result = nipals_lpls(Xr,Y,Z, a_max,alpha=xz_alpha,mean_ctr=mean_ctr)
+
+
+# ##### PLS ONLY, CHECK FOR SIMILARITY BETWEEN W and Z #######
+if bevel_check:
+    Xr = DXr.asarray()
+    Y = DY.asarray()
+    from pylab import figure, scatter, xlabel, subplot,xticks,yticks
+    Xrcc = Xr - Xr.mean(0) - Xr.mean(1)[:,newaxis] + Xr.mean()
+    Zcc = Z - Z.mean(0) - Z.mean(1)[:,newaxis] + Z.mean()
+    Yc = Y - Y.mean(0)
+    xy_pls_result = pls(Xrcc, Yc, a_max)
+    xz_pls_result = pls(Xrcc.T, Zcc.T, a_max)
+    # check for linearity between scores of xz-result and W of xy-result
+    Wxy = xy_pls_result['W']
+    Txz = xz_pls_result['T']
+    figure()
+    n = 0
+    for i in range(a_max):
+        w = Wxy[:,i]
+        for j in range(a_max):
+            n += 1
+            t = Txz[:,j]
+            r2 = stats.corrcoef(w, t)[0,-1]
+            subplot(a_max, a_max, n)
+            scatter(w, t)
+            xticks([])
+            yticks([])
+            xlabel('(Wxy(%d), Tzx(%d)), r2: %.1f ' %(i+1,j+1,r2))
+# ####### LPLSR ######## 
+
+if save_calc and not data_cached:
+    print "Saving calculations"
+    import cPickle
+    fh =  open("g2go_s.pkl", "w")
+    cPickle.dump(gene2goterms, fh)
+    fh.close()
+    dataset.write_ftsv(open('Xs.ftsv', 'w'), DXr, decimals=7)
+    dataset.write_ftsv(open('Ysg.ftsv', 'w'), DY, decimals=7)
+    dataset.write_ftsv(open('Yspy.ftsv', 'w'), DYr, decimals=7)
+    dataset.write_ftsv(open('Zs.ftsv', 'w'), DZ, decimals=7)
+
+def read_calc():
+    import cPickle
+    fh =  open("g2go_s.pkl")
+    gene2goterms = cPickle.load(fh)
+    fh.close()
+    DXr = dataset.read_ftsv('Xu.ftsv')
+    DY = dataset.read_ftsv('Yu.ftsv')
+    DYr = dataset.read_ftsv('Ydu.ftsv')
+    DZ = dataset.read_ftsv('Zu.ftsv')
+    return DXr, DY, DYr, DZ, gene2goterms
+
+    
+print "LPLSR ..."
+lpls_result = nipals_lpls(Xr,Y,Z, a_max,alpha=xz_alpha, center_axis=center_axis, zorth=zorth)
 globals().update(lpls_result)
 
 # Correlation loadings
@@ -197,46 +423,75 @@ cadz,Rz,rssz = correlation_loadings(Z.T, W, L)
 
 # Prediction error
 if calc_rmsep:
-    rmsep , yhat, class_error = lpls_val(Xr, Y, Z, a_max, alpha=xz_alpha,mean_ctr=mean_ctr)
+    rmsep , yhat, class_error = pyblm.crossvalidation.lpls_val(Xr, Y, Z, a_max, alpha=xz_alpha,center_axis=center_axis, nsets=nsets,zorth=zorth)
-
 
+Alpha = arange(0.0, 1.01, .05)
 if alpha_check:
-    Alpha = arange(0.01, 1, .1)
     Rmsep,Yhat, CE = [],[],[]
     for a in Alpha:
         print "alpha %f" %a
-        rmsep , yhat, ce = lpls_val(Xr, Y, Z, a_max, alpha=a,mean_ctr=mean_ctr,nsets=nsets)
+        rmsep_a , yhat, ce = pyblm.lpls_val(Xr, Y, Z, a_max, alpha=a,
-        Rmsep.append(rmsep.copy())
+                                          center_axis=center_axis,nsets=nsets,
-        #Yhat.append(yhat.copy())
+                                          zorth=zorth)
-        #CE.append(ce.copy())
+        Rmsep.append(rmsep_a.copy())
+        Yhat.append(yhat.copy())
+        CE.append(ce.copy())
     Rmsep = asarray(Rmsep)
-    #Yhat = asarray(Yhat)
+    Yhat = asarray(Yhat)
     #CE = asarray(CE)
 
+random_alpha_check = True
+if random_alpha_check:
+    n_zrand = 100
+    RMS,YHAT, CEE = [],[],[]
+    zindex = arange(Z.shape[1])
+    for ii in range(n_zrand):
+        zind_rand = zindex.copy()
+        random.shuffle(zind_rand)
+        Zrand = Z[:,zind_rand]
+        #Alpha = arange(0.0, 1.1, .25)
+        Rmsep_r,Yhat_r, CE_r = [],[],[]
+        for a in Alpha:
+            print "Iter: %d  alpha %.2f" %(ii, a)
+            rmsep , yhat, ce = pyblm.lpls_val(Xr, Y, Zrand, a_max, alpha=a,center_axis=center_axis,nsets=nsets, zorth=zorth)
+            Rmsep_r.append(rmsep.copy())
+            Yhat_r.append(yhat.copy())
+            CE_r.append(ce.copy())
+        RMS.append(Rmsep_r)
+        YHAT.append(Yhat_r)
+        CEE.append(CE_r)
+    RMS = asarray(RMS)
+    YHAT = asarray(YHAT)
+    CEE = asarray(CEE)
 
 # Significance Hotellings T
-#Wx, Wz = lpls_jk(Xr, Y, Z, aopt, mean_ctr=mean_ctr, xz_alpha=xz_alpha, nsets=nsets)
+calc_qvals = True
-#Ws = W*apply_along_axis(norm, 0, T)
+if not calc_qvals:
-#tsqx = cx_stats.hotelling(Wx, Ws[:,:aopt], alpha=w_alpha)
+    Wx, Wz = pyblm.crossvalidation.lpls_jk(Xr, Y, Z, aopt, center_axis=center_axis, xz_alpha=xz_alpha, nsets=nsets)
-#tsqz = cx_stats.hotelling(Wz, L[:,:aopt], alpha=0)
+    Ws = W*apply_along_axis(norm, 0, T)
+    Ws = Ws[:,:aopt]
+    cal_tsq_x = pyblm.statistics.hotelling(Wx, Ws[:,:aopt], alpha=w_alpha)
+    Ls = L*apply_along_axis(norm, 0, K) 
+    cal_tsq_z = pyblm.statistics.hotelling(Wz, Ls[:,:aopt], alpha=0.01)
 
 # qvals
-cal_tsq_z, pert_tsq_z, cal_tsq_x, pert_tsq_x = cx_stats.lpls_qvals(Xr, Y, Z, aopt=aopt, zx_alpha=xz_alpha, n_iter=n_iter)
 
-qvalz = cx_stats.fdr(cal_tsq_z, pert_tsq_z, 'median')
+if calc_qvals:
-qvalx = cx_stats.fdr(cal_tsq_x, pert_tsq_x, 'median')
+    cal_tsq_z, pert_tsq_z, cal_tsq_x, pert_tsq_x = pyblm.lpls_qvals(Xr, Y, Z, aopt=aopt, zx_alpha=xz_alpha, n_iter=n_iter, nsets=nsets)
+
+    qvalz = pyblm.statistics._fdr(cal_tsq_z, pert_tsq_z, median)
+    qvalx = pyblm.statistics._fdr(cal_tsq_x, pert_tsq_x, median)
 
 
-# p-values, set-enrichment analysis
+    # p-values, set-enrichment analysis
-active_genes_ids = where(qvalx < qval_cutoff)[0]
+    active_genes_ids = where(qvalx < qval_cutoff)[0]
-active_genes = [name for i,name in enumerate(gene_ids) if i in active_genes_ids]
+    active_genes = [name for i,name in enumerate(gene_ids) if i in active_genes_ids]
-active_universe = gene_ids
+    active_universe = gene_ids
-gsea_result, gsea_params= rpy_go.gene_GO_hypergeo_test(genelist=active_genes,universe=active_universe,chip=chip,pval_cutoff=1.0,cond=False,test_direction="over")
+    gsea_result, gsea_params= rpy_go.gene_GO_hypergeo_test(genelist=active_genes,universe=active_universe,chip=chip,pval_cutoff=1.0,cond=False,test_direction="over")
-active_goterms_ids = where(qvalz < qval_cutoff)[0]
+    active_goterms_ids = where(qvalz < qval_cutoff)[0]
-active_goterms = [name for i,name in enumerate(terms) if i in active_goterms_ids]
+    active_goterms = [name for i,name in enumerate(terms) if i in active_goterms_ids]
-
-gsea_t2p = dict(zip(gsea_result['GOBPID'], gsea_result['Pvalue']))
     
+    gsea_t2p = dict(zip(gsea_result['GOBPID'], gsea_result['Pvalue']))
 
 
 
@@ -247,22 +502,35 @@ from scipy import where
 dg = plots_lpls.dag(terms, "bp")
 pos = None
 
-figure(300)
+if calc_qvals:
-subplot(2,1,1)
+    figure(300)
-pos = plots_lpls.plot_dag(dg, node_color=cal_tsq_z, pos=pos, nodelist=terms)
+    subplot(2,1,1)
-subplot(2,1,2)
+    pos = plots_lpls.plot_dag(dg, node_color=cal_tsq_z, pos=pos, nodelist=terms)
-pos = plot_dag(dg, node_color=qvalz, pos=pos, nodelist=terms)
+    ax = gca()
-
+    colorbar(ax.collections[0])
+    xlabel('q values')
+    xticks([])
+    yticks([])
+    subplot(2,1,2)
+    pos = plot_dag(dg, node_color=qvalz, pos=pos, nodelist=terms)
+    ax = gca()
+    colorbar(ax.collections[0])
+    xlabel('T2 values')
+else:
+    figure(300)
+    subplot(2,1,1)
+    pos = plots_lpls.plot_dag(dg, pos=pos, nodelist=terms)
 
 if calc_rmsep:
-    figure(1) #rmsep
+    figure(190) #rmsep
-    bar_w = .2
+    
-    bar_col = 'rgb'*5
+    bar_col = 'rgbcmyk'*2
     m = Y.shape[1]
+    bar_w = 1./(m + 2.)
     for a in range(m):
-        bar(arange(a_max)+a*bar_w+.1, rmsep[:,a], width=bar_w, color=bar_col[a])
+        bar(arange(a_max)+a*bar_w+.1, rmsep[a,:], width=bar_w, color=bar_col[a])
     ylim([rmsep.min()-.05, rmsep.max()+.05]) 
-    title('RMSEP: Y(%s)' %Y.get_name())
+    title('RMSEP: Y(%s)' %DY.get_name())
 
 #figure(2)
 #for a in range(m):
@@ -270,26 +538,28 @@ if calc_rmsep:
 #ylim([class_error.min()-.05, class_error.max()+.05])
 #title('Classification accuracy')
 
-figure(3) # Hyploid correlations
+figure(5) # Hyploid correlations
+pc1 = 2
+pc2 = 3
 tsqz = cal_tsq_z
 tsqx = cal_tsq_x
-tsqz_s = 250*tsqz/tsqz.max()
+tsqz_s = 550*tsqz/tsqz.max()
 td = rpy_go.goterm2desc(terms)
 tlabels = [td[i] for i in terms]
-keep = where(qvalz<0.01)[0]
+#keep = tsqz.argsort()[:100]
-k_Rz = Rz[keep,:]
+#k_Rz = Rz[keep,:]
-k_tsqz_s = tsqz_s[keep]
+#k_tsqz_s = tsqz_s[keep]
-k_tsq = tsqz[keep]
+#k_tsq = tsqz[keep]
-k_tlabels = [name for i,name in enumerate(tlabels) if i in keep]
+#k_tlabels = [name for i,name in enumerate(tlabels) if i in keep]
-plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz_s, c=tsqz, zorder=5, expvar=evz, ax=None,alpha=.5,labels=None)
+plot_corrloads(Rz, pc1=pc1, pc2=pc2, s=tsqz_s, c=tsqz, zorder=6, expvar=evz, ax=None,alpha=.9,labels=None)
 #plot_corrloads(k_Rz, pc1=0, pc2=1, s=k_tsqz_s, c=k_tsqz, zorder=5, expvar=evz, ax=None,alpha=.5,labels=None)
 ax = gca()
-yglabels = DYg.get_identifiers(DYg.get_dim_name()[1], sorted=True)
+ylabels = DYg.get_identifiers(DYg.get_dim_name()[1], sorted=True)
-ylabels = DY.get_identifiers(DY.get_dim_name()[1], sorted=True)
+#ylabels = DYr.get_identifiers(DYr.get_dim_name()[1], sorted=True)
-blabels = yglabels[:]
+#blabels = yglabels[:]
-blabels.append(ylabels[0])
+#blabels.append(ylabels[0])
-plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', marker='s', zorder=5, expvar=evy, ax=ax,labels=None,alpha=.9)
+plot_corrloads(Ry, pc1=pc1, pc2=pc2, s=350, c='g', marker='s', zorder=7, expvar=evy, ax=ax,labels=ylabels,alpha=1.0, drawback=False)
-plot_corrloads(Rx, pc1=0, pc2=1, s=5, c='k', zorder=1, expvar=evx, ax=ax)
+plot_corrloads(Rx, pc1=pc1, pc2=pc2, s=3, c=(.6,.6,.6), alpha=1, zorder=4, expvar=evx, ax=ax, drawback=False, faceted=False)
 
 
 figure(4)
@@ -299,7 +569,7 @@ plot_corrloads(Rx, pc1=0, pc2=1, s=tsqx/2.0, c='b', zorder=5, expvar=evx, ax=ax)
 # title('X correlation')
 subplot(222)
 ax = gca()
-plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', zorder=5, expvar=evy, ax=ax)
+plot_corrloads(Ry, pc1=0, pc2=1, s=250, c='g', zorder=5, expvar=evy, ax=ax)
 #title('Y correlation')
 subplot(223)
 ax = gca()
@@ -312,4 +582,73 @@ plot(evz, 'r', label='Z', linewidth=2)
 legend(loc=2)
 ylabel('Explained variance')
 xlabel('Component')
+xticks((arange(len(evx))), [str(int(i+1)) for i in arange(len(evx))])
 show()
+
+
+figure(19)
+#subplot(1,2,1)
+# RMS : (n_rand_iter, n_alpha, nvarY, a_max)
+# Rmsep : (n_alpha, nvarY, a_max)
+
+rms = RMS[:,:,:,aopt] # looking at solution at aopt
+m_rms = rms.mean(2) # mean over all y-variables
+mm_rms = m_rms.mean(0) # mean over iterations
+std_mrms = m_rms.std(0) # standard deviation over iterations
+
+rms_t = Rmsep[:,:,aopt]
+m_rms_t = rms_t.mean(1)
+xax = arange(mm_rms.shape[0])
+std2_lim_down =  mm_rms - 1.*std_mrms
+std2_lim_up =  mm_rms + 1.*std_mrms
+xx = r_[xax, xax[::-1]] 
+yy = r_[std2_lim_down, std2_lim_up[::-1]]
+fill(xx, yy, fc='.9')
+plot(mm_rms, '--r', lw=1.5, label='Perm. mean')
+plot(std2_lim_down, 'b--')
+plot(std2_lim_up, 'b--', label='Perm. 2*std')
+plot(m_rms_t, 'g', lw=1.5, label='True')
+#c_ylim = ylim()
+#ylim(c_ylim[0], c_ylim[1]-1)
+alpha_ind = linspace(0, Alpha.shape[0]-1, 11)
+xticks(alpha_ind, ['%.1f' %a for a in arange(0,1.01, .1)])
+xlabel(r'$\alpha$')
+ylabel('mean error')
+leg = legend(loc=2)
+# delete fill from legend
+del leg.texts[-1]
+del leg.legendHandles[-1]
+# delete one of the std legends
+del leg.texts[1]
+del leg.legendHandles[1]
+
+klass = True
+
+if klass:
+    figure(20)
+    # subplot(1,2,1)
+    # RMS : (n_rand_iter, n_alpha, nvarY, a_max)
+    # Rmsep : (n_alpha, nvarY, a_max)
+
+    cee = CEE[:,:,aopt,:] # looking at solution at aopt
+    m_cee = cee.mean(-1) # mean over all y-variables
+    mm_cee = m_cee.mean(0) # mean over iterations
+    std_cee = m_cee.std(0) # standard deviation over iterations
+    CE = asarray(CE)
+    cee_t = CE[:,:,aopt]
+    m_cee_t = cee_t.mean(1)
+    xax = arange(mm_cee.shape[0])
+    std2_lim_down =  mm_cee - 2*std_cee
+    std2_lim_up =  mm_cee + 2*std_cee
+    xx = r_[xax, xax[::-1]] 
+    yy = r_[std2_lim_down, std2_lim_up[::-1]]
+    fill(xx, yy, fc='.9')
+    plot(mm_cee, '--r', lw=1.5)
+    plot(std2_lim_down, 'b--')
+    plot(std2_lim_up, 'b--')
+    plot(m_cee_t, 'g', lw=1.5)
+    c_ylim = ylim()
+    ylim = ylim(c_ylim[0], .2)
+    xticks(xax, [str(a)[:3] for a in Alpha])
+    xlabel(r'$\alpha$')
+    ylabel('mean error')