Comment mettre en œuvre la R du p.ajuster en Python

J'ai une liste de valeurs de p et je voudrais calculer l'ajuster les valeurs de p pour les comparaisons multiples pour le RAD. Dans R, je peux utiliser:

pval <- read.csv("my_file.txt",header=F,sep="\t")
pval <- pval[,1]
FDR <- p.adjust(pval, method= "BH")
print(length(pval[FDR<0.1]))
write.table(cbind(pval, FDR),"pval_FDR.txt",row.names=F,sep="\t",quote=F )

Comment puis-je mettre en œuvre ce code en Python? Voici mon feable tentative en Python, avec l'aide de Google:

pvalue_list [2.26717873145e-10, 1.36209234286e-11 , 0.684342083821...] # my pvalues
pvalue_lst = [v.r['p.value'] for v in pvalue_list]
p_adjust = R.r['p.adjust'](R.FloatVector(pvalue_lst),method='BH')
for v in p_adjust:
    print v

Le code ci-dessus déclenche une AttributeError: 'float' object has no attribute 'r' erreur. Quelqu'un peut-il aider à mon problème? Merci d'avance pour l'aide!

InformationsquelleAutor drbunsen | 2011-09-16
  1. 16

    Si vous souhaitez être sûr de ce que vous obtenez à partir de R, vous pouvez aussi indiquer que vous souhaitez utiliser la fonction dans le package R 'stats':

    from rpy2.robjects.packages import importr
from rpy2.robjects.vectors import FloatVector

stats = importr('stats')

p_adjust = stats.p_adjust(FloatVector(pvalue_list), method = 'BH')
    • Merci pour l'aide! Lorsque j'exécute le code, Python renvoie une ImportError: No module named packages erreur. Aucune idée de quel est le problème? Je suis en cours d'exécution de la R 2.13.1.
    • Je dirais que vous utilisez une version obsolète de rpy2. Essayez rpy2.__version__ en cas de doute. Le courant est 2.2.2.
    • Yep, cela fonctionne pour moi avec R 2.2. x. Malheureusement, je suis bloqué avec l'utilisation de R 2.13.1 sur un serveur distant. Toutes les suggestions?
    • hmmm... je me réfère à rpy2 version, pas de R versions. Demander une mise à niveau de rpy2 à vos administrateurs système ou de mise à niveau pour vous-même (pensez à utiliser le paquet Python 'virtualenv" pour créer vos propres Python).
    • Désolé pour la confusion. J'ai mal lu vos commentaires. J'ai mis à jour mon local rpy2 à 2,2 x et votre code a fonctionné. Merci beaucoup pour l'aide!
    InformationsquelleAutor lgautier
  2. 15

    Cette question est un peu vieux, mais il y a de multiples comparaison des corrections disponibles dans statsmodels pour Python. Nous avons

    http://statsmodels.sourceforge.net/devel/generated/statsmodels.sandbox.stats.multicomp.multipletests.html#statsmodels.sandbox.stats.multicomp.multipletests

    • Merci à vous, continuez le bon travail avec statsmodels!
    • Salut, une petite question sur le multipletests? Comment cette fonction de prendre soin de NaN valeurs dans la liste des valeurs de p lors de l'utilisation de BH? Il semble qu'il suppose que toutes les valeurs de p sont finis, c'est bien cela?
    InformationsquelleAutor jseabold
  3. 9

    Voici la fonction que j'utilise:

    def correct_pvalues_for_multiple_testing(pvalues, correction_type = "Benjamini-Hochberg"):                
"""                                                                                                   
consistent with R - print correct_pvalues_for_multiple_testing([0.0, 0.01, 0.029, 0.03, 0.031, 0.05, 0.069, 0.07, 0.071, 0.09, 0.1]) 
"""
from numpy import array, empty                                                                        
pvalues = array(pvalues) 
n = float(pvalues.shape[0])                                                                           
new_pvalues = empty(n)
if correction_type == "Bonferroni":                                                                   
new_pvalues = n * pvalues
elif correction_type == "Bonferroni-Holm":                                                            
values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]                                      
values.sort()
for rank, vals in enumerate(values):                                                              
pvalue, i = vals
new_pvalues[i] = (n-rank) * pvalue                                                            
elif correction_type == "Benjamini-Hochberg":                                                         
values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]                                      
values.sort()
values.reverse()                                                                                  
new_values = []
for i, vals in enumerate(values):                                                                 
rank = n - i
pvalue, index = vals                                                                          
new_values.append((n/rank) * pvalue)                                                          
for i in xrange(0, int(n)-1):  
if new_values[i] < new_values[i+1]:                                                           
new_values[i+1] = new_values[i]                                                           
for i, vals in enumerate(values):
pvalue, index = vals
new_pvalues[index] = new_values[i]                                                                                                                  
return new_pvalues
    • Excellente solution. Je l'ai porté sous python 3 et l'a placé sur un référentiel sur github. Si vous souhaitez m'ajouter votre nom à l'avis de droit d'auteur, merci de me fournir par MP.
    InformationsquelleAutor emre
  4. 5

    À l'aide de Python numpy bibliothèque, sans appel à R à tous, voici un raisonnablement efficace de la mise en œuvre de la BH méthode:

    import numpy as np
def p_adjust_bh(p):
"""Benjamini-Hochberg p-value correction for multiple hypothesis testing."""
p = np.asfarray(p)
by_descend = p.argsort()[::-1]
by_orig = by_descend.argsort()
steps = float(len(p)) / np.arange(len(p), 0, -1)
q = np.minimum(1, np.minimum.accumulate(steps * p[by_descend]))
return q[by_orig]

    (Basé sur le R code BondedDust posté)

    • Devrait être float(len(p)), sinon ça sera la division entière
    InformationsquelleAutor Eric Talevich
  5. 2

    (Je sais que ce n'est pas la réponse... juste essayer d'être utile.) Le BH code dans R p.ajustez est juste:

    BH = {
i <- lp:1L   # lp is the number of p-values
o <- order(p, decreasing = TRUE) # "o" will reverse sort the p-values
ro <- order(o)
pmin(1, cummin(n/i * p[o]))[ro]  # n is also the number of p-values
}
    InformationsquelleAutor 42-
  6. 1

    Vieille question, mais voici une traduction de la R RAD code en python (qui est probablement assez inefficaces):

    def FDR(x):
"""
Assumes a list or numpy array x which contains p-values for multiple tests
Copied from p.adjust function from R  
"""
o = [i[0] for i in sorted(enumerate(x), key=lambda v:v[1],reverse=True)]
ro = [i[0] for i in sorted(enumerate(o), key=lambda v:v[1])]
q = sum([1.0/i for i in xrange(1,len(x)+1)])
l = [q*len(x)/i*x[j] for i,j in zip(reversed(xrange(1,len(x)+1)),o)]
l = [l[k] if l[k] < 1.0 else 1.0 for k in ro]
return l
    InformationsquelleAutor Chrismit
  7. 0

    Bien, pour obtenir votre code de travail, je suppose que quelque chose de ce genre:

    import rpy2.robjects as R
pvalue_list = [2.26717873145e-10, 1.36209234286e-11 , 0.684342083821...] # my pvalues
p_adjust = R['p.adjust'](R.FloatVector(pvalue_list),method='BH')
for v in p_adjust:
print v

    Si p est.ajustez est assez simple, vous pouvez l'écrire en Python, de manière à éviter la nécessité de remettre en R. Et si vous voulez l'utiliser beaucoup, vous pouvez faire un simple wrapper Python:

    def adjust_pvalues(pvalues, method='BH'):
return R['p.adjust'](R.FloatVector(pvalues), method=method)
    InformationsquelleAutor Thomas K