ups

scripts/lpls/plots_lpls.py

@@ -5,7 +5,7 @@ import networkx as nx
 def plot_corrloads(R, pc1=0,pc2=1,s=20, c='b', zorder=5,expvar=None,ax=None,drawback=True, labels=None):
     """ Correlation loading plot."""
 
-    # backgorund
+    # background
     if ax==None or drawback==True:
         radius = 1
         center = (0,0)

scripts/lpls/rpy_go.py

@@ -2,7 +2,7 @@
 import scipy
 import rpy
 silent_eval = rpy.with_mode(rpy.NO_CONVERSION, rpy.r)
-
+import collections
 
 def goterms_from_gene(genelist, ontology='BP', garbage=None):
     """ Returns the go-terms from a specified genelist (Entrez id).
@@ -18,7 +18,7 @@ def goterms_from_gene(genelist, ontology='BP', garbage=None):
               "IDA" : "inferred from direct assay",
               "IEP" : "inferred from expression pattern",
               "IEA" : "inferred from electronic annotation",
-              "TAS" : "traceable author statement",
+               "TAS" : "traceable author statement",
               "NAS" : "non-traceable author statement",
               "ND" : "no biological data available",
               "IC" : "inferred by curator"
@@ -167,3 +167,47 @@ def gene_GO_hypergeo_test(genelist,universe="entrezUniverse",ontology="BP",chip
     result = rpy.r.summary(rpy.r.hyperGTest(params))
     
     return rpy.r.summary(result), params
+
+def data_aff2loc_hgu133a(X, aff_ids, verbose=False):
+    aff_ids = scipy.asarray(aff_ids)
+    if verbose:
+        print "\nNumber of probesets in affy list: %s" %len(aff_ids)
+    import rpy
+    rpy.r.library("hgu133a")
+    trans_table = rpy.r.as_list(rpy.r.hgu133aENTREZID)
+    if verbose:
+        print "Number of entrez ids: %d" %(scipy.asarray(trans_table.values())>0).sum()
+    enz2aff = collections.defaultdict(list)
+    #aff2enz = collections.defaultdict(list)
+    for aff, enz in trans_table.items():
+        if int(enz)>0 and (aff in aff_ids):
+            enz2aff[enz].append(aff)
+            #aff2enz[aff].append(enz)
+    if verbose:
+        print "\nNumber of translated entrez ids: %d" %len(enz2aff)
+    aff2ind = dict(zip(aff_ids, scipy.arange(len(aff_ids))))
+    var_x = X.var(0)
+    new_data = []
+    new_ids = []
+    m = 0
+    s = 0
+    for enz, aff_id_list in enz2aff.items():
+        index = [aff2ind[aff_id] for aff_id in aff_id_list]
+        if len(index)>1:
+            m+=1
+            if verbose:
+                pass
+                #print "\nEntrez id: %s has %d probesets" %(enz, len(index))
+                #print index
+            xsub = X[:,index]
+            choose_this = scipy.argmax(xsub.var(0))
+            new_data.append(xsub[:,choose_this].ravel())
+        else:
+            s+=1
+            new_data.append(X[:,index].ravel())
+        new_ids.append(enz)
+    if verbose:
+        print "Ids with multiple probesets: %d" %m
+        print "Ids with unique probeset: %d" %s
+    X = scipy.asarray(new_data).T
+    return X, new_ids

scripts/lpls/run_smoker.py

@@ -11,41 +11,30 @@ import cx_stats
 from plots_lpls import plot_corrloads
 
 ######## DATA ##########
-# full smoker data
+use_data='uma'
-DX = dataset.read_ftsv(open("../../data/smokers-full/Smokers.ftsv"))
+if use_data=='smoker':
-DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv"))
+    # full smoker data
-Y = DY.asarray().astype('d')
+    DX = dataset.read_ftsv(open("../../data/smokers-full/Smokers.ftsv"))
-# select subset genes by SAM
+    DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv"))
-rpy.r.library("siggenes")
+    Y = DY.asarray().astype('d')
-rpy.r.library("qvalue")
+    gene_ids = DX.get_identifiers('gene_ids', sorted=True)
-data = DX.asarray().T
+elif use_data=='scherf':
-# data = data[:100,:]
+    DX = dataset.read_ftsv(open("../../data/scherf/Scherf.ftsv"))
-rpy.r.assign("data", data)
+    DY = dataset.read_ftsv(open("../../data/scherf/Yd.ftsv"))
-cl = dot(DY.asarray(), diag([1,2,3])).sum(1)
+    Y = DY.asarray().astype('d')
-rpy.r.assign("cl", cl)
+    gene_ids = DX.get_identifiers('gene_ids', sorted=True)
-rpy.r.assign("B", 20)
+elif use_data=='staunton':
-# Perform a SAM analysis.
+    pass
-print "Starting SAM"
+elif use_data=='uma':
-sam = rpy.r('sam.out<-sam(data=data,cl=cl,B=B,rand=123)')
+    DX = dataset.read_ftsv(open("../../data/uma/X133.ftsv"))
-print "SAM done"
+    DY = dataset.read_ftsv(open("../../data/uma/Yg133.ftsv"))
-# Compute the q-values of the genes.
+    Y = DY.asarray().astype('d')
-qq = rpy.r('qobj<-qvalue(sam.out@p.value)')
+    gene_ids = DX.get_identifiers('gene_ids', sorted=True)
-qvals = asarray(qq['qvalues'])
-# cut off
-cutoff = 0.05
-index = where(qvals<cutoff)[0]
 
-# Subset data
+# Use only subset defined on GO
-X = DX.asarray()
+ontology = 'BP'
-Xr = X[:,index]
-gene_ids = DX.get_identifiers('gene_ids', index)
-print "\nWorking on subset with %s genes " %len(gene_ids)
-#gene2ind = {}
-#for i, gene in enumerate(gene_ids):
-#    gene2ind[gene] = i
-
-### Build GO data ####
 print "\n\nFiltering genes by Go terms "
+
 gene2goterms = rpy_go.goterms_from_gene(gene_ids)
 all_terms = set()
 for t in gene2goterms.values():
@@ -55,61 +44,72 @@ print "\nNumber of go-terms: %s" %len(terms)
 # update genelist
 gene_ids = gene2goterms.keys()
 print "\nNumber of genes: %s" %len(gene_ids)
+
+# use subset based on SAM or IQR
+subset = 'm'
+if subset=='sam':
+    # select subset genes by SAM
+    rpy.r.library("siggenes")
+    rpy.r.library("qvalue")
+    data = DX.asarray().T
+    # data = data[:100,:]
+    rpy.r.assign("data", data)
+    cl = dot(DY.asarray(), diag([1,2,3])).sum(1)
+    rpy.r.assign("cl", cl)
+    rpy.r.assign("B", 20)
+    # Perform a SAM analysis.
+    print "Starting SAM"
+    sam = rpy.r('sam.out<-sam(data=data,cl=cl,B=B,rand=123)')
+    print "SAM done"
+    # Compute the q-values of the genes.
+    qq = rpy.r('qobj<-qvalue(sam.out@p.value)')
+    qvals = asarray(qq['qvalues'])
+    # cut off
+    cutoff = 0.001
+    index = where(qvals<cutoff)[0]
+    # Subset data
+    X = DX.asarray()
+    #Xr = X[:,index]
+    gene_ids = DX.get_identifiers('gene_ids', index)
+    print "\nWorking on subset with %s genes " %len(gene_ids)
+else:
+    # noimp (smoker data is prefiltered)
+    pass
+
+
 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)')
 print "\nCalculating term-term similarity matrix"
 
-rpytmat1 = rpy.with_mode(rpy.NO_CONVERSION, rpy.r.getTermSim)(terms, method=meth,verbose=False)
+rpytmat = rpy.with_mode(rpy.NO_CONVERSION, rpy.r.getTermSim)(terms, method=meth,verbose=False)
-tmat1 = rpy.r.assign("haha", rpytmat1)
+tmat = rpy.r.assign("haha", rpytmat)
-
+print "\n Calculating Z matrix"
-# check if all terms where found
+Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim="OA",term_sim=meth)
-nanindex = where(isnan(tmat1[:,0]))[0]
-if len(nanindex)>0:
-    raise valueError("NANs in tmat")
-
-# Z-matrix
-#Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean)
-#Z = Z.T
-Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat1,go_term_sim="OA",term_sim=meth)
-
-
-#### do another
-#meth = methods[4]
-#rpytmat = rpy.with_mode(rpy.NO_CONVERSION, rpy.r.getTermSim)(terms, method=meth,verbose=False)
-#tmat = rpy.r.assign("haha", rpytmat)
-
-# check if all terms where found
-#nanindex = where(isnan(tmat[:,0]))[0]
-#if len(nanindex)>0:
-#    raise valueError("NANs in tmat")
-
-# Z-matrix
-#Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean)
-#Z = Z.T
-#Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim="OA",term_sim=meth)
-
-
 
 # update data (X) matrix
-#newind = [gene2ind[gene] for gene in gene_ids]
 newind = DX.get_indices('gene_ids', gene_ids)
-Xr = X[:,newind]
+Xr = DX.asarray()[:,newind]
-#new_gene_ids = asarray(gene_ids)[newind]
 
 
 ######## LPLSR ########
 print "LPLSR ..."
 a_max = 5
 aopt = 3
-alpha=.4
+xz_alpha = .5
-mean_ctr = [2, 0, 1]
+w_alpha = .1
+mean_ctr = [2, 0, 2]
+
+# standardize Z?
+sdtz = False
+if sdtz:
+    Z = Z/Z.std(0)
+
 T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(Xr,Y,Z, a_max,
-                                                        alpha=alpha,
+                                                        alpha=xz_alpha,
                                                         mean_ctr=mean_ctr)
 
 # Correlation loadings
@@ -117,24 +117,22 @@ dx,Rx,rssx = correlation_loadings(Xr, T, P)
 dx,Ry,rssy = correlation_loadings(Y, T, Q)
 cadz,Rz,rssz = correlation_loadings(Z.T, W, L)
 # Prediction error
-rmsep , yhat, class_error = cv_lpls(Xr, Y, Z, a_max, alpha=alpha,mean_ctr=mean_ctr)
+rmsep , yhat, class_error = cv_lpls(Xr, Y, Z, a_max, alpha=xz_alpha,mean_ctr=mean_ctr)
-alpha_check=True
+alpha_check=False
 if alpha_check:
     Alpha = arange(0.01, 1, .1)
     Rmsep,Yhat, CE = [],[],[]
     for a in Alpha:
-        rmsep , yhat, ce = cv_lpls(Xr, Y, Z, a_max, alpha=alpha)
+        rmsep , yhat, ce = cv_lpls(Xr, Y, Z, a_max, alpha=xz_alpha,mean_ctr=mean_ctr)
         Rmsep.append(rmsep)
         Yhat.append(yhat)
         CE.append(ce)
     Rmsep = asarray(Rmsep)
     Yhat = asarray(Yhat)
     CE = asarray(CE)
-    
-
 
 # Significance Hotellings T
-Wx, Wz, Wy, = jk_lpls(Xr, Y, Z, aopt)
+Wx, Wz, Wy, = jk_lpls(Xr, Y, Z, aopt, mean_ctr=mean_ctr,alpha=w_alpha)
 Ws = W*apply_along_axis(norm, 0, T)
 tsqx = cx_stats.hotelling(Wx, Ws[:,:aopt])
 tsqz = cx_stats.hotelling(Wz, L[:,:aopt])
@@ -157,7 +155,8 @@ ylim([class_error.min()-.05, class_error.max()+.05])
 title('Classification accuracy')
 
 figure(3) # Hypoid correlations
-plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz/10.0, c='b', zorder=5, expvar=evz, ax=None)
+tsqz_s = 250*tsqz/tsqz.max()
+plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz_s, c='b', zorder=5, expvar=evz, ax=None)
 ax = gca()
 ylabels = DY.get_identifiers('_status', sorted=True)
 plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', zorder=5, expvar=evy, ax=ax,labels=ylabels)