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:
- Se ordena el dataframe en base a la columna
ThickStart
- 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)
- Se usa esa columna como índice.
- 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:
- En el dataframe
aux
iliar hacemos que el índice sea de nuevo un entero, para lo que hay que convertir los datetimes a timestamps.
- Creamos el dataframe
result
copiando el dataframe df
original, pero haciendo que el index sea la columna ThickStart (no hace falta ordenarla)
- 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.
- En
result
quitamos el índice para que vuelva a ser el original, 0, 1, 2... y así la columna ThickStart reaparezca
- 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).