1

Dada una tabla como esta en pandas: Tabla pandas

Necesito crear una columna que sea una media de la columna PercentMeth, pero que solo coja los valores dentro del rango de ThickStart +/- 50. Es decir si yo tengo un ThickStart de 1000191, que haga la media solo de los valores PercentMeth dados desde ThickStart 1000141 hasta 1000241.Los valores de ThickStart no son una serie seguida, es decir no es 1,2,3,4 si no que puede ser 1,2,4,9.

Yo he probado algo así pero me devuelve el valor PercentMeth de cada fila.

introducir la descripción de la imagen aquí

Adjunto un archivo ejemplo de datos de la web GenomeBrowser.

http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/wgEncodeHaibMethylRrbsA549Dm002p7dHaibSitesRep1.bed.gz

Y parte del código que tengo hasta ahora.

import pandas as pd

df = pd.read_csv('prueba2.bed', sep='\t', comment='t', header=None)
header = ['Chrom', 'ChromStart', 'ChromEnd', 'Name', 'Score', 'Strand', 'ThickStart', 'ThickEnd', 'ItemRgb', 'ReadCount', 'PercentMeth', 'BlockStarts']

df.columns = header[:len(df.columns)]

df.head()

def per_meth(fila):
    
    s = fila["ThickStart"] - 50
    e = fila["ThickStart"] + 50
    a = fila["PercentMeth"]
    b = fila["ThickStart"]
    c = 1
    
    for t in range(s, e):
        if b == t:
            a = a + ["PercentMeth"]
            c = c+1
            
        return (a/c)
´´´

Gracias.


5
  • 1
    No podemos hacer pruebas para dar con la solución si no proporcionas unos datos de ejemplo, ya sea un enlace a un archivo .csv con los datos, o bien un extracto de tu dataframe en la propia pregunta (pero no una imagen, por favor, no obligues a quien te ayude a copiar a mano los números). El código como imagen tampoco es muy apropiado, mejor si lo insertas como texto.
    – abulafia
    el 9 feb. 2022 a las 9:44
  • Agregando a lo de Abulafia, recuerda que en realidad no es nuestra obligación ayudarte ni nos pagan por ello, todo se basa en el deseo de ayudar, la pasion por programar y aumentar entre todos nuestros conocimientos. Si pones imagenes de código o dataframes en lugar de texto formateado, algunos miembros de la comunidad que sepan la respuesta no querrán ayudarte por tener que transcribir texto de tus imagenes... Otros motivos aquí: es.meta.stackoverflow.com/a/3976/158880
    – Dante S.
    el 9 feb. 2022 a las 11:40
  • Una cosa no me queda clara. Imagina que tu columna thickstart tiene los valores consecutivos desde 1000 hasta 2000 (el df tiene entoncs 1000 filas). ¿Cuántas filas querrías que tuviera el resultado? ¿Tambien 1000 o menos? Por ejemplo, la primera fila tendría la media de todos los thickstart comprendidos entre 950 y 1050, por tanto de las primeras 50 filas de df. Pero la siguiente fila tenía thikcstart 1001, entonces ¿debería contener la media de las 51 primeras filas del df otra vez? ¿O por haber formado parte ya del primer resultado se omiten las 50 primeras filas?
    – abulafia
    el 10 feb. 2022 a las 11:08
  • Dicho de otra forma.... la media esa proviene de "agrupar" (o resamplear) el thikcstart de 100 en 100 (50 por detrás y 50 por delante), o es una "rolling window" de 100?
    – abulafia
    el 10 feb. 2022 a las 11:08
  • Muy buenas, la siguiente se tendría que volver a calcular, es decir la primera calcularía la muestra 1000 desde 950 hasta 1050 y la siguiente la 1001 de 951 a 1051. Sería un rolling window
    – user269662
    el 10 feb. 2022 a las 12:14

1 respuesta 1

0

El problema que planteas es bastante complejo.

Tu enfoque inicial de usar apply() podría en principio ser válido, pero la función que aplicas estaba mal. La función debería recorrer el dataframe completo, y no una sola fila, buscando en él todas las filas que tengan en ThickStart valores en el rango especificado (el rango sí dependerá de qué fila estés procesando), para calcular el promedio de esos valores y retornarlo.

Es decir, sería así:

def per_meth(fila):
    b = fila["ThickStart"]
    mean = df[df.ThickStart.between(b-50, b+50)].PercentMeth.mean()
    return mean

Observa cómo, aunque la función recibe una fila, ésta sólo la usa para extraer el valor de ThickStart (que queda guardado en b), pero para calcular la media usa en cambio df, que no recibe como parámetro sino que es la variable global que contiene el dataframe completo. Usando .between() seleccionamos sólo aquellas filas cuyo ThickStart está dentro del rango buscado y hallamos la media.

Aunque lo anterior funciona para dataframes no muy grandes es tremendamente ineficiente y esto se debe a que para cada fila del dataframe se llamará a esa función, y en cada llamada se tiene que volver a procesar el dataframe completo. La complejidad algorítmica de esto crece con el cuadrado de N (siendo N el tamaño del dataframe). Para un dataframe de 1000 entradas esto termina rápido, pero si en vez de 1000 son 10000 (diez veces más) tardará cien veces más.

Tu dataframe es enorme (más de un millón de filas, concretamente 1464031) y este mecanismo por tanto requeriría un tiempo prohibitivo. Solo por curiosidad, intenté ejecutarlo pero me aburrí de esperar y lo aborté cuando llevaba más de una hora (y vete a saber si quizás habría tardado días).

Otro enfoque más eficiente

Hay que evitar recorrer el dataframe en cada iteración. Hay que evitar la complejidad cuadrática.

Se me ocurrió, ya que el problema es al final una rolling window, hacer uso del método .rolling() de pandas, pero el problema es que este método tiene limitaciones:

  • Se aplica sólo sobre el índice del dataframe (y no sobre una columna arbitraria como ThickStart)
  • El índice ha de estar ordenado crecientemente (y esta columna no lo está)
  • El índice ha de ser de tipo datetime

No obstante, se puede crear un dataframe auxiliar que cumpla estos objetivos, para lo cual:

  1. Se ordena el dataframe en base a la columna ThickStart
  2. Se trata la columna ThickStart como si fuese una estampa de tiempos, por ejemplo un contador de segundos. Eso permite convertir todos sus valores al tipo datetime. Las fechas y horas que salen son irrelevantes, pero gracias a esto podremos hacer una "rolling window" de 100s centrada (y así tendrá 50 por detrás y 50 por delante)
  3. Se usa esa columna como índice.
  4. Se aplica la ventana rodante, y se usa para calcular la media que querías, asignándolo a otra columna del dataframe auxiliar.

En código:

# 1. Ordenar ascendentemente por ThickStart
aux = df.sort_values(by="ThickStart")

# 2. Crear otra columna auxilar con el ThickStart convertido a datetime
aux["aux"] = pd.to_datetime(aux.ThickStart, unit="s")

# 3. Usar esta columna auxiliar como índice
aux = aux.set_index("aux")

# 4. Calcular la media de la ventana rodante
aux["Meth(+/- 50)"] = aux.rolling("100s", center=True).PercentMeth.mean()

Ya casi lo tenemos. Pero el resultado de esto no es tu dataframe original, pues en aux las filas están ordenadas por ThickStart, de menor a mayor. Este sería el aspecto del resultado en aux:

                     Chrom  ChromStart  Score  ...   ThickEnd  PercentMeth  Meth(+/- 50)
aux                                            ...                                      
1970-01-01 00:51:19   chrM        3079     71  ...       3080            6      4.500000
1970-01-01 00:51:34   chrM        3094     71  ...       3095            3      4.500000
1970-01-01 01:55:18  chr17        6918      7  ...       6919           86     97.250000
1970-01-01 01:55:26  chr17        6926      7  ...       6927          100     97.555556
1970-01-01 01:55:28  chr17        6928      7  ...       6929          100     97.555556
...                    ...         ...    ...  ...        ...          ...           ...
1977-11-24 17:05:44   chr1   249239144     24  ...  249239145           88     97.500000
1977-11-24 17:05:46   chr1   249239146     24  ...  249239147          100     97.500000
1977-11-24 17:05:56   chr1   249239156     24  ...  249239157           96     97.500000
1977-11-24 17:06:02   chr1   249239162     24  ...  249239163          100     97.500000
1977-11-24 17:06:05   chr1   249239165     24  ...  249239166          100     97.500000

El índice son fechas y horas que como dije no tienen relevancia pero han servido para encontrar la media que se buscaba, que está en la última columna.

Sólo nos queda asignar los resultados que acabamos de encontrar al dataframe original. El problema es que las filas no aparecen en el mismo orden, así que toca hacer alguna trampa adicional para lograrlo. En concreto:

  1. En el dataframe auxiliar hacemos que el índice sea de nuevo un entero, para lo que hay que convertir los datetimes a timestamps.
  2. Creamos el dataframe result copiando el dataframe df original, pero haciendo que el index sea la columna ThickStart (no hace falta ordenarla)
  3. Ahora tanto result como aux usan un índice común (aunque no estén en el mismo orden eso no importa) gracias a lo cual podemos asignar una nueva columna en result que sea una copia de la columna de resultados en aux. La asignación se hará en base al index y no en base al lugar que ocupan las filas.
  4. En result quitamos el índice para que vuelva a ser el original, 0, 1, 2... y así la columna ThickStart reaparezca
  5. Como consecuencia del paso 4 las columnas no están en el mismo orden que al inicio, así que las reordenamos.

En código:

import numpy as np

# 1. Hacer que el índice de aux sean enteros (la división por un millón
#    es necesaria porque los timestamps vienen en nanosegundos)
aux.index = aux.index.values.astype(np.int64)//1_000_000_000

# 2. Creación del dataframe de resultados que use como índice Thickstart
result = df.set_index("ThickStart")

# 3. Asignación del resultado
result["Meth(+/- 50)"] = aux.drop_duplicates("ThickStart")["Meth(+/- 50)"]

# 4. Eliminación del índice
result = result.reset_index()

# 5. Reorganización de las columnas
result = result[header[:-1] + [ "Meth(+/- 50)"]]

Y esto es lo que sale en result:

        Chrom  ChromStart  Score Strand  ...  ThickStart   ThickEnd  PercentMeth  Meth(+/- 50)
0        chr1   100002989      3      +  ...   100002989  100002990           33     33.000000
1        chr1   100003153     19      -  ...   100003153  100003154            5      5.000000
2        chr1     1000170     26      +  ...     1000170    1000171           27     27.111111
3        chr1     1000190     26      +  ...     1000190    1000191           54     33.100000
4        chr1     1000191     22      -  ...     1000191    1000192           41     39.181818
...       ...         ...    ...    ...  ...         ...        ...          ...           ...
1464026  chrY     9930940     39      -  ...     9930940    9930941           79     65.000000
1464027  chrY     9982886    352      +  ...     9982886    9982887           35      3.055556
1464028  chrY     9987654    154      +  ...     9987654    9987655            5      5.000000
1464029  chrY     9987800     46      +  ...     9987800    9987801            4      2.500000
1464030  chrY     9987801    100      -  ...     9987801    9987802            1      2.500000

Nota

No puedo garantizar que la solución propuesta (que por cierto se ejecuta en un par de segundos) salga lo mismo que la solución inicial (que tardaba horas o días), pues no las puedo comparar, al no tener la segunda.

Hice una prueba sobre un dataframe mucho más pequeño (1000 filas) y aunque la mayoría de valores salían idénticos, algunos de ellos no lo eran.

Tras investigar un poco encontré que ello se debe a que la ventana de 100 no coincide exactamente en tamaño con lo que produce .between(), ya que éste último incluye ambos extremos, y por tanto equivaldría a una ventana de 101 valores y no de 100. Cambiando en el código anterior donde pone "100s" por "101s" ya salen idénticos. O también añadiendo un parámetro a .between() para que excluya el extremo derecho.

No obstante hay que tener claro qué quieres hacer exactamente. En tu código original usabas range(s,e) y la función range() dejaría fuera el extremo derecho, por lo que en realidad sí eran 100 valores (50 por detrás y 49 por delante). Si quieres absoluta precisión en el resultado habría que aclarar este detalle acerca de si -50/+50 incluye o no a los extremos (y por tanto la ventana es de 101 y no de 100).

1
  • Parece que funciona perfectamente, probare con varios archivos para acabar de comprobarlo. Agradecerte la respuesta, la verdad aparte de resolver esta muy bien explicado. Muchas gracias.
    – user269662
    el 14 feb. 2022 a las 8:20

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.