r_in_python.py

from subprocess import *

def phyper(q, m, n, k):
    # Get an array of values to run
    rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
    runMe = []
    for i in range(len(q)):
        runMe.append('phyper('+str(q[i])+','+str(m[i])+','+str(n[i])+','+str(k[i])+',lower.tail=F)')
    runMe = '\n'.join(runMe)+'\n'
    out = rProc.communicate(runMe)
    return [line.strip().split(' ')[1] for line in out[0].strip().split('\n') if line]
    
# Get first principle component
def firstPrincipalComponent(matrix):
    ""
    Cacluate the first prinicipal component of a gene expression matrix.
    Input: Expression matrix gene (rows) x conditions (columns), expects that the
        python equivalent will be an array of arrays of expression values (rows).
    Returns: Array of first prinicipal component and variance explained by first
        principal component.
    ""
    # Fire up R
    rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
    runMe = []
    # Make the data into an R matrix
    rbind = 'm1 = rbind('
    rbind += ','.join(['c('+','.join([str(i) for i in row])+')' for row in matrix])
    rbind += ')'
    runMe.append(rbind)
    # Compute first principal component
    runMe.append('tmp.pc = try(princomp(t(m1)),TRUE)')
    runMe.append('pc.1 = NA')
    runMe.append('var.exp = NA')
    runMe.append('if(!class(tmp.pc)==\'try-error\') {')
    runMe.append('    pc.1 = tmp.pc$scores[,1]')
    runMe.append('    var.exp = ((tmp.pc$sdev^2)/sum(tmp.pc$sdev^2))[1]')
    runMe.append('}')
    runMe.append('var.exp')
    runMe.append('pc.1')
    runMe = '\n'.join(runMe)+'\n'
    out = rProc.communicate(runMe)
    # Process output
    splitUp = out[0].strip().split('\n')
    splitUp.pop(0) # Get rid of header dealie
    varExp = float(splitUp.pop(0).strip())
    pc1 = []
    for r1 in splitUp:
        pc1 += [float(i) for i in r1.split(' ') if i and (not i.count('[')==1)]
    # Return output
    return [pc1, varExp]
    
# Get a correlation p-value from R
def correlation(a1, a2):
    """
    Calculate the correlation coefficient and p-value between two variables.
    Input: Two arrays of float or integers.
    Returns: Corrleation coefficient and p-value.
    """
    # Fire up R
    rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
    runMe = []
    # Make the data into an R matrix
    runMe.append('c1 = cor.test(c('+','.join([str(i) for i in a1])+'),c('+','.join([str(i) for i in a2])+'))')
    runMe.append('c1$estimate')
    runMe.append('c1$p.value')
    runMe = '\n'.join(runMe)+'\n'
    out = rProc.communicate(runMe)
    # Process output
    splitUp = out[0].strip().split('\n')
    rho = float(splitUp[1])
    pValue = float((splitUp[2].split(' '))[1])
    return [rho, pValue]

# Compute survival p-value from R
def survival(survival, dead, pc1, age):
    """
    Calculate the survival correlation coefficient and p-value between two variables.
    Input: Four arrays of float or integers.
    Returns:
    """
    # Fire up R
    rProc = Popen('R --no-save --slave', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE)
    runMe = []
    # Make the data into an R matrix
    runMe.append('library(survival)')
    runMe.append('s1 = c('+','.join([str(i) for i in survival])+')')
    runMe.append('dead1 = c('+','.join(['\''+str(i)+'\'' for i in dead])+')')
    runMe.append('pc1 = c('+','.join([str(i) for i in pc1])+')')
    runMe.append('age1 = c('+','.join([str(i) for i in age])+')')
    runMe.append('scph1 = summary(coxph(Surv(s1,dead1==\'DEAD\') ~ pc1))')
    runMe.append('scph2 = summary(coxph(Surv(s1,dead1==\'DEAD\') ~ pc1 + age1))')
    runMe.append('scph1$coef[1,4]')
    runMe.append('scph1$coef[1,5]')
    runMe.append('scph2$coef[1,4]')
    runMe.append('scph2$coef[1,5]')
    runMe = '\n'.join(runMe)+'\n'
    out = rProc.communicate(runMe)
    # Process output
    splitUp = out[0].strip().split('\n')
    z1 = float((splitUp[0].split(' '))[1])
    pValue1 = float((splitUp[1].split(' '))[1])
    z2 = float((splitUp[2].split(' '))[1])
    pValue2 = float((splitUp[3].split(' '))[1])
    return [[z1, pValue1], [z2, pValue2]]

# To test the survival function
survival([10,20,30,10,15], ['DEAD','ALIVE','ALIVE','DEAD','ALIVE'], [0.1,0.25,0.4,0.6,0.9], [25,35,45,55,65])