Initial import
This commit is contained in:
parent
dd04e28a62
commit
7ee7aa968a
|
|
@ -0,0 +1,409 @@
|
||||||
|
import sys
|
||||||
|
from pylab import *
|
||||||
|
import matplotlib
|
||||||
|
from scipy import *
|
||||||
|
from scipy.linalg import inv,norm
|
||||||
|
|
||||||
|
sys.path.append("/home/flatberg/fluents/fluents/lib")
|
||||||
|
import select_generators
|
||||||
|
sys.path.remove("/home/flatberg/fluents/fluents/lib")
|
||||||
|
|
||||||
|
|
||||||
|
def nipals_lpls(X, Y, Z, a_max, alpha=.7, mean_ctr=[2, 0, 1], verbose=True):
|
||||||
|
""" L-shaped Partial Least Sqaures Regression by the nipals algorithm.
|
||||||
|
|
||||||
|
(X!Z)->Y
|
||||||
|
:input:
|
||||||
|
X : data matrix (m, n)
|
||||||
|
Y : data matrix (m, l)
|
||||||
|
Z : data matrix (n, o)
|
||||||
|
|
||||||
|
:output:
|
||||||
|
T : X-scores
|
||||||
|
W : X-weights/Z-weights
|
||||||
|
P : X-loadings
|
||||||
|
Q : Y-loadings
|
||||||
|
U : X-Y relation
|
||||||
|
L : Z-scores
|
||||||
|
K : Z-loads
|
||||||
|
B : Regression coefficients X->Y
|
||||||
|
b0: Regression coefficient intercept
|
||||||
|
evx : X-explained variance
|
||||||
|
evy : Y-explained variance
|
||||||
|
evz : Z-explained variance
|
||||||
|
|
||||||
|
:Notes:
|
||||||
|
|
||||||
|
"""
|
||||||
|
if mean_ctr:
|
||||||
|
xctr, yctr, zctr = mean_ctr
|
||||||
|
X, mnX = center(X, xctr)
|
||||||
|
Y, mnY = center(Y, xctr)
|
||||||
|
Z, mnZ = center(Z, zctr)
|
||||||
|
|
||||||
|
varX = pow(X, 2).sum()
|
||||||
|
varY = pow(Y, 2).sum()
|
||||||
|
varZ = pow(Z, 2).sum()
|
||||||
|
|
||||||
|
m, n = X.shape
|
||||||
|
k, l = Y.shape
|
||||||
|
u, o = Z.shape
|
||||||
|
|
||||||
|
# initialize
|
||||||
|
U = empty((k, a_max))
|
||||||
|
Q = empty((l, a_max))
|
||||||
|
T = empty((m, a_max))
|
||||||
|
W = empty((n, a_max))
|
||||||
|
P = empty((n, a_max))
|
||||||
|
K = empty((o, a_max))
|
||||||
|
L = empty((u, a_max))
|
||||||
|
B = empty((a_max, n, l))
|
||||||
|
b0 = empty((a_max, m, l))
|
||||||
|
var_x = empty((a_max,))
|
||||||
|
var_y = empty((a_max,))
|
||||||
|
var_z = empty((a_max,))
|
||||||
|
|
||||||
|
for a in range(a_max):
|
||||||
|
if verbose:
|
||||||
|
print "\n Working on comp. %s" %a
|
||||||
|
u = Y[:,:1]
|
||||||
|
diff = 1
|
||||||
|
MAX_ITER = 100
|
||||||
|
lim = 1e-7
|
||||||
|
niter = 0
|
||||||
|
while (diff>lim and niter<MAX_ITER):
|
||||||
|
niter += 1
|
||||||
|
u1 = u.copy()
|
||||||
|
w = dot(X.T, u)
|
||||||
|
w = w/sqrt(dot(w.T, w))
|
||||||
|
l = dot(Z, w)
|
||||||
|
k = dot(Z.T, l)
|
||||||
|
k = k/sqrt(dot(k.T, k))
|
||||||
|
w = alpha*k + (1-alpha)*w
|
||||||
|
w = w/sqrt(dot(w.T, w))
|
||||||
|
t = dot(X, w)
|
||||||
|
c = dot(Y.T, t)
|
||||||
|
c = c/sqrt(dot(c.T, c))
|
||||||
|
u = dot(Y, c)
|
||||||
|
diff = abs(u1 - u).max()
|
||||||
|
if verbose:
|
||||||
|
print "Converged after %s iterations" %niter
|
||||||
|
tt = dot(t.T, t)
|
||||||
|
p = dot(X.T, t)/tt
|
||||||
|
q = dot(Y.T, t)/tt
|
||||||
|
l = dot(Z, w)
|
||||||
|
U[:,a] = u.ravel()
|
||||||
|
W[:,a] = w.ravel()
|
||||||
|
P[:,a] = p.ravel()
|
||||||
|
T[:,a] = t.ravel()
|
||||||
|
Q[:,a] = q.ravel()
|
||||||
|
L[:,a] = l.ravel()
|
||||||
|
K[:,a] = k.ravel()
|
||||||
|
|
||||||
|
X = X - dot(t, p.T)
|
||||||
|
Y = Y - dot(t, q.T)
|
||||||
|
Z = (Z.T - dot(w, l.T)).T
|
||||||
|
|
||||||
|
var_x[a] = pow(X, 2).sum()
|
||||||
|
var_y[a] = pow(Y, 2).sum()
|
||||||
|
var_z[a] = pow(Z, 2).sum()
|
||||||
|
B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
|
||||||
|
b0[a] = mnY - dot(mnX, B[a])
|
||||||
|
# variance explained
|
||||||
|
evx = 100.0*(1 - var_x/varX)
|
||||||
|
evy = 100.0*(1 - var_y/varY)
|
||||||
|
evz = 100.0*(1 - var_z/varZ)
|
||||||
|
|
||||||
|
return T, W, P, Q, U, L, K, B, b0, evx, evy, evz
|
||||||
|
|
||||||
|
def svd_lpls(X, Y, Z, a_max, alpha=.7, mean_ctr=[2, 0, 1], verbose=True):
|
||||||
|
"""
|
||||||
|
NB: In the works ...
|
||||||
|
L-shaped Partial Least Sqaures Regression by the svd algorithm.
|
||||||
|
|
||||||
|
(X!Z)->Y
|
||||||
|
:input:
|
||||||
|
X : data matrix (m, n)
|
||||||
|
Y : data matrix (m, l)
|
||||||
|
Z : data matrix (n, o)
|
||||||
|
|
||||||
|
:output:
|
||||||
|
T : X-scores
|
||||||
|
W : X-weights/Z-weights
|
||||||
|
P : X-loadings
|
||||||
|
Q : Y-loadings
|
||||||
|
U : X-Y relation
|
||||||
|
L : Z-scores
|
||||||
|
K : Z-loads
|
||||||
|
B : Regression coefficients X->Y
|
||||||
|
b0: Regression coefficient intercept
|
||||||
|
evx : X-explained variance
|
||||||
|
evy : Y-explained variance
|
||||||
|
evz : Z-explained variance
|
||||||
|
|
||||||
|
:Notes:
|
||||||
|
Not quite there ,,,,,,,,,,,,,,
|
||||||
|
|
||||||
|
"""
|
||||||
|
if mean_ctr:
|
||||||
|
xctr, yctr, zctr = mean_ctr
|
||||||
|
X, mnX = center(X, xctr)
|
||||||
|
Y, mnY = center(Y, xctr)
|
||||||
|
Z, mnZ = center(Z, zctr)
|
||||||
|
|
||||||
|
varX = pow(X, 2).sum()
|
||||||
|
varY = pow(Y, 2).sum()
|
||||||
|
varZ = pow(Z, 2).sum()
|
||||||
|
|
||||||
|
m, n = X.shape
|
||||||
|
k, l = Y.shape
|
||||||
|
u, o = Z.shape
|
||||||
|
|
||||||
|
# initialize
|
||||||
|
U = empty((k, a_max))
|
||||||
|
Q = empty((l, a_max))
|
||||||
|
T = empty((m, a_max))
|
||||||
|
W = empty((n, a_max))
|
||||||
|
P = empty((n, a_max))
|
||||||
|
K = empty((o, a_max))
|
||||||
|
L = empty((u, a_max))
|
||||||
|
var_x = empty((a_max,))
|
||||||
|
var_y = empty((a_max,))
|
||||||
|
var_z = empty((a_max,))
|
||||||
|
|
||||||
|
for a in range(a_max):
|
||||||
|
if verbose:
|
||||||
|
print "\n Working on comp. %s" %a
|
||||||
|
xyz = dot(dot(Z,X.T),Y)
|
||||||
|
u,s,vt = linalg.svd(xyz, 0)
|
||||||
|
w = u[:,o]
|
||||||
|
t = dot(X, w)
|
||||||
|
tt = dot(t.T, t)
|
||||||
|
p = dot(X.T, t)/tt
|
||||||
|
q = dot(Y.T, t)/tt
|
||||||
|
l = dot(Z.T, w)
|
||||||
|
W[:,a] = w.ravel()
|
||||||
|
P[:,a] = p.ravel()
|
||||||
|
T[:,a] = t.ravel()
|
||||||
|
Q[:,a] = q.ravel()
|
||||||
|
L[:,a] = l.ravel()
|
||||||
|
K[:,a] = k.ravel()
|
||||||
|
X = X - dot(t, p.T)
|
||||||
|
Y = Y - dot(t, q.T)
|
||||||
|
Z = (Z.T - dot(w, l.T)).T
|
||||||
|
var_x[a] = pow(X, 2).sum()
|
||||||
|
var_y[a] = pow(Y, 2).sum()
|
||||||
|
var_z[a] = pow(Z, 2).sum()
|
||||||
|
B = dot(dot(W, inv(dot(P.T, W))), Q.T)
|
||||||
|
b0 = mnY - dot(mnX, B)
|
||||||
|
# variance explained
|
||||||
|
evx = 100.0*(1 - var_x/varX)
|
||||||
|
evy = 100.0*(1 - var_y/varY)
|
||||||
|
evz = 100.0*(1 - var_z/varZ)
|
||||||
|
return T, W, P, Q, U, L, K, B, b0, evx, evy, evz
|
||||||
|
|
||||||
|
|
||||||
|
def lplsr(X, Y, Z, a_max, mean_ctr=[2,0,1]):
|
||||||
|
""" Haralds LPLS.
|
||||||
|
"""
|
||||||
|
if mean_ctr!=None:
|
||||||
|
xctr, yctr, zctr = mean_ctr
|
||||||
|
X, mnX = center(X, xctr)
|
||||||
|
Y, mnY = center(Y, yctr)
|
||||||
|
Z, mnZ = center(Z, zctr)
|
||||||
|
|
||||||
|
varX = pow(X, 2).sum()
|
||||||
|
varY = pow(Y, 2).sum()
|
||||||
|
varZ = pow(Z, 2).sum()
|
||||||
|
m, n = X.shape
|
||||||
|
k, l = Y.shape
|
||||||
|
u, o = Z.shape
|
||||||
|
|
||||||
|
# initialize
|
||||||
|
Wy = empty((l, a_max))
|
||||||
|
Py = empty((l, a_max))
|
||||||
|
Ty = empty((m, a_max))
|
||||||
|
Tz = empty((o, a_max))
|
||||||
|
Wz = empty((u, a_max))
|
||||||
|
Pz = empty((u, a_max))
|
||||||
|
var_x = empty((a_max,))
|
||||||
|
var_y = empty((a_max,))
|
||||||
|
var_z = empty((a_max,))
|
||||||
|
|
||||||
|
# residuals
|
||||||
|
Ey = Y.copy()
|
||||||
|
Ez = Z.copy()
|
||||||
|
Ex = X.copy()
|
||||||
|
for i in range(a_max):
|
||||||
|
YtXZ = dot(Ey.T, dot(X, Ez.T))
|
||||||
|
U, S, V = linalg.svd(YtXZ)
|
||||||
|
wy = U[:,0]
|
||||||
|
print wy
|
||||||
|
wz = V[0,:]
|
||||||
|
ty = dot(Ey, wy)
|
||||||
|
tz = dot(Ez.T, wz)
|
||||||
|
py = dot(Ey.T, ty)/dot(ty.T,ty)
|
||||||
|
pz = dot(Ez, tz)/dot(tz.T,tz)
|
||||||
|
Wy[:,i] = wy
|
||||||
|
Wz[:,i] = wz
|
||||||
|
Ty[:,i] = ty
|
||||||
|
Tz[:,i] = tz
|
||||||
|
Py[:,i] = py
|
||||||
|
Pz[:,i] = pz
|
||||||
|
Ey = Ey - outer(ty, py.T)
|
||||||
|
Ez = (Ez.T - outer(tz, pz.T)).T
|
||||||
|
var_y[i] = pow(Ey, 2).sum()
|
||||||
|
var_z[i] = pow(Ez, 2).sum()
|
||||||
|
|
||||||
|
tyd = apply_along_axis(norm, 0, Ty)
|
||||||
|
tzd = apply_along_axis(norm, 0, Tz)
|
||||||
|
Tyu = Ty/tyd
|
||||||
|
Tzu = Tz/tzd
|
||||||
|
C = dot(dot(Tyu.T, X), Tzu)
|
||||||
|
for i in range(a_max):
|
||||||
|
Ex = Ex - dot(dot(Ty[:,:i+1],C[:i+1,:i+1]), Tz[:,:i+1].T)
|
||||||
|
var_x[i] = pow(Ex,2).sum()
|
||||||
|
# variance explained
|
||||||
|
print "var_x:"
|
||||||
|
print var_x
|
||||||
|
print "varX total:"
|
||||||
|
print varX
|
||||||
|
|
||||||
|
evx = 100.0*(1 - var_x/varX)
|
||||||
|
evy = 100.0*(1 - var_y/varY)
|
||||||
|
evz = 100.0*(1 - var_z/varZ)
|
||||||
|
|
||||||
|
return Ty, Tz, Wy, Wz, Py, Pz, C, Ey, Ez, Ex, evx, evy, evz
|
||||||
|
|
||||||
|
def bifpls(X, Y, Z, a_max, alpha):
|
||||||
|
"""Swedssihsh LPLS by nipals.
|
||||||
|
"""
|
||||||
|
u = X[:,0]
|
||||||
|
Ey = Y.copy()
|
||||||
|
Ez = Z.copy()
|
||||||
|
for i in range(100):
|
||||||
|
w = dot(X.T,u)
|
||||||
|
w = w/vnorm(w)
|
||||||
|
t = dot(X, w)
|
||||||
|
q = dot(Ey, t.T)/dot(t.T,t)
|
||||||
|
qnorm = vnorm(q)
|
||||||
|
q = q/qnorm
|
||||||
|
v = dot(Ez, q)
|
||||||
|
s = dot(Ez.T, v)/dot(v.T,v)
|
||||||
|
v = v*vnorm(s)
|
||||||
|
s = s/vnorm(s)
|
||||||
|
c = qnorm*(alpha*q + (1-alpha)*s)
|
||||||
|
u = dot(Ey, c)/dot(s.T,s)
|
||||||
|
p = dot(X.T, t)/dot(t.T,t)
|
||||||
|
v2 = dot(Ez, s)/dot(s.T,s)
|
||||||
|
Ey = Ey - dot(t, p.T)
|
||||||
|
Ez = Ez - dot(v2, c.T)
|
||||||
|
# variance explained
|
||||||
|
evx = 100.0*(1 - var_x/varX)
|
||||||
|
evy = 100.0*(1 - var_y/varY)
|
||||||
|
evz = 100.0*(1 - var_z/varZ)
|
||||||
|
|
||||||
|
def center(a, axis):
|
||||||
|
# 0 = col center, 1 = row center, 2 = double center
|
||||||
|
# -1 = nothing
|
||||||
|
if axis==-1:
|
||||||
|
return a
|
||||||
|
elif axis==0:
|
||||||
|
mn = a.mean(0)
|
||||||
|
return a - mn, mn
|
||||||
|
elif axis==1:
|
||||||
|
mn = a.mean(1)[:,newaxis]
|
||||||
|
return a - mn , mn
|
||||||
|
elif axis==2:
|
||||||
|
mn = a.mean(0) + a.mean(1)[:,newaxis] - a.mean()
|
||||||
|
return a - mn, mn
|
||||||
|
else:
|
||||||
|
raise IOError("input error: axis must be in [-1,0,1,2]")
|
||||||
|
|
||||||
|
def correlation_loadings(D, T, P, test=True):
|
||||||
|
""" Returns correlation loadings.
|
||||||
|
|
||||||
|
:input:
|
||||||
|
- D: [nsamps, nvars], data (non-centered data)
|
||||||
|
- T: [nsamps, a_max], Scores
|
||||||
|
- P: [nvars, a_max], Loadings
|
||||||
|
:ouput:
|
||||||
|
- Rloads: [nvars, a_max], Correlation loadings
|
||||||
|
- rmseVars: [nvars], scaling coeff. for each var in D
|
||||||
|
|
||||||
|
:notes:
|
||||||
|
- FIXME: Calculation is not valid .... using corrceof instead
|
||||||
|
"""
|
||||||
|
nsamps, nvars = D.shape
|
||||||
|
nsampsT, a_max = T.shape
|
||||||
|
nvarsP, a_maxP = P.shape
|
||||||
|
if nsamps!=nsampsT: raise IOError("D/T mismatch")
|
||||||
|
if a_max!=a_maxP: raise IOError("a_max mismatch")
|
||||||
|
if nvars!=nvarsP: raise IOError("D/P mismatch")
|
||||||
|
|
||||||
|
#init
|
||||||
|
Rloads = empty((nvars, a_max), 'd')
|
||||||
|
stdvar = stats.std(D, 0)
|
||||||
|
rmseVars = sqrt(nsamps-1)*stdvar
|
||||||
|
|
||||||
|
# center
|
||||||
|
D = D - D.mean(0)
|
||||||
|
TT = diag(dot(T.T, T))
|
||||||
|
sTT = sqrt(TT)
|
||||||
|
for a in range(a_max):
|
||||||
|
Rloads[:,a] = sTT[a]*P[:,a]/rmseVars
|
||||||
|
R = empty_like(Rloads)
|
||||||
|
for a in range(a_max):
|
||||||
|
for k in range(nvars):
|
||||||
|
r = corrcoef(D[:,k], T[:,a])
|
||||||
|
R[k,a] = r[0,1]
|
||||||
|
#Rloads = R
|
||||||
|
return Rloads, R, rmseVars
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def cv_lpls(X, Y, Z, a_max=2, nsets=None,alpha=.5):
|
||||||
|
"""Performs crossvalidation to get generalisation error in lpls"""
|
||||||
|
cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=True,index_out=True)
|
||||||
|
k, l = Y.shape
|
||||||
|
Yhat = empty((a_max,k,l), 'd')
|
||||||
|
for i, (xcal,xi,ycal,yi,ind) in enumerate(cv_iter):
|
||||||
|
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z,
|
||||||
|
a_max=a_max,
|
||||||
|
alpha=alpha,
|
||||||
|
mean_ctr=[0,0,1],
|
||||||
|
verbose=False)
|
||||||
|
for a in range(a_max):
|
||||||
|
Yhat[a,ind,:] = b0[a][0][0] + dot(xi, B[a])
|
||||||
|
Yhat_class = zeros_like(Yhat)
|
||||||
|
for a in range(a_max):
|
||||||
|
for i in range(k):
|
||||||
|
Yhat_class[a,i,argmax(Yhat[a,i,:])]=1.0
|
||||||
|
class_err = 100*((Yhat_class+Y)==2).sum(1)/Y.sum(0).astype('d')
|
||||||
|
sep = (Y - Yhat)**2
|
||||||
|
rmsep = sqrt(sep.mean(1))
|
||||||
|
return rmsep, Yhat, class_err
|
||||||
|
|
||||||
|
def jk_lpls(X, Y, Z, a_max, nsets=None, alpha=.5):
|
||||||
|
cv_iter = select_generators.pls_gen(X, Y, n_blocks=nsets,center=True,index_out=False)
|
||||||
|
m, n = X.shape
|
||||||
|
k, l = Y.shape
|
||||||
|
o, p = Z.shape
|
||||||
|
if nsets==None:
|
||||||
|
nsets = m
|
||||||
|
WWx = empty((nsets, n, a_max), 'd')
|
||||||
|
WWz = empty((nsets, o, a_max), 'd')
|
||||||
|
WWy = empty((nsets, l, a_max), 'd')
|
||||||
|
for i, (xcal,xi,ycal,yi) in enumerate(cv_iter):
|
||||||
|
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(xcal,ycal,Z,
|
||||||
|
a_max=a_max,
|
||||||
|
alpha=alpha,
|
||||||
|
mean_ctr=[0,0,1],
|
||||||
|
verbose=False)
|
||||||
|
WWx[i,:,:] = W
|
||||||
|
WWz[i,:,:] = L
|
||||||
|
WWy[i,:,:] = Q
|
||||||
|
print "Q"
|
||||||
|
print Q
|
||||||
|
|
||||||
|
return WWx, WWz, WWy
|
||||||
|
|
@ -0,0 +1,38 @@
|
||||||
|
import pylab
|
||||||
|
import matplotlib
|
||||||
|
|
||||||
|
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
|
||||||
|
if ax==None or drawback==True:
|
||||||
|
radius = 1
|
||||||
|
center = (0,0)
|
||||||
|
c100 = matplotlib.patches.Circle(center,radius=radius,
|
||||||
|
facecolor='gray',
|
||||||
|
alpha=.1,
|
||||||
|
zorder=1)
|
||||||
|
c50 = matplotlib.patches.Circle(center, radius=radius/2.0,
|
||||||
|
facecolor='gray',
|
||||||
|
alpha=.1,
|
||||||
|
zorder=2)
|
||||||
|
ax = pylab.gca()
|
||||||
|
ax.add_patch(c100)
|
||||||
|
ax.add_patch(c50)
|
||||||
|
ax.axhline(lw=1.5,color='k')
|
||||||
|
ax.axvline(lw=1.5,color='k')
|
||||||
|
|
||||||
|
# corrloads
|
||||||
|
ax.scatter(R[:,pc1], R[:,pc2], s=s, c=c,zorder=zorder)
|
||||||
|
ax.set_xlim([-1,1])
|
||||||
|
ax.set_ylim([-1,1])
|
||||||
|
if expvar!=None:
|
||||||
|
xstring = "Comp: %d expl.var: %.1f " %(pc1+1, expvar[pc1])
|
||||||
|
pylab.xlabel(xstring)
|
||||||
|
ystring = "Comp: %d expl.var.: %.1f " %(pc2+1, expvar[pc2])
|
||||||
|
pylab.ylabel(ystring)
|
||||||
|
if labels:
|
||||||
|
assert(len(labels)==R.shape[0])
|
||||||
|
for name, r in zip(labels, R):
|
||||||
|
ax.text(r[pc1], r[pc2], " " + name)
|
||||||
|
#pylab.show()
|
||||||
|
|
@ -0,0 +1,110 @@
|
||||||
|
""" Module for Gene ontology related functions called in R"""
|
||||||
|
import scipy
|
||||||
|
|
||||||
|
import rpy
|
||||||
|
silent_eval = rpy.with_mode(rpy.NO_CONVERSION, rpy.r)
|
||||||
|
|
||||||
|
def get_term_sim(termlist, method = "JiangConrath", verbose=False):
|
||||||
|
"""Returns the similariy matrix between go-terms.
|
||||||
|
|
||||||
|
Arguments:
|
||||||
|
termlist: character vector of GO terms
|
||||||
|
method: one of
|
||||||
|
("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
|
||||||
|
verbose: print out various information or not
|
||||||
|
"""
|
||||||
|
_methods = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
|
||||||
|
assert(method in _methods)
|
||||||
|
assert(termlist[0][:2]=='GO')
|
||||||
|
rpy.r.library("GOSim")
|
||||||
|
return rpy.r.getTermSim(termlist, method = method, verbose = verbose)
|
||||||
|
|
||||||
|
def get_gene_sim(genelist, similarity='OA',
|
||||||
|
distance="Resnick"):
|
||||||
|
rpy.r.library("GOSim")
|
||||||
|
rpy.r.assign("ids", genelist)
|
||||||
|
silent_eval('a<-getGeneSim(ids)', verbose=FALSE)
|
||||||
|
|
||||||
|
def goterms_from_gene(genelist, ontology=['BP'], garbage = ['IEA', 'ISS', 'ND']):
|
||||||
|
""" Returns the go-terms from a specified genelist (Entrez id).
|
||||||
|
|
||||||
|
"""
|
||||||
|
rpy.r.library("GO")
|
||||||
|
_CODES = {"IMP" : "inferred from mutant phenotype",
|
||||||
|
"IGI" : "inferred from genetic interaction",
|
||||||
|
"IPI" :"inferred from physical interaction",
|
||||||
|
"ISS" : "inferred from sequence similarity",
|
||||||
|
"IDA" : "inferred from direct assay",
|
||||||
|
"IEP" : "inferred from expression pattern",
|
||||||
|
"IEA" : "inferred from electronic annotation",
|
||||||
|
"TAS" : "traceable author statement",
|
||||||
|
"NAS" : "non-traceable author statement",
|
||||||
|
"ND" : "no biological data available",
|
||||||
|
"IC" : "inferred by curator"
|
||||||
|
}
|
||||||
|
_ONTOLOGIES = ['BP', 'CC', 'MF']
|
||||||
|
assert(scipy.all([(code in _CODES) for code in garbage]))
|
||||||
|
assert(scipy.all([(ont in _ONTOLOGIES) for ont in ontology]))
|
||||||
|
|
||||||
|
goterms = {}
|
||||||
|
for gene in genelist:
|
||||||
|
goterms[gene] = []
|
||||||
|
info = rpy.r('GOENTREZID2GO[["' + str(gene) + '"]]')
|
||||||
|
#print info
|
||||||
|
if info:
|
||||||
|
for term, desc in info.items():
|
||||||
|
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):
|
||||||
|
ngenes = len(gene_ids)
|
||||||
|
nterms = len(term_ids)
|
||||||
|
gene2indx = {}
|
||||||
|
for i,id in enumerate(gene_ids):
|
||||||
|
gene2indx[id]=i
|
||||||
|
term2indx = {}
|
||||||
|
for i,id in enumerate(term_ids):
|
||||||
|
term2indx[id]=i
|
||||||
|
#G = scipy.empty((nterms, ngenes),'d')
|
||||||
|
G = []
|
||||||
|
newindex = []
|
||||||
|
for gene, terms in goterms.items():
|
||||||
|
g_ind = gene2indx[gene]
|
||||||
|
if len(terms)>0:
|
||||||
|
t_ind = []
|
||||||
|
newindex.append(g_ind)
|
||||||
|
for term in terms:
|
||||||
|
if term2indx.has_key(term): t_ind.append(term2indx[term])
|
||||||
|
print t_ind
|
||||||
|
subsim = tmat[t_ind, :]
|
||||||
|
gene_vec = scipy.apply_along_axis(func, 0, subsim)
|
||||||
|
G.append(gene_vec)
|
||||||
|
|
||||||
|
return scipy.asarray(G), newindex
|
||||||
|
|
||||||
|
def goterm2desc(gotermlist):
|
||||||
|
"""Returns the go-terms description keyed by go-term
|
||||||
|
"""
|
||||||
|
rpy.r.library("GO")
|
||||||
|
term2desc = {}
|
||||||
|
for term in gotermlist:
|
||||||
|
try:
|
||||||
|
desc = rpy.r('Term(GOTERM[["' +str(term)+ '"]])')
|
||||||
|
term2desc[str(term)] = desc
|
||||||
|
except:
|
||||||
|
raise Warning("Description not found for %s\n Mapping incomplete" %term)
|
||||||
|
return term2desc
|
||||||
|
|
||||||
|
def parents_dag(go_terms, ontology=['BP']):
|
||||||
|
""" Returns a list of lists representation of a GO DAG parents of goterms."""
|
||||||
|
try:
|
||||||
|
rpy.r.library("GOstats")
|
||||||
|
except:
|
||||||
|
raise ImportError, "Gostats"
|
||||||
|
assert(go_terms[0][:3]=='GO:')
|
||||||
|
|
||||||
|
# go valid namespace
|
||||||
|
go_env = {'BP':rpy.r.BPPARENTS, 'MF':rpy.r.MFPARENTS, 'CC': rpy.r.CCPARENTS}
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -0,0 +1,141 @@
|
||||||
|
import sys
|
||||||
|
import rpy
|
||||||
|
from pylab import gca, figure, subplot
|
||||||
|
from scipy import *
|
||||||
|
from lpls import *
|
||||||
|
import rpy_go
|
||||||
|
sys.path.append("../../fluents") # home of dataset
|
||||||
|
sys.path.append("../../fluents/lib") # home of cx_stats
|
||||||
|
import dataset
|
||||||
|
import cx_stats
|
||||||
|
from plots_lpls import plot_corrloads
|
||||||
|
|
||||||
|
######## DATA ##########
|
||||||
|
# full smoker data
|
||||||
|
DX = dataset.read_ftsv(open("../../data/smokers-full/Xfull.ftsv"))
|
||||||
|
DY = dataset.read_ftsv(open("../../data/smokers-full/Yg.ftsv"))
|
||||||
|
Y = DY.asarray()
|
||||||
|
# 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", 100)
|
||||||
|
# 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
|
||||||
|
co = 0.001
|
||||||
|
index = where(qvals<0.01)[0]
|
||||||
|
|
||||||
|
# Subset data
|
||||||
|
X = DX.asarray()
|
||||||
|
Xr = X[:,index]
|
||||||
|
gene_ids = DX.get_identifiers('gene_ids', index)
|
||||||
|
|
||||||
|
### Build GO data ####
|
||||||
|
|
||||||
|
print "Go terms ..."
|
||||||
|
goterms = rpy_go.goterms_from_gene(gene_ids)
|
||||||
|
terms = set()
|
||||||
|
for t in goterms.values():
|
||||||
|
terms.update(t)
|
||||||
|
terms = list(terms)
|
||||||
|
rpy.r.library("GOSim")
|
||||||
|
# Go-term similarity matrix
|
||||||
|
methods = ("JiangConrath","Resnik","Lin","CoutoEnriched","CoutoJiangConrath","CoutoResnik","CoutoLin")
|
||||||
|
meth = methods[0]
|
||||||
|
print "Term-term similarity matrix (method = %s)" %meth
|
||||||
|
if meth=="CoutoEnriched":
|
||||||
|
rpy.r('setEnrichmentFactors(alpha=0.1,beta=0.5)')
|
||||||
|
tmat = rpy.r.getTermSim(terms, verbose=False, method=meth)
|
||||||
|
# check if all terms where found
|
||||||
|
nanindex = where(isnan(tmat[:,0]))[0]
|
||||||
|
keep=[]
|
||||||
|
has_miss = False
|
||||||
|
if len(nanindex)>0:
|
||||||
|
has_miss = True
|
||||||
|
print "Some terms missing in similarity matrix"
|
||||||
|
keep = where(isnan(tmat[:,0])!=True)[0]
|
||||||
|
print "Number of nans: %d" %len(nanindex)
|
||||||
|
tmat_new = tmat[:,keep][keep,:]
|
||||||
|
new_terms = [i for ind,i in enumerate(terms) if ind in keep]
|
||||||
|
bad_terms = [i for ind,i in enumerate(terms) if ind not in keep]
|
||||||
|
# update go-term dict
|
||||||
|
for gene,trm in goterms.items():
|
||||||
|
for t in trm:
|
||||||
|
if t in bad_terms:
|
||||||
|
trm.remove(t)
|
||||||
|
if len(trm)==0:
|
||||||
|
print "Removing gene: %s" %gene
|
||||||
|
goterms[gene]=trm
|
||||||
|
terms = new_terms
|
||||||
|
tmat = tmat_new
|
||||||
|
# Z-matrix
|
||||||
|
# func (min, max, median, mean, etc),
|
||||||
|
# func decides on the representation of gene-> goterm when multiple
|
||||||
|
# goterms exist for one gene
|
||||||
|
Z, newind = rpy_go.genego_matrix(goterms, tmat, gene_ids, terms,func=mean)
|
||||||
|
Z = Z.T
|
||||||
|
# update X matrix (no go-terms available)
|
||||||
|
Xr = Xr[:,newind]
|
||||||
|
gene_ids = asarray(gene_ids)[newind]
|
||||||
|
|
||||||
|
|
||||||
|
######## LPLSR ########
|
||||||
|
print "LPLSR ..."
|
||||||
|
a_max = 5
|
||||||
|
aopt = 2
|
||||||
|
alpha=.5
|
||||||
|
T, W, P, Q, U, L, K, B, b0, evx, evy, evz = nipals_lpls(Xr,Y,Z, a_max, alpha)
|
||||||
|
|
||||||
|
# Correlation loadings
|
||||||
|
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)
|
||||||
|
|
||||||
|
# Significance Hotellings T
|
||||||
|
Wx, Wz, Wy, = jk_lpls(Xr, Y, Z, aopt)
|
||||||
|
tsqx = cx_stats.hotelling(Wx,W[:,:aopt])
|
||||||
|
tsqz = cx_stats.hotelling(Wz,L[:,:aopt])
|
||||||
|
|
||||||
|
|
||||||
|
## plots ##
|
||||||
|
figure(1) #rmsep
|
||||||
|
#bar()
|
||||||
|
figure(2) # Hypoid correlations
|
||||||
|
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)
|
||||||
|
|
||||||
|
figure(3)
|
||||||
|
subplot(221)
|
||||||
|
ax = gca()
|
||||||
|
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)
|
||||||
|
#title('Y correlation')
|
||||||
|
subplot(223)
|
||||||
|
ax = gca()
|
||||||
|
plot_corrloads(Rz, pc1=0, pc2=1, s=tsqz/10.0, c='r', zorder=5, expvar=evz, ax=ax)
|
||||||
|
#title('Z correlation')
|
||||||
|
subplot(224)
|
||||||
|
plot(arange(len(evx)), evx, 'b', label='X', linewidth=2)
|
||||||
|
plot(evy, 'g', label='Y', linewidth=2)
|
||||||
|
plot(evz, 'r', label='Z', linewidth=2)
|
||||||
|
legend(loc=2)
|
||||||
|
ylabel('Explained variance')
|
||||||
|
xlabel('Component')
|
||||||
|
show()
|
||||||
|
|
@ -0,0 +1,44 @@
|
||||||
|
|
||||||
|
def smdb_annot(orflist=None, input_fname='registry.genenames.tab', output_fname='yeast.annot'):
|
||||||
|
|
||||||
|
"""Reads registry.genenames.tab from the Stanford yeast
|
||||||
|
microarray database.
|
||||||
|
|
||||||
|
Available from:
|
||||||
|
ftp://genome-ftp.stanford.edu/pub/yeast/data_download/gene_registry/registry.genenames.tab
|
||||||
|
|
||||||
|
input: orf -- list of orfs (open reading frames)
|
||||||
|
file -- (optional) file to fetch info from
|
||||||
|
|
||||||
|
registry.genames contains:
|
||||||
|
|
||||||
|
0 = Locus name
|
||||||
|
1 = Other name
|
||||||
|
2 = Description
|
||||||
|
3 = Gene product
|
||||||
|
4 = Phenotype
|
||||||
|
5 = ORF name
|
||||||
|
6 = SGDID
|
||||||
|
|
||||||
|
NB! Other name, Gene product and Phenotype may have more
|
||||||
|
than one mapping. These are separated by |
|
||||||
|
|
||||||
|
Output: writes an annotation file
|
||||||
|
|
||||||
|
"""
|
||||||
|
outfile = open(output_fname, 'w')
|
||||||
|
header = "Orf\tLocus_id\tOther_name\tDescription\tGene_product\tPhenotype\tSGD_ID\n"
|
||||||
|
outfile.write(header)
|
||||||
|
text = open(input_fname, 'r').read().splitlines()
|
||||||
|
for line in text:
|
||||||
|
els = line.split('\t')
|
||||||
|
orf_name = els.pop(5)
|
||||||
|
if orf_name!='': # we dont care about non-named orfs
|
||||||
|
if orflist and orf_name not in orflist:
|
||||||
|
break
|
||||||
|
for e in els:
|
||||||
|
if e !='':
|
||||||
|
outfile.write(str(e) + "\t")
|
||||||
|
else:
|
||||||
|
outfile.write("NA")
|
||||||
|
f.write("\n")
|
||||||
Reference in New Issue