Projects/laydi
Projects
/
laydi
Archived
7
0
Fork 0
Code Issues Pull Requests Activity

Irrelevant play

This commit is contained in:
Arnar Flatberg 2007-08-02 11:19:16 +00:00
parent 2d419a9862
commit 438e7cb918
3 changed files with 38 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
scripts/lpls/lpls.py
View File

@ -365,7 +365,7 @@ def correlation_loadings(D, T, P, test=True):
def cv_lpls(X, Y, Z, a_max=2, nsets=None,alpha=.5): def cv_lpls(X, Y, Z, a_max=2, nsets=None,alpha=.5, mean_ctr=[2,0,1]):
"""Performs crossvalidation to get generalisation error in lpls""" """Performs crossvalidation to get generalisation error in lpls"""
cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=False,index_out=True) cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=False,index_out=True)
k, l = Y.shape k, l = Y.shape
@ -374,7 +374,7 @@ def cv_lpls(X, Y, Z, a_max=2, nsets=None,alpha=.5):
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z, T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z,
a_max=a_max, a_max=a_max,
alpha=alpha, alpha=alpha,
mean_ctr=[2,0,1], mean_ctr=mean_ctr,
verbose=False) verbose=False)
for a in range(a_max): for a in range(a_max):
Yhat[a,ind,:] = b0[a][0][0] + dot(xi, B[a]) Yhat[a,ind,:] = b0[a][0][0] + dot(xi, B[a])
@ -387,7 +387,7 @@ def cv_lpls(X, Y, Z, a_max=2, nsets=None,alpha=.5):
rmsep = sqrt(sep.mean(1)) rmsep = sqrt(sep.mean(1))
return rmsep, Yhat, class_err return rmsep, Yhat, class_err
def jk_lpls(X, Y, Z, a_max, nsets=None, alpha=.5): def jk_lpls(X, Y, Z, a_max, nsets=None, alpha=.5, mean_ctr=[2,0,1]):
cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=False,index_out=False) cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=False,index_out=False)
m, n = X.shape m, n = X.shape
k, l = Y.shape k, l = Y.shape
@ -401,7 +401,7 @@ def jk_lpls(X, Y, Z, a_max, nsets=None, alpha=.5):
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z, T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z,
a_max=a_max, a_max=a_max,
alpha=alpha, alpha=alpha,
mean_ctr=[2,0,1], mean_ctr=mean_ctr,
scale='loads', scale='loads',
verbose=False) verbose=False)
WWx[i,:,:] = W WWx[i,:,:] = W

17
scripts/lpls/rpy_go.py
View File

@ -152,7 +152,18 @@ def parents_dag(go_terms, ontology=['BP']):
edge_dict[nn] = [head] edge_dict[nn] = [head]
return edge_dict return edge_dict
def gene_GO_hypergeo_test(genelist, universe, ontology = ['BP']): def gene_GO_hypergeo_test(genelist,universe="entrezUniverse",ontology="BP",chip = "hgu133a",pval_cutoff=0.01,cond=False,test_direction="over"):
pvals = geneGoHyperGeoTest(entrezGeneIds, lib=None, ontology=ontology[0], universe=universe) #assert(scipy.alltrue([True for i in genelist if i in universe]))
return pvals universeGeneIds=universe
params = rpy.r.new("GOHyperGParams",
geneIds=genelist,
annotation="hgu133a",
ontology=ontology,
pvalueCutoff=pval_cutoff,
conditional=cond,
testDirection=test_direction
)
result = rpy.r.summary(rpy.r.hyperGTest(params))
return rpy.r.summary(result), params

39
scripts/lpls/run_smoker.py
View File

@ -12,7 +12,7 @@ from plots_lpls import plot_corrloads
######## DATA ########## ######## DATA ##########
# full smoker data # full smoker data
DX = dataset.read_ftsv(open("../../data/smokers-full/Xfull.ftsv")) DX = dataset.read_ftsv(open("../../data/smokers-full/Smokers.ftsv"))
DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv")) DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv"))
Y = DY.asarray().astype('d') Y = DY.asarray().astype('d')
# select subset genes by SAM # select subset genes by SAM
@ -32,7 +32,7 @@ print "SAM done"
qq = rpy.r('qobj<-qvalue(sam.out@p.value)') qq = rpy.r('qobj<-qvalue(sam.out@p.value)')
qvals = asarray(qq['qvalues']) qvals = asarray(qq['qvalues'])
# cut off # cut off
cutoff = 2 cutoff = 0.05
index = where(qvals<cutoff)[0] index = where(qvals<cutoff)[0]
# Subset data # Subset data
@ -58,7 +58,7 @@ print "\nNumber of genes: %s" %len(gene_ids)
rpy.r.library("GOSim") rpy.r.library("GOSim")
# Go-term similarity matrix # Go-term similarity matrix
methods = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin") methods = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
meth = methods[3] meth = methods[0]
print "Term-term similarity matrix (method = %s)" %meth print "Term-term similarity matrix (method = %s)" %meth
if meth=="CoutoEnriched": if meth=="CoutoEnriched":
rpy.r('setEnrichmentFactors(alpha=0.1,beta=0.5)') rpy.r('setEnrichmentFactors(alpha=0.1,beta=0.5)')
@ -75,23 +75,23 @@ if len(nanindex)>0:
# Z-matrix # Z-matrix
#Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean) #Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean)
#Z = Z.T #Z = Z.T
Z1 = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat1,go_term_sim="OA",term_sim=meth) Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat1,go_term_sim="OA",term_sim=meth)
#### do another #### do another
meth = methods[4] #meth = methods[4]
rpytmat = 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)
tmat = rpy.r.assign("haha", rpytmat) #tmat = rpy.r.assign("haha", rpytmat)
# check if all terms where found # check if all terms where found
nanindex = where(isnan(tmat[:,0]))[0] #nanindex = where(isnan(tmat[:,0]))[0]
if len(nanindex)>0: #if len(nanindex)>0:
raise valueError("NANs in tmat") # raise valueError("NANs in tmat")
# Z-matrix # Z-matrix
#Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean) #Z, newind = rpy_go.genego_matrix(terms, tmat, gene_ids, terms,func=mean)
#Z = Z.T #Z = Z.T
Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim="OA",term_sim=meth) #Z = rpy_go.genego_sim(gene2goterms,gene_ids,terms,rpytmat,go_term_sim="OA",term_sim=meth)
@ -105,17 +105,20 @@ Xr = X[:,newind]
######## LPLSR ######## ######## LPLSR ########
print "LPLSR ..." print "LPLSR ..."
a_max = 5 a_max = 5
aopt = 2 aopt = 3
alpha=.6 alpha=.4
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(Xr,Y,Z, a_max, alpha) mean_ctr = [2, 0, 1]
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(Xr,Y,Z, a_max,
alpha=alpha,
mean_ctr=mean_ctr)
# Correlation loadings # Correlation loadings
dx,Rx,rssx = correlation_loadings(Xr, T, P) dx,Rx,rssx = correlation_loadings(Xr, T, P)
dx,Ry,rssy = correlation_loadings(Y, T, Q) dx,Ry,rssy = correlation_loadings(Y, T, Q)
cadz,Rz,rssz = correlation_loadings(Z.T, W, L) cadz,Rz,rssz = correlation_loadings(Z.T, W, L)
# Prediction error # Prediction error
rmsep , yhat, class_error = cv_lpls(Xr, Y, Z, a_max, alpha=alpha) rmsep , yhat, class_error = cv_lpls(Xr, Y, Z, a_max, alpha=alpha,mean_ctr=mean_ctr)
alpha_check=False alpha_check=True
if alpha_check: if alpha_check:
Alpha = arange(0.01, 1, .1) Alpha = arange(0.01, 1, .1)
Rmsep,Yhat, CE = [],[],[] Rmsep,Yhat, CE = [],[],[]
@ -127,8 +130,6 @@ if alpha_check:
Rmsep = asarray(Rmsep) Rmsep = asarray(Rmsep)
Yhat = asarray(Yhat) Yhat = asarray(Yhat)
CE = asarray(CE) CE = asarray(CE)
figure(200)
@ -158,7 +159,7 @@ title('Classification accuracy')
figure(3) # Hypoid correlations figure(3) # Hypoid correlations
plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz/10.0, c='b', zorder=5, expvar=evz, ax=None) plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz/10.0, c='b', zorder=5, expvar=evz, ax=None)
ax = gca() ax = gca()
ylabels = DY.get_identifiers('_cat', sorted=True) 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) plot_corrloads(Ry, pc1=0, pc2=1, s=150, c='g', zorder=5, expvar=evy, ax=ax,labels=ylabels)
figure(3) figure(3)