Una propuesta novedosa para la estimación de la hora de la muerte a partir de datos post mortem de expresiones de genes

Tutores: Yolanda Larriba González, Miguel Alejandro Fernández Temprano Facultad de Ciencias Trabajo Fin de Grado Grado en Estadística Una propuesta novedosa para la estimación de la hora de la muerte a partir de datos post mortem de expresiones de genes Autor: Carlota María Prieto Martínez 1 Resumen El estudio del patrón de expresiones de los genes puede ser un indicador importante del estado de las células y de problemas de salud. Los genes circadianos, asociados al ciclo sueño-vigilia, presentan un patrón de expresión rítmico u oscilatorio, y regulan funciones biológicas básicas como la respiración o digestión. Alteraciones en los patrones de expresión de estos genes están relacionadas con patologías. Debido al riesgo o coste que supone la obtención de este tipo de datos, es habitual trabajar con datos de expresión post mortem para los que el momento en el que se tomaron las muestras es desconocido. Este trabajo propone una metodología novedosa, basada en inferencia con restricciones, para la estimación del orden temporal de expresiones de genes post mortem y su análisis a partir del ajuste de modelos paramétricos (Cosinor, FMM) y no paramétricos de señal oscilatoria. Por último, propone medidas de error y bondad de ajuste para comparar y validar los métodos utilizados. Se obtienen resultados de interés tanto desde el punto de vista metodológico (mejor comportamiento del nuevo orden temporal propuesto) como biológico (propuesta de posibles nuevos genes cíclicos). Abstract Gene expression pattern analysis is a marker of the cell states and diseases. Circadian genes, related to day-night cycle, display rhythmic or up-down-up patterns and govern basic biological functions. Pattern analysis is key to identify pathologies. However, in practice, gene expression data are not easy to collect since it may suppose a risk for health or could be expensive. Hence, post-mortem expression data, for which the moment at which samples were taken is unknown, are commonly used in practice. This work proposes a novel methodology, based or order restricted inference, to estimate the temporal order among postmortem gene expression data. Moreover, non-parametric and parametric models (Cosinor, FMM) for oscillatory signals are fitted to analyze these data. Finally, error and goodness of fit measures are used to compare and validate the methods employed. Interesting results are obtained from the methodological point of view (better behaviour of the new proposed temporal order) and from the biological (proposal of possible new cyclic genes). 2 Contenido 1. Introducción .................................................................................................... 3 1.1 Objetivos ................................................................................................... 6 1.2 Estructura del documento ......................................................................... 6 1.3 Asignaturas relacionadas .......................................................................... 7 1.4 Ampliación de materia ............................................................................... 7 2. Metodología .................................................................................................... 8 2.1 Métodos de estimación del orden temporal............................................... 9 2.1.1 Método naive: orden según el ToD ................................................... 10 2.1.2 Método ORI (Inferencia con Restricciones de Orden) ...................... 10 2.2 Modelos de señal oscilatoria ................................................................... 12 2.2.1 Modelo no paramétrico ..................................................................... 12 2.2.2 Modelos paramétricos ....................................................................... 12 2.3 Medidas de calidad de los modelos ........................................................ 16 3. Datos ............................................................................................................ 17 4. Resultados ................................................................................................... 20 4.1 Resultados globales ................................................................................ 21 4.2 Resultados para los genes core .............................................................. 22 4.3 Posibles nuevos genes rítmicos .............................................................. 40 5. Conclusiones ................................................................................................ 44 Referencias ...................................................................................................... 46 Anexo A: Código de generación de órdenes ................................................... 47 Anexo B: Ajuste de los modelos ....................................................................... 50 Anexo C: Funciones y librerías ......................................................................... 54 3 1. Introducción La bioestadística ha adquirido un lugar relevante en los últimos años gracias a los continuos avances en diversas áreas y campos biomédicos. Esta ciencia es una rama de la estadística que se ocupa de los problemas planteados dentro de las ciencias de la vida, como la biología o la medicina. Tan relevante es la bioestadística, que gran parte de la evidencia en salud está construida en base a ésta. (Castro, 2019) El crecimiento de los métodos cuantitativos en las ciencias biomédicas ha hecho de esta disciplina un elemento clave en muchas áreas entre las que están los ensayos clínicos y la cronobiología, área en la que se estudian los fenómenos fisiológicos cíclicos, o ritmos biológicos, en los seres vivos. Para este estudio es fundamental la medición de las expresiones de los genes en los diferentes tejidos de los organismos vivos. Los recientes avances en biotecnología están permitiendo realizar esta medición de la expresión de los genes de forma más segura, más sencilla y barata, produciendo un mayor rendimiento de la información (Shetty, 2020). La información obtenida con estas nuevas técnicas tiene gran potencial para el estudio de su expresión y de los sistemas regulatorios. Un paso clave en estos estudios es detectar genes involucrados en diversos procesos cuya expresión presenta patrones característicos. El análisis de este tipo de datos, desde el punto de vista estadístico supone un reto, no es trivial. (Colleen, A., & Doherty , 2010) Un gen es una unidad molecular que codifica un producto funcional específico, por ejemplo, una proteína. También son los responsables de transmitir información a la descendencia del organismo. Llamamos expresión de un gen al proceso mediante el cual la información codificada en un gen se utiliza para dirigir el montaje de una molécula de proteína. Este proceso está estrictamente regulado y permite que una célula responda a cambios de su entorno, controla la síntesis de proteínas y estas a su vez los distintos procesos biológicos. (J., y otros, 2004) El nivel de expresión de los genes es dinámico y los patrones que estos siguen son fundamentales para entender los procesos biológicos, desde la inflamación hasta el envejecimiento. El estudio de los patrones de las expresiones de los genes se empieza a ver muy útil para diagnosticar enfermedades como el cáncer. (Roth, 2002) 4 En concreto, en este trabajo nos centraremos en los genes cuya expresión sigue un ritmo circadiano, conocidos como genes circadianos. El ritmo circadiano es un ciclo natural de cambios físicos, mentales y de comportamiento que experimenta un individuo en un ciclo de 24 horas. Este ritmo permite a un organismo adaptar su fisiología con anticipación entre la noche y el día provocando oscilaciones en un conjunto diverso de procesos biológicos y así poder preparar al individuo para responder a condiciones ambientales predecibles (Zhang, Lahens, Ballance, Hughes, & Hogenesch, 2014); lo que se traduce en un patrón de expresión oscilatorio en el que hay un único máximo y un único mínimo en el ciclo de expresión, que se denominará patrón up-down- up, coincidiendo el máximo de la expresión con el momento del día (ciclo) en el que el gen realiza la acción biológica. Este tipo de señales que se conocen como señal circular y cuya definición formal se establecerá más adelante, subyacen en la gran mayoría de ritmos biológicos. La expresión de los genes es muy relevante puesto que es un indicador del estado de las células y se puede asociar con la aparición de enfermedades como la neurodegeneración, depresión, trastornos metabólicos, etc. (Liu, Gershon, & Kelsoe, 2017) Los ritmos circadianos anormales también pueden estar relacionados con la obesidad, la diabetes, la depresión, el trastorno bipolar, el trastorno afectivo estacional y los trastornos del sueño. (Zhang, Lahens, Ballance, Hughes, & Hogenesch, 2014) Por tanto, el análisis de los patrones de las expresiones de los genes circadianos puede suponer un gran avance para la detección de algunas enfermedades. Sin embargo, aun cuando los procesos para su obtención han mejorado sensiblemente en los últimos tiempos, la obtención de los mismos sigue siendo complicada y supone un riesgo para la salud. (Valadares, Gorki, Liebold , & Hoenicka, 2017). Por ello es habitual que el análisis de expresiones de genes se efectúe a partir de expresiones post mortem, es decir, de personas ya fallecidas. Pero en este orden de circunstancias es relevante tener en cuenta que la muerte clínica, paro cardíaco, no implica el paro de las funciones biológicas y por tanto de la expresión de los genes. Esto supone, en la práctica, que la estimación del momento de la muerte sea desconocida o provenga de estimaciones imprecisas. Como ejemplo de esta cuestión, en el panel izquierdo de la Figura 1 se muestra la representación, ordenado según el orden dado por el momento de fallecimiento estimado (ToD), de la expresión del gen NPRL2, gen típicamente circadiano, que, al contrario de lo esperado, presenta un patrón bastante alejado de ser rítmico. Por el contrario, en el panel derecho se muestra otra estimación del orden temporal que se ajusta de forma más adecuada al ritmo circadiano esperado para el gen NPRL2. Por lo tanto, parece claro que esa estimación hecha de forma directa a partir del momento estimado del fallecimiento puede ser mejorable y que se requiere de un mejor procedimiento de estimación de orden temporal, previo al análisis de expresiones circadianas a partir de modelos de señales oscilatorias. 5 Así pues, uno de los objetivos de este trabajo consiste en analizar y comparar distintos métodos de estimación del orden temporal y modelos para señales oscilatorias para el análisis de datos de expresión de genes post mortem. Hay que señalar también que, una vez ordenados los datos, existen diferentes modelos para ajustar una señal rítmica a los datos. Los modelos que se van a considerar en este trabajo, y que se describen con detalle en la sección de metodología, son el clásico modelo Cosinor (Cornelissen, 2014) y otros dos modelos más flexibles, un modelo no paramétrico específicamente diseñado en (Larriba , Rueda , Fernández , & Peddada , 2019) para el ajuste de una señal updown-up, y un modelo paramétrico, definido en (Rueda , Larriba, & Peddada , 2019), que permite, a diferencia de Cosinor, el ajuste de modelos asimétricos y cuyos parámetros son fácilmente interpretables. Estos modelos, juntamente con las correspondientes medidas de bondad de ajuste, permitirán valorar el desempeño de los órdenes propuestos en este trabajo a la hora de describir los datos analizados. Los datos que se van a considerar en este trabajo, y que se describirán posteriormente con más detalle, corresponden a 104 individuos sanos y 46 individuos esquizofrénicos y se han obtenido de CommonMind Consortium (CMC). Estos datos se han considerado también en (Seney , y otros, 2019) donde se estudian los posibles cambios en la ritmicidad de algunos genes debidos a la esquizofrenia. Dado que en ese trabajo se considera como punto de partida el orden dado por el momento estimado de fallecimiento, orden que se pretende mejorar, otro de los objetivos que se persigue en este TFG es el de confirmar, o no, las conclusiones obtenidas en ese otro estudio en lo que se refiere a los cambios de ritmicidad producidos por la enfermedad. Figura 1: Comparación de dos órdenes temporales para el mismo gen NPRL2. En el panel izquierdo la estimación ha sido realizada con ToD. En el panel derecho la estimación ha sido realizada con el método ORI que se propone en este trabajo. 6 1.1 Objetivos Los objetivos de este trabajo, que se han esbozado durante la descripción anterior, pueden enumerarse de la forma siguiente: • Estimación del orden temporal para conjuntos de datos de expresiones de genes en los que este orden es desconocido por tratarse de datos de expresiones post mortem, mediante metodologías alternativas a la estimación directa de la hora de muerte. • Ajuste de las expresiones de los genes en base a patrones rítmicos mediante modelos paramétricos y no paramétricos. • Valoración global de los órdenes estimados y los ajustes obtenidos en los modelos anteriores. • Análisis de los resultados obtenidos mediante la metodología anterior en un conjunto de datos recientemente utilizado en la literatura para estudiar la influencia de la esquizofrenia en los patrones rítmicos. Valoración de las hipótesis establecidas en la literatura a la luz de esta nueva metodología. • Búsqueda de posibles nuevos genes con patrones rítmicos no observados en metodologías previas. 1.2 Estructura del documento La memoria de este trabajo fin de grado se compone de los siguientes capítulos: Metodología: En este capítulo se desarrolla la metodología utilizada en este trabajo. Primeramente, se describen los métodos de estimación del orden temporal, a continuación, los modelos que se utilizan para los ajustes de los datos y finalmente las medidas utilizadas para la valoración de los ajustes. Datos: Se describen los conjuntos de datos utilizados, de dónde se han obtenido y los tratamientos previos que se han efectuado en los mismos. Resultados: En este capítulo se describen, comparan y valoran los resultados obtenidos de acuerdo a los diferentes órdenes y modelos empleados en el trabajo. Asimismo, se valoran las hipótesis establecidas en la literatura sobre la influencia de la esquizofrenia en el carácter rítmico de los genes y se ofrecen hipótesis sobre nuevos posibles genes rítmicos sugeridos por la metodología. Conclusiones: Se resumen los resultados obtenidos y se sugieren posibles líneas de desarrollo futuro a partir de dichos resultados. Anexos: Se incluye el código R utilizado para la elaboración del trabajo. 7 1.3 Asignaturas relacionadas La relación entre las técnicas utilizadas en este trabajo y las diferentes asignaturas del grado es la siguiente: • Minería de datos: en ella se estudia el tratamiento previo que se debe realizar a los datos. • Computación estadística: esta asignatura aporta una base sólida de R, lenguaje utilizado para la realización de este trabajo. • Modelos Lineales y Modelos Estadísticos Avanzados: en estas asignaturas se imparten las bases para la comprensión de modelos estadísticos. • Modelos de Investigación Operativa y Algoritmos y Computación: en estas asignaturas se estudia y se dan soluciones al TSP (Problema del viajante). 1.4 Ampliación de materia Para poder desarrollar este trabajo fin de grado ha sido necesario el estudio de los siguientes contenidos adicionales a lo estudiado en el grado. • Modelos de señal oscilatoria: estos modelos han sido utilizados para realizar el ajuste de las expresiones de los genes. • Método ORI: método propuesto para la estimación del orden temporal. 8 2. Metodología En este capítulo se presenta la definición de señal circular, los métodos de estimación de orden temporal, los modelos de señales oscilatorias ajustados y las métricas aplicadas en este trabajo. Sea 𝑥𝑗= (𝑥1𝑗,𝑥2𝑗,…,𝑥𝑛𝑗)’ con 𝑗= 1…𝑝 los datos de expresión del gen j para los instantes de tiempo 𝑖= 1…𝑛, siendo p el número de genes y n el número de instantes de tiempo. Sea X el conjunto de datos de los p genes. Se dice que 𝑥𝑗 sigue un modelo de señal circular si 𝑥𝑗=µ+ 𝜀𝑗 (1) donde µ es una señal circular tal y como se define a continuación. Definición 1: Definición de señal circular. (Caso U<L) Una señal µ es una señal circular si y solo si µ ∈ 𝐶= ⋃ 𝐶𝐿𝑈 𝐿,𝑈 (2) donde L =𝑎𝑟𝑔𝑚𝑖𝑛1…𝑛 𝜇𝑖 , U =𝑎𝑟𝑔𝑚𝑎𝑥1…𝑛 𝜇𝑖 , 𝐿,𝑈 ∈{1,…,𝑛} 𝑦 ∁𝐿𝑈={𝜇 𝜖 ℝ𝑛: 𝜇1 ≤⋯≤𝜇𝑈≥⋯≥ 𝜇𝐿≤⋯≤𝜇𝑛≤𝜇1 }. Definición 2: Definición de orden circular Se dice que una señal 𝝓 sigue un orden Circular en el espacio Circular si y solo si ϕ ∈ 𝐶0 = {𝜙 ∈ [0,2𝜋)𝑛 ∶ 𝜙1⪯ · · · ⪯ 𝜙𝑛 ⪯ 𝜙1} (3) Donde ≼ se lee como “es seguido por”. En este caso decimos que ϕ sigue un orden circular. La siguiente fórmula presenta la equivalencia entre una señal circular en el espacio euclídeo y una señal circular en el espacio circular: 15 Otros parámetros importantes del modelo son los picos y sus tiempos que se calculan: 𝑡𝑈= α + 2atan (1 ω tan(−β/2)) (11) 𝑡𝐿= α+ 2atan (1 ω tan(𝜋−β) 2)) (12) Y los valores de la señal en esos puntos son: 𝑍𝑈= 𝑀 + 𝐴 (13) 𝑍𝐿= 𝑀− 𝐴 (14) Figura 5: Influencia de los parámetros β ω en el modelo FMM con M=0, A=1 y α =0. Imagen obtenida de (Pérez, 2020) 16 Relación entre el modelo Cosinor y el modelo FMM El modelo Cosinor es un caso particular del modelo FMM. Esto sucede cuando 𝜔=1 donde 𝜑 = 𝛽 − 𝛼. Esto se puede ver de manera directa: ϕ(t)= β+ 2arctan(ωtan(t − α 2))= β + 2arctan(tan(t − α 2)) =t +β − α (15) 2.3 Medidas de calidad de los modelos Por último, se presentan las distintas medidas utilizadas para la validación y comparación de los distintos métodos de estimación del orden y modelos de señales oscilatorias ajustados usados en el trabajo. Medida de error Para los modelos como medida del error se ha calculado el MSE (Error Cuadrático Medio). MSEi=∑(Xi − Xi) 2 ni=1 n (16) 𝑀𝑆𝐸=∑𝑀𝑆𝐸𝑖 𝑃 𝑖=1𝑝 (17) Donde 𝑋𝑖  es el valor ajustado por el modelo para los datos de la expresión del gen 𝑖, 𝑋𝑖 los datos de la expresión del gen 𝑖. Medida de bondad de ajuste Y como medida de bondad de ajuste se ha calculado el R2. 𝑅2=1 − ∑(𝑋𝑖− 𝑋𝑖 )2 𝑛 𝑖=1 ∑(𝑋𝑖− 𝑋)2 𝑛 𝑖=1 (18) Donde 𝑋𝑖  es el valor ajustado por el modelo para los datos de la expresión del gen 𝑖, 𝑋𝑖 los datos de la expresión del gen 𝑖, 𝑋 es el valor medio. 17 3. Datos Al tratarse de pacientes, los datos están anonimizados, y cada individuo tiene asociado un identificador, de forma que, para conseguir relacionar un individuo con los valores de la expresión de sus genes, conocer la hora de su muerte (ToD) y saber si es control o esquizofrénico se deben cruzar tales identificadores en los ficheros correspondientes. Una vez efectuado este tratamiento se tienen dos conjuntos de datos, con la expresión de los genes y el ToD, uno para individuos controles y otro para individuos esquizofrénicos. El conjunto de datos de individuos control está formado por 104 observaciones de las cuales, como puede verse en la Tabla 1, el 78% (81) son hombres y el 22% (23) son mujeres. El 83% (86) de la totalidad de las observaciones pertenecen a personas de raza blanca, el 16% (17) a personas de raza negra y el 1% (1) a personas asiáticas. El 59% (61) de las observaciones vienen de PIT (Pittsburgh) y el 41% (43) de MSSM (Mt. Sinai School of Medicine). Estos individuos tienen una edad media de 48.4 años. Tabla 1: Distribución de los individuos control por sexo, raza, lugar y edad. Por otro lado, el conjunto de datos de pacientes esquizofrénicos está formado por 46 individuos. Como puede verse en la Tabla 2, el 70% (32) son hombres y el 30% (14) son mujeres. El 74% (34) son de raza blanca y el 26% (12) son de raza negra. El 48% (22) proviene de PIT y el 52% (22) de MSSM. Su edad media es de 50.1 años. Tabla 2: Distribución de los individuos esquizofrénicos por sexo, raza, lugar y edad Control Sexo Raza Lugar Edad Hombre Mujer Blanca Negra Asiática PIT MSSM 78% (81) 22% (23) 83% (86) 16% (17) 1% (1) 59% (61) 41% (43) 48.4 Esquizofrénico Sexo Raza Lugar Edad Hombre Mujer Blanca Negra Asiática PIT MSSM 70% (32) 30% (14) 74% (34) 26% (12) 0% (0) 48% (22) 52% (22) 50.1 18 Todos estos individuos han sido seleccionados en base a tres criterios recogidos en el trabajo (Seney , y otros, 2019): 1. Sujetos para los que la estimación del ToD es conocida y su deceso se desencadenó de forma repentina. 2. Sujetos menores de 65 años. 3. Sujetos con un intervalo post mortem inferior a las 30 horas. De todos estos individuos se recoge la expresión de 13914 genes en su hora de muerte. En la práctica se trabaja con dos matrices: de individuos sanos o control, 𝑀13914𝑥104 con 13914 filas de genes y 104 columnas de individuos y de individuos esquizofrénicos 𝑆𝐶𝑍13914𝑥46 con 13914 filas de genes y 46 columnas de individuos. Tratamiento de los datos Se precisa de dos rescalados sobre los datos para permitir la comparación simultánea de modelos y órdenes, así como para adecuarse a los requisitos de los modelos ajustados. El primero de ellos se realiza sobre los datos de las expresiones de los genes para poder compararlos, es un escalado al intervalo [−1,1]. La transformación utilizada para este escalado en el eje de las ordenadas es la siguiente: 𝑑𝑎𝑡𝑜𝐸𝑠𝑐𝑎𝑙𝑎𝑑𝑜= y − ymin ymax−ymin (19) En segundo lugar, se realiza un escalado de los tiempos ToD, es decir, en el eje de abscisas. Los datos originales vienen dados en el intervalo [−6,18], tiempos habituales en la escala Zeitgeber. Sin embargo, tanto el modelo Cosinor como el modelo FMM requieren que estos tiempos vengan dados en el intervalo [0,2𝜋]. Para este segundo escalado se ha aplicado la siguiente fórmula: 𝑑𝑎𝑡𝑜𝐸𝑠𝑐𝑎𝑙𝑎𝑑𝑜=2𝜋 x − xmin xmax−xmin=x+ 6 12 𝜋 (20) Nótese que el orden ORI que se ha descrito anteriormente no proporciona un instante de tiempo exacto para cada uno de los individuos, sino que solamente proporciona una ordenación entre los mismos. En consecuencia, para poder 19 estimar los modelos paramétricos, se establecen dos conjuntos equiespaciados de datos en el intervalo [0,2𝜋], teniendo en cuenta los tamaños muestrales de cada conjunto de datos. Concretamente, para los individuos control se toman 104 valores equiespaciados en el intervalo [0,2𝜋], mientras que para los esquizofrénicos se toman 46 en ese mismo intervalo. 20 4. Resultados Una vez obtenidos los datos de los individuos se comienza con su procesamiento y análisis. En esta sección mostraremos los principales resultados recogidos en la sección de objetivos del trabajo. Generación de los órdenes A continuación, se describe cómo se han obtenido las estimaciones de los órdenes ToD y ORI propuestos en el trabajo. Para calcular el orden ToD se ordenan los individuos según su hora de muerte de forma ascendente. A partir del método ORI generaremos dos órdenes. El primero de ellos se estimará a partir de un subconjunto de genes, que llamaremos genes core. Se dice que un gen es core si las proteínas que sintetiza son esenciales en la generación y regulación de ritmos circadianos. Los genes core, presentan patrones de expresión rítmicos. Para el segundo de los órdenes se emplean todos los genes de los datos. Para los individuos control los genes core utilizados para calcular el orden son: ARNTL, NPAS2, CLOCK, NFIL3, CRY1, NR1D1, BHLHE41, NR1D2, DBP, CIART, PER1, PER3, TEF, HLF, CRY2, PER2; Esta selección contiene genes core muy conocidos en la biología, está basado en (Wu, 2020). Sin embargo, para los individuos esquizofrénicos el subconjunto de genes utilizado es diferente ya que la biología circadiana en esta patología es distinta, en el sentido de que genes rítmicos en pacientes sanos pueden dejar de serlo en pacientes esquizofrénicos y genes que presentan ritmicidad en esquizofrénicos pueden dejar de presentarla en sanos. Para estos el conjunto de genes core seleccionado es: CIART, WNT10B, LAMB3, OPRL1, CYB561, HDAC8, NIM1K, DUBR, KRT17P1, EBP, PGBD2, CTSK, ZBTB22, NPRL2, IFT122, NFATC4, RNF112, VOPP1, NIT1, USF1. Esta selección de genes se ha hecho en base a los genes que presentan más ritmicidad para los individuos esquizofrénicos de acuerdo con (Seney , y otros, 2019) Finalmente, se calcula tanto para individuos control como para esquizofrénicos el orden ORI con todos los genes. Ajuste de los modelos Una vez obtenidos los órdenes se procede al análisis de los datos de expresión de los genes con el modelo no paramétrico (NP), el modelo Cosinor y el modelo FMM. 21 Para ello, debe tenerse en cuenta que el modelo Cosinor y FMM necesitan que los instantes temporales estén en el intervalo [0,2𝜋]. En el caso del orden ToD deberemos escalar las horas con la fórmula (20). Para los órdenes ORI se requiere que los vectores de instantes temporales mencionados al final de la sección anterior estén también en el intervalo [0,2𝜋]. 4.1 Resultados globales En esta sección se van a considerar los resultados globales obtenidos para cada uno de los tres órdenes y modelos considerados, teniendo en cuenta todos los genes disponibles, tanto para el grupo de control como para los pacientes de esquizofrenia. Para las valoraciones que se hacen en este apartado se ha considerado como medida de comparación el MSE puesto que permite una mejor valoración de la globalidad de los resultados. Las Tablas 3 y 4 contienen los valores globales de MSE para los conjuntos de datos de control y de pacientes esquizofrénicos respectivamente. En cada una de las tablas están los valores para los tres órdenes considerados (ORI, ORI Reducido y ToD) y los tres modelos que se han ajustado con esos órdenes (NP, FMM y Cosinor). Puede verse en las tablas que en todos los casos el modelo NP tiene un MSE menor que los otros dos modelos. Esto es consecuencia de que este modelo no paramétrico impone menos restricciones a los ajustes. Además, para el modelo FMM siempre se obtienen valores de MSE menores que para el modelo Cosinor como consecuencia de que el modelo Cosinor es un caso particular del modelo FMM. Es claro entonces que estos resultados eran esperables como consecuencia de la definición de cada uno de los modelos. Tabla 3 : MSE para los datos control para los tres órdenes valorados y los tres modelos considerados. Tabla 4: MSE para los datos de esquizofrénicos para los tres órdenes valorados y los tres modelos considerados. Control MSE ORI ORI Reducido ToD NP 0.0645 0.0827 0.1164 FMM 0.0937 0.1118 0.1343 Cosinor 0.1056 0.1217 0.1455 Esquizofrénicos MSE ORI ORI Reducido ToD NP 0.0636 0.0937 0.1247 FMM 0.1052 0.1329 0.1564 Cosinor 0.1251 0.1599 0.1864 22 Una segunda observación, de más interés que la anterior, es que el orden ToD es globalmente peor, en lo que se refiere al MSE, que los otros dos como puede observarse para ambos conjuntos de datos y para todos los modelos. Este resultado también era esperable puesto que los órdenes ORI se han construido con el objetivo de minimizar una función de error relacionada con el MSE del modelo NP. De todos hay que notar que, a diferencia de lo que ocurre con los modelos, no había seguridad total en este punto puesto que no se resuelve el problema directo de minimización sino uno relacionado con él mediante el TSP. También es interesante comentar las diferencias observadas en los resultados de los órdenes ORI y ORI Reducido. Se ve en las tablas como el MSE para el orden ORI es inferior al del ORI Reducido. Esto puede ser consecuencia de que el orden ORI está establecido a partir del conjunto completo de genes mientras que ORI Reducido está generado a partir de los genes core, cuyo carácter rítmico es bien conocido en la literatura. Esto quiere decir que el orden ORI no es necesariamente mejor que el ORI Reducido puesto que al utilizar todos los genes en la generación del orden se puede estar produciendo un efecto de “sobreajuste” y como consecuencia resultar que los genes core, con patrones de expresión claramente rítmicos, no aparezcan como tal en el orden ORI. Este es un posible efecto que se estudiará posteriormente para decidir entre ambos órdenes ORI. Por último, se observa que en general los resultados para los pacientes esquizofrénicos son, en términos de MSE, peores que para los controles. Esta observación parece estar en concordancia con las hipótesis establecidas en (Seney , y otros, 2019) que sugieren una pérdida de ritmicidad en los pacientes esquizofrénicos. 4.2 Resultados para los genes core En esta sección se describen los resultados obtenidos para los genes core de los dos conjuntos de datos utilizados (control y esquizofrénicos) para los distintos órdenes y modelos considerados en este trabajo. Para la comparación de los modelos se ha utilizado el R2 como medida de bondad de ajuste de cada uno de los modelos considerados ya que, de este modo, se actúa como, por ejemplo, en los modelos de regresión y las comparaciones son más intuitivas. La Tabla 5 contiene los resultados relativos a los genes core para el grupo de control. Los genes están ordenados en orden decreciente considerando el valor obtenido para R2 en el orden ORI Reducido para el modelo FMM. Puede observarse, en primer lugar, que los valores de R2 correspondientes al orden ORI Reducido son mejores (excepto para CRY2 y NPAS2) que los correspondientes al orden ORI. Esta observación confirma que el orden ORI Reducido ajusta mejor los genes core que, como ya se ha comentado, son 23 aquellos cuya ritmicidad es conocida en la literatura y nos permite afirmar que parece conveniente utilizar el orden ORI Reducido en el resto de las conclusiones en lugar del orden ORI generado utilizando todos los genes disponibles en el conjunto de datos. En cuanto al orden ToD como ya se comentó en los resultados globales, salvo un par de excepciones muy puntuales (genes PER2 y ARNTL para el modelo Cosinor), los resultados de los ajustes son peores que los obtenidos por el orden ORI Reducido. Genes Core Control R2 ORI ORI Reducido ToD Gen NP FMM Cosinor NP FMM Cosinor NP FMM Cosinor NR1D1 0,625 0,417 0,392 0,871 0,684 0,653 0,391 0,225 0,208 CIART 0,177 0,074 0,030 0,827 0,676 0,541 0,585 0,494 0,479 PER1 0,637 0,431 0,420 0,850 0,654 0,613 0,438 0,288 0,282 PER3 0,588 0,423 0,337 0,793 0,620 0,523 0,413 0,275 0,251 DBP 0,561 0,386 0,379 0,746 0,596 0,562 0,387 0,267 0,236 BHLHE41 0,620 0,448 0,253 0,759 0,552 0,530 0,301 0,122 0,054 NR1D2 0,549 0,332 0,231 0,687 0,454 0,350 0,333 0,178 0,167 CRY2 0,652 0,438 0,384 0,643 0,431 0,341 0,387 0,196 0,187 CRY1 0,402 0,223 0,152 0,560 0,389 0,355 0,420 0,276 0,225 NPAS2 0,748 0,537 0,462 0,534 0,361 0,132 0,204 0,069 0,010 PER2 0,521 0,338 0,307 0,546 0,331 0,169 0,413 0,262 0,216 HLF 0,541 0,325 0,239 0,548 0,328 0,309 0,199 0,086 0,020 TEF 0,382 0,204 0,184 0,506 0,327 0,303 0,396 0,217 0,174 NFIL3 0,434 0,241 0,069 0,609 0,303 0,195 0,390 0,252 0,181 CLOCK 0,756 0,522 0,396 0,540 0,290 0,267 0,205 0,114 0,018 ARNTL 0,295 0,137 0,007 0,476 0,289 0,178 0,417 0,242 0,236 Tabla 5: Valores de R2 para los ajustes de los genes core en el grupo de control para los diferentes órdenes y modelos considerados. A continuación, se presentan tres ejemplos de genes core del grupo de control para los que se ve con cierta claridad como el orden ORI Reducido proporciona ajustes mejores que el orden ToD utilizado en (Seney , y otros, 2019). El ajuste proporcionado por el modelo NP está en verde, el de Cosinor en rojo y el que arroja FMM en azul. 24 Es destacable mencionar que los tres genes que se incluyen como ejemplo, PER3, CRY2 y CLOCK, juegan un papel decisivo en la estructura molecular del reloj circadiano en mamíferos, ubicado en la base del núcleo supraquiasmático. A nivel molecular, las interacciones entre estos genes permiten "orquestar"/sincronizar y regular las distintas funciones biológicas que se suceden a lo largo del día en los distintos tejidos y órganos (Panda , y otros, 2002). Gen CLOCK Figura 6: Gen CLOCK con el orden ORI Reducido. Figura 7 : Gen CLOCK con el orden ToD. ORI Reducido Modelo R2 NP 0,540 FMM 0,290 Cosinor 0,267 ToD Modelo R2 NP 0,205 FMM 0,114 Cosinor 0,018 31 Figura 16: Parámetros  y  del modelo FMM para los genes core del grupo de control en el conjunto de control y en el de pacientes con esquizofrenia. Figura 17: Parámetros A del modelo FMM y valor de tU estimado con este mismo modelo para los genes core del grupo de control en el conjunto de control y en el de pacientes con esquizofrenia. 32 La Figura 17 muestra los valores del parámetro A (amplitud) del modelo FMM contra los valores estimados con este mismo modelo para 𝑡𝑈 (momento de máxima expresión del gen) para las mismas condiciones descritas para la Figura 16. Se observa como los valores de A son, en general, mayores o iguales para el conjunto de pacientes con esquizofrenia y que también cambian los valores estimados para 𝑡𝑈 siendo, en general, más bajos para el grupo de los pacientes control. La Figura 18 permite la comparación entre el grupo de control y el de pacientes con esquizofrenia, de los valores del parámetro  y de R2 para el modelo FMM para los genes core del conjunto de datos de control. Se observa como los valores de R2 son en general más bajos para los pacientes de esquizofrenia. Figura 18: Parámetro  y valor de R2 del modelo FMM para los genes core del grupo de control en el conjunto de control y en el de pacientes con esquizofrenia. 33 Finalmente, la Figura 19 permite la comparación entre el grupo de control y los pacientes esquizofrénicos de los valores de R2 del conjunto de los genes core del grupo de control. Se observa una disminución en general de este valor en el grupo de los pacientes esquizofrénicos ya que la mayor parte de los valores están por encima de la diagonal del gráfico. Esta disminución es esperable puesto que los genes que se están considerando son los genes core del grupo de control. Figura 19: Valores de R2 del modelo FMM para los genes core del grupo de control en el conjunto de control y en el de pacientes con esquizofrenia. A continuación, se ofrecen el mismo tipo de gráficos anteriores para los genes core del grupo de esquizofrenia con el objetivo de valorar los cambios que aparecen en estos genes entre los grupos de control y de pacientes esquizofrénicos. La Tabla 8 contiene la identificación de los genes core del grupo de pacientes esquizofrénicos utilizada en los gráficos que se ofrecen. 34 Equivalencia gen – letra Genes core esquizofrenia Gen Letra CIART a WNT10B b LAMB3 c OPRL1 d CYB561 e HDAC8 f NIM1K g DUBR h KRT17P1 i EBP j PGBD2 k CTSK l ZBTB22 m NPRL2 n IFT122 o NFATC4 p RNF112 q VOPP1 r NIT1 s USF1 t Tabla 8: Identificación de los genes core del grupo de esquizofrenia En la Figura 20 se consideran los parámetros  y  del modelo FMM para los genes core del conjunto de pacientes con esquizofrenia. En este caso no se observan tendencias claras, con los genes para el grupo de esquizofrenia ocupando en general la parte central del gráfico. 35 Figura 20: Parámetros  y  del modelo FMM para los genes core del grupo de esquizofrenia en el conjunto de control y en el de pacientes con esquizofrenia. La Figura 21 muestra los valores del parámetro A del modelo FMM contra los valores estimados con este mismo modelo para 𝑡𝑈 para las mismas condiciones descritas para la Figura 20. Se observa como los valores de A son, en general, más bajos para el conjunto control y que también cambian los valores estimados para 𝑡𝑈 siendo, en general, más bajos para el grupo de los pacientes con esquizofrenia. 36 Figura 21: Parámetros A del modelo FMM y valor de 𝑡𝑈 estimado con este mismo modelo para los genes core del grupo de esquizofrenia en el conjunto de control y en el de pacientes con esquizofrenia. La figura 22 permite la comparación entre el grupo de control y el de pacientes con esquizofrenia, de los valores del parámetro  y de R2 para el modelo FMM para los genes core del conjunto de datos de esquizofrenia. Se observa como los valores de R2 son en general más bajos para los pacientes del grupo de control. Hay que tener en cuenta que esto es lógico puesto que en este gráfico se están considerando los genes core para el grupo de esquizofrenia. 37 . La Figura 23 permite la comparación entre el grupo de control y los pacientes esquizofrénicos de los valores de R2 del conjunto de los genes core del grupo de esquizofrenia. Se observa una disminución en general de este valor en el grupo control ya que la mayor parte de los valores están por debajo de la diagonal del gráfico. Al igual que en el gráfico anterior, esta disminución es esperable puesto que los genes que se están considerando son los genes core del grupo de esquizofrenia. Figura 22: Parámetro  y valor de R2 del modelo FMM para los genes core del grupo de esquizofrenia en el conjunto de control y en el de pacientes con esquizofrenia. 38 Figura 23: Valores de R2 del modelo FMM para los genes core del grupo de esquizofrenia en el conjunto de control y en el de pacientes con esquizofrenia. Como últimos análisis, en las Figuras 24 y 25 se ofrecen comparaciones entre los valores 𝑡𝑈 obtenidos mediante el modelo Cosinor y el modelo FMM. En la Figura 24 aparecen los resultados correspondientes a los genes y el conjunto de control, mientras que en la 25 se observan los correspondientes al grupo de pacientes esquizofrénicos y los genes core de este mismo. Conviene recordar que este valor tiene interés desde el punto de vista biológico puesto que corresponde al momento de máxima expresión del gen que es cuando realiza su función biológica. En la Figura 24 no se observan grandes diferencias en la estimación de 𝑡𝑈 para los genes del grupo de control, pero en la figura 25 se observa que sí que existen diferencias en las estimaciones importantes de este valor en algunos de los genes (los más alejados de la diagonal) para el caso de los genes core correspondientes al grupo de esquizofrenia. 39 Figura 24: Valores de 𝑡𝑈 estimados a partir de los modelos FMM y Cosinor para los genes core del grupo de control en el mismo conjunto de datos. Figura 25: Valores de 𝑡𝑈 estimados a partir de los modelos FMM y Cosinor para los genes core del grupo de esquizofrenia en el conjunto de pacientes con esquizofrenia. 40 4.3 Posibles nuevos genes rítmicos En esta sección se ofrecen, a modo de ejemplo, los resultados obtenidos para algunos genes que no están en este momento catalogados como rítmicos pero que, con el nuevo procedimiento de estimación de orden considerado en este trabajo y con el modelo FMM parecen tener un patrón claramente rítmico. Pensamos que puede ser interesante, desde el punto de vista biológico, investigar las funciones de los genes de este tipo que aparecen gracias a la metodología considerada en este trabajo. Gen FOXRED2 Control ORI Reducido ToD Figura 26: Gen FOXRED2 con el orden ORI Reducido y el orden ToD para individuos control. Esquizofrenia ORI Reducido ToD Figura 27: Gen FOXRED2 con el orden ORI Reducido y el orden ToD para individuos esquizofrénicos. 47 Anexo A: Código de generación de órdenes ###################### Grupo de Control ##################### #Fichero de individuos y genes Filtered_log2CPM.csv inData<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/Filtered_log2CPM.csv", header=TRUE, sep=",")[,-1] inData<-as.matrix(inData) rownames(inData)<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/Filtered_log2CPM.csv", header=TRUE, sep=",")[,1] # Fichero con la hora de muerte data1Seney2019.xlsx deathData <- read.xlsx(file="C:/Users/carlo/Desktop/TFG Estadistica/data1Seney2019.xlsx", 1,header=TRUE, sep=",") #Fichero con todas las caracteristicas CMC_MSSM-Penn-Pitt_Clinicalv2.csv clinicalData<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/CMC_MSSM- Penn-Pitt_Clinicalv2.csv", header=TRUE, sep=";")#[,1] clinicalData<-as.matrix(clinicalData) controlDeath<- subset(deathData, Dx == "Control") controlData<-merge(controlDeath, clinicalData, by="Individual_ID") #Devuelve los indices indControlSubjects<- match(controlData[,"DLPFC_RNA_Sequencing_Sample_ID"],colnames(inData)) dataControlSubjects<-inData[,indControlSubjects] controlData<- as.matrix(controlData) ####################### ORI ########################### ORI_order<-TSP_Euc_v7(datos=dataControlSubjects,dist_type="Euclidea", datos2p=dataControlSubjects,pesosE=FALSE,pesosMSE=TRUE, pesosAdd=FALSE,unPeriodo=TRUE, pesos=rep(1,ncol(dataControlSubjects)), centrar=FALSE,intentos=25,onlyHeuristica2=FALSE) indORIData<-dataControlSubjects[,ORI_order[[1]]] dim(dataControlSubjects) ncol(dataControlSubjects) time<-rescale(rep(1:104,1),to=c(0, 2 * pi)) periodo<-24 #################### ORI Reducido ######################### geneNames<- c("ARNTL","NPAS2","CLOCK","NFIL3","CRY1","NR1D1","BHLHE41","NR1D2","DBP","CIAR T","PER1","PER3","TEF","HLF","CRY2","PER2") geneSelectionData<- dataControlSubjects[geneNames,] ORI_order_Reduced<-TSP_Euc_v7(datos=geneSelectionData,dist_type="Euclidea", datos2p=dataControlSubjects,pesosE=FALSE,pesosMSE=TRUE, pesosAdd=FALSE,unPeriodo=TRUE, pesos=rep(1,ncol(gn)), 48 centrar=FALSE,intentos=15,onlyHeuristica2=FALSE) indORIDataReduced<-dataControlSubjects[,ORI_order_Reduced[[1]]] time<-rescale(rep(1:104,1),to=c(0, 2 * pi)) periodo<-24 ##################### Orden ZT – ToD ###################### orderZT<- controlData[order(as.numeric(controlData[,"ZeitgeberTime"]) ,decreasing=FALSE),] indOrder<-match(orderZT[ ,"DLPFC_RNA_Sequencing_Sample_ID" ],colnames(inData)) indOrderData<-inData[,indOrder] timeZT<-as.numeric((orderZT[,"ZeitgeberTime"]) ###################### Grupo SCZ ######################### #Fichero de individuos y genes Filtered_log2CPM.csv inDataSCZ<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/Filtered_log2CPM.csv", header=TRUE, sep=",")[,-1] inDataSCZ<-as.matrix(inDataSCZ) rownames(inDataSCZ)<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/Filtered_log2CPM.csv", header=TRUE, sep=",")[,1] #Fichero con todas las caracteristicas CMC_MSSM-Penn-Pitt_Clinicalv2.csv clinicalDataSCZ<-read.csv(file="C:/Users/carlo/Desktop/TFG Estadistica/CMC_MSSM-Penn-Pitt_Clinicalv2.csv", header=TRUE, sep=";")#[,1] dim(clinicalDataSCZ) # Fichero con la hora de muerte y el ID ZTSCZ <- read.xlsx(file="C:/Users/carlo/Desktop/TFG Estadistica/tod_sz_control.xlsx", 1,header=TRUE, sep=",") dim(ZTSCZ)# 92 5 ZTSCZ<- subset(ZTSCZ, Dx == "SCZ") dim(ZTSCZ)# 46 5 SCZDeath<- subset(clinicalDataSCZ, Dx == "SCZ") dim(SCZDeath)# 275 18 #Devuelve los indices indSCZSubjects<- match(SCZDeath[,"DLPFC_RNA_Sequencing_Sample_ID"],colnames(inDataSCZ), na.omit) dataSCZSubjects<-inDataSCZ[,indSCZSubjects] dataSCZSubjects<-as.data.frame(dataSCZSubjects) dataSCZSubjects <- dataSCZSubjects[,colSums(is.na(dataSCZSubjects))<nrow(dataSCZSubjects)] dim(dataSCZSubjects) dataSCZSubjects<-as.matrix(dataSCZSubjects) ############################ ORI ######################### ORI_order_SCZ<-TSP_Euc_v7(datos=dataSCZSubjects,dist_type="Euclidea", datos2p=dataSCZSubjects,pesosE=FALSE,pesosMSE=TRUE, pesosAdd=FALSE,unPeriodo=TRUE, pesos=rep(1,ncol(gn)), centrar=FALSE,intentos=25,onlyHeuristica2=FALSE) ORI_order_SCZ<-ORI_order_ReducedSCZ 49 indORIDataSCZ<-dataSCZSubjects[,ORI_order_SCZ[[1]]] timeSCZ<-rescale(rep(1:46,1),to=c(0, 2 * pi)) ########################### ORI Reducido 20 ################ geneNames<-c("CIART", "WNT10B","LAMB3", "OPRL1", "CYB561", "HDAC8", "NIM1K", "DUBR", "KRT17P1", "EBP", "PGBD2", "CTSK","ZBTB22","NPRL2","IFT122","NFATC4","RNF112","VOPP1","NIT1","USF1") geneSelectionDataSCZ<- dataSCZSubjects[geneNames,] geneSelectionDataSCZ<-data.matrix(geneSelectionDataSCZ) ORI_order_ReducedSCZ<- TSP_Euc_v7(datos=geneSelectionDataSCZ,dist_type="Euclidea", datos2p=dataSCZSubjects,pesosE=FALSE,pesosMSE=TRUE, pesosAdd=FALSE,unPeriodo=TRUE, pesos=rep(1,ncol(gn)), centrar=FALSE,intentos=25,onlyHeuristica2=FALSE) indORIDataReducedSCZ<-dataSCZSubjects[,ORI_order_ReducedSCZ[[1]]] timeSCZ<-rescale(rep(1:46,1),to=c(0, 2 * pi)) ##################### Orden ZT - ToD ####################### mergeData<- merge(orderZTSCZ,SCZDeath, by="Individual_ID") mergeData<- mergeData[order(as.numeric(mergeData[,"ZeitgeberTime"]) ,decreasing=FALSE),] indicesOrderSCZ<-match(mergeData[ ,"DLPFC_RNA_Sequencing_Sample_ID" ],colnames(dataSCZSubjects)) indOrderSCZ<-dataSCZSubjects[,indicesOrderSCZ] timeZTSCZ<-mergeData[,"ZeitgeberTime"] 50 Anexo B: Ajuste de los modelos ################## GRUPO CONTROL ##################### #################### ORI ########################### indORIData<-dataControlSubjects[,orderORIControl] errorORI<- calculoError(indORIData ,nrow(indORIData),104, time, 24)#, "NP_ORI_Control", "Cos_ORI_Control", "Cos_ORI_Control" ) colnames(errorORI$error)<- c("gen","ORI_Control_NP_Error","ORI_Control_Cos_Error", "ORI_Control_FMM_Error") colnames(errorORI$R2)<-c("gen","ORI_Control_NP_R2","ORI_Control_Cos_R2", "ORI_Control_FMM_R2") colnames(errorORI$parametros_Cosinor)<- c("gen","ORI_Control_Cosinor_M","ORI_Control_Cosinor_A", "ORI_Control_Cosinor_phi") colnames(errorORI$parametros_FMM)<- c("gen","ORI_Control_FMM_M","ORI_Control_FMM_A", "ORI_Control_FMM_alpha","ORI_Control_FMM_beta","ORI_Control_FMM_omega") colnames(errorORI$pico_Cosinor)<- c("gen","ORI_Control_Cosinor_ZL","ORI_Control_Cosinor_ZU", "ORI_Control_Cosinor_TL", "ORI_Control_Cosinor_TU") colnames(errorORI$pico_FMM)<- c("gen","ORI_Control_FMM_ZL","ORI_Control_FMM_ZU", "ORI_Control_FMM_TL", "ORI_Control_FMM_TU") ####################### ORI Reducido ###################### indORIDataReduced<-dataControlSubjects[,orderORIReducedControl] errorORIReduced<- calculoError(indORIDataReduced ,nrow(indORIDataReduced),104, time, 24)#, "NP_ORI_Reducido_Control", "Cos_ORI_Reducido_Control", "FMM_ORI_Reducido_Control") colnames(errorORIReduced$error)<- c("gen","ORI_Reducido_Control_NP_Error","ORI_Reducido_Control_Cos_Error", "ORI_Reducido_Control_FMM_Error") colnames(errorORIReduced$R2)<- c("gen","ORI_Reducido_Control_NP_R2","ORI_Reducido_Control_Cos_R2", "ORI_Reducido_Control_FMM_R2") colnames(errorORIReduced$parametros_Cosinor)<- c("gen","ORI_Reducido_Control_Cosinor_M","ORI_Reducido_Control_Cosinor_A", "ORI_Reducido_Control_Cosinor_phi") colnames(errorORIReduced$parametros_FMM)<- c("gen","ORI_Reducido_Control_FMM_M","ORI_Reducido_Control_FMM_A", "ORI_Reducido_Control_FMM_alpha","ORI_Reducido_Control_FMM_beta","ORI_Reducido _Control_FMM_omega") colnames(errorORIReduced$pico_Cosinor)<- c("gen","ORI_Reducido_Control_Cosinor_ZL","ORI_Reducido_Control_Cosinor_ZU", "ORI_Reducido_Control_Cosinor_TL", "ORI_Reducido_Control_Cosinor_TU") colnames(errorORIReduced$pico_FMM)<- c("gen","ORI_Reducido_Control_FMM_ZL","ORI_Reducido_Control_FMM_ZU", "ORI_Reducido_Control_FMM_TL", "ORI_Reducido_Control_FMM_TU") 51 ###################### ZT - ToD ####################### indOrderData<-inData[,indOrder] errorZT<- calculoError(indOrderData ,nrow(indOrderData),104, escalado(timeZT), 24)#, "NP_ZT_Control", "Cos_ZT_Control", "FMM_ZT_Control") colnames(errorZT$error)<-c("gen","ZT_Control_NP_Error","ZT_Control_Cos_Error", "ZT_Control_FMM_Error") colnames(errorZT$R2)<-c("gen","ZT_Control_NP_R2","ZT_Control_Cos_R2", "ZT_Control_FMM_R2") colnames(errorZT$parametros_Cosinor)<- c("gen","ZT_Control_Cosinor_M","ZT_Control_Cosinor_A", "ZT_Control_Cosinor_phi") colnames(errorZT$parametros_FMM)<- c("gen","ZT_Control_FMM_M","ZT_Control_FMM_A", "ZT_Control_FMM_alpha","ZT_Control_FMM_beta","ZT_Control_FMM_omega") colnames(errorZT$pico_Cosinor)<- c("gen","ZT_Control_Cosinor_ZL","ZT_Control_Cosinor_ZU", "ZT_Control_Cosinor_TL", "ZT_Control_Cosinor_TU") colnames(errorZT$pico_FMM)<-c("gen","ZT_Control_FMM_ZL","ZT_Control_FMM_ZU", "ZT_Control_FMM_TL", "ZT_Control_FMM_TU") ################## FICHERO ###################### error_Control<-merge(errorORI$error, merge(errorORIReduced$error,errorZT$error, by="gen" ),by="gen") R2_Control<-merge(errorORI$R2, merge(errorORIReduced$R2,errorZT$R2, by="gen" ),by="gen") Cosinor_Parametros_Control<-merge(errorORI$parametros_Cosinor, merge(errorORIReduced$parametros_Cosinor,errorZT$parametros_Cosinor, by="gen" ),by="gen") FMM_Parametros_Control<-merge(errorORI$parametros_FMM, merge(errorORIReduced$parametros_FMM,errorZT$parametros_FMM, by="gen" ),by="gen") Cosinor_Pico_Control<-merge(errorORI$pico_Cosinor, merge(errorORIReduced$pico_Cosinor,errorZT$pico_Cosinor, by="gen" ),by="gen") FMM_Pico_Control<-merge(errorORI$pico_FMM, merge(errorORIReduced$pico_FMM,errorZT$pico_FMM, by="gen" ),by="gen") error_R2_Control<- merge(error_Control,R2_Control,by="gen" ) parametros_Control<- merge(Cosinor_Parametros_Control,FMM_Parametros_Control,by="gen" ) picos_Control<- merge(Cosinor_Pico_Control,FMM_Pico_Control,by="gen" ) total_Control<- merge(error_R2_Control,merge(parametros_Control,picos_Control , by="gen" ),by="gen" ) ################ GRUPO SCZ ################### ################### ORI ########################### indORIDataSCZ<-dataSCZSubjects[,orderORISCZ] errorORISCZ<- calculoError(indORIDataSCZ ,nrow(indORIDataSCZ),46, timeSCZ, 24) colnames(errorORISCZ$error)<-c("gen","ORI_SCZ_NP_Error","ORI_SCZ_Cos_Error", "ORI_SCZ_FMM_Error") colnames(errorORISCZ$R2)<-c("gen","ORI_SCZ_NP_R2","ORI_SCZ_Cos_R2", "ORI_SCZ_FMM_R2") colnames(errorORISCZ$parametros_Cosinor)<- c("gen","ORI_SCZ_Cosinor_M","ORI_SCZ_Cosinor_A", "ORI_SCZ_Cosinor_phi") 52 colnames(errorORISCZ$parametros_FMM)<-c("gen","ORI_SCZ_FMM_M","ORI_SCZ_FMM_A", "ORI_SCZ_FMM_alpha","ORI_SCZ_FMM_beta","ORI_SCZ_FMM_omega") colnames(errorORISCZ$pico_Cosinor)<- c("gen","ORI_SCZ_Cosinor_ZL","ORI_SCZ_Cosinor_ZU", "ORI_SCZ_Cosinor_TL", "ORI_SCZ_Cosinor_TU") colnames(errorORISCZ$pico_FMM)<-c("gen","ORI_SCZ_FMM_ZL","ORI_SCZ_FMM_ZU", "ORI_SCZ_FMM_TL", "ORI_SCZ_FMM_TU") ############## ORI Reducido ###################### indORIDataReducedSCZ<-dataSCZSubjects[,orderORIReducedSCZ] errorORIReducedSCZ<- calculoError(indORIDataReducedSCZ ,nrow(indORIDataReducedSCZ),46, timeSCZ, 24) #, "NP_ORI_Reducido_Control", "Cos_ORI_Reducido_Control", "FMM_ORI_Reducido_Control") colnames(errorORIReducedSCZ$error)<- c("gen","ORI_Reducido_SCZ_NP_Error","ORI_Reducido_SCZ_Cos_Error", "ORI_Reducido_SCZ_FMM_Error") colnames(errorORIReducedSCZ$R2)<- c("gen","ORI_Reducido_SCZ_NP_R2","ORI_Reducido_SCZ_Cos_R2", "ORI_Reducido_SCZ_FMM_R2") colnames(errorORIReducedSCZ$parametros_Cosinor)<- c("gen","ORI_Reducido_SCZ_Cosinor_M","ORI_Reducido_SCZ_Cosinor_A", "ORI_Reducido_SCZ_Cosinor_phi") colnames(errorORIReducedSCZ$parametros_FMM)<- c("gen","ORI_Reducido_SCZ_FMM_M","ORI_Reducido_SCZ_FMM_A", "ORI_Reducido_SCZ_FMM_alpha","ORI_Reducido_SCZ_FMM_beta","ORI_Reducido_SCZ_FMM _omega") colnames(errorORIReducedSCZ$pico_Cosinor)<- c("gen","ORI_Reducido_SCZ_Cosinor_ZL","ORI_Reducido_SCZ_Cosinor_ZU", "ORI_Reducido_SCZ_Cosinor_TL", "ORI_Reducido_SCZ_Cosinor_TU") colnames(errorORIReducedSCZ$pico_FMM)<- c("gen","ORI_Reducido_SCZ_FMM_ZL","ORI_Reducido_SCZ_FMM_ZU", "ORI_Reducido_SCZ_FMM_TL", "ORI_Reducido_SCZ_FMM_TU") ####################### ZT - ToD ###################### indOrderSCZ<-dataSCZSubjects[,indicesOrderSCZ] timeZTSCZ<-mergeData[,"ZeitgeberTime"] errorZTSCZ<- calculoError(indOrderSCZ ,nrow(indOrderSCZ),46, escalado(timeZTSCZ), 24) #, "NP_ZT_Control", "Cos_ZT_Control", "FMM_ZT_Control") colnames(errorZTSCZ$error)<-c("gen","ZT_SCZ_NP_Error","ZT_SCZ_Cos_Error", "ZT_SCZ_FMM_Error") colnames(errorZTSCZ$R2)<-c("gen","ZT_SCZ_NP_R2","ZT_SCZ_Cos_R2", "ZT_SCZ_FMM_R2") colnames(errorZTSCZ$parametros_Cosinor)<- c("gen","ZT_SCZ_Cosinor_M","ZT_SCZ_Cosinor_A", "ZT_SCZ_Cosinor_phi") colnames(errorZTSCZ$parametros_FMM)<-c("gen","ZT_SCZ_FMM_M","ZT_SCZ_FMM_A", "ZT_SCZ_FMM_alpha","ZT_SCZ_FMM_beta","ZT_SCZ_FMM_omega") colnames(errorZTSCZ$pico_Cosinor)<- c("gen","ZT_SCZ_Cosinor_ZL","ZT_SCZ_Cosinor_ZU", "ZT_SCZ_Cosinor_TL", "ZT_SCZ_Cosinor_TU") colnames(errorZTSCZ$pico_FMM)<-c("gen","ZT_SCZ_FMM_ZL","ZT_Control_FMM_ZU", "ZT_SCZ_FMM_TL", "ZT_SCZ_FMM_TU") 53 ######################## FICHERO ################## error_SCZ<-merge(errorORISCZ$error, merge(errorORIReducedSCZ$error,errorZTSCZ$error, by="gen" ),by="gen") R2_SCZ<-merge(errorORISCZ$R2, merge(errorORIReducedSCZ$R2,errorZTSCZ$R2, by="gen" ),by="gen") Cosinor_Parametros_SCZ<-merge(errorORISCZ$parametros_Cosinor, merge(errorORIReducedSCZ$parametros_Cosinor,errorZTSCZ$parametros_Cosinor, by="gen" ),by="gen") FMM_Parametros_SCZ<-merge(errorORISCZ$parametros_FMM, merge(errorORIReducedSCZ$parametros_FMM,errorZTSCZ$parametros_FMM, by="gen" ),by="gen") Cosinor_Pico_SCZ<-merge(errorORISCZ$pico_Cosinor, merge(errorORIReducedSCZ$pico_Cosinor,errorZTSCZ$pico_Cosinor, by="gen" ),by="gen") FMM_Pico_SCZ<-merge(errorORISCZ$pico_FMM, merge(errorORIReducedSCZ$pico_FMM,errorZTSCZ$pico_FMM, by="gen" ),by="gen") error_R2_SCZ<- merge(error_SCZ,R2_SCZ,by="gen" ) parametros_SCZ<- merge(Cosinor_Parametros_SCZ,FMM_Parametros_SCZ,by="gen" ) picos_SCZ<- merge(Cosinor_Pico_SCZ,FMM_Pico_SCZ,by="gen" ) total_SCZ<- merge(error_R2_SCZ,merge(parametros_SCZ,picos_SCZ , by="gen" ),by="gen" ) 54 Anexo C: Funciones y librerías library(Iso) library(TSP) library(foreach) library(iterators) library(parallel) library(doParallel) library(bigstatsr) library(scales) library(xlsx) library(FMM) library(ggplot2) library(cosinor2) library(ggplot2) library(gridExtra) library(grid) library("writexl") library("grid") library("ggplotify") cosinor.lm <- function(formula, period = 12, dato, na.action = na.omit){ # build time tranformations Terms <- terms(formula, specials = c("time", "amp.acro")) stopifnot(attr(Terms, "specials")$time != 1) varnames <- get_varnames(Terms) timevar <- varnames[attr(Terms, "specials")$time - 1] xx<-cos(dato[,timevar]) zz<-sin(dato[,timevar]) fit<-lm((dato[,1])~xx+zz) Mest<-fit$coefficients[1] bb<-fit$coefficients[2] gg<-fit$coefficients[3] phiEst<-atan2(-gg,bb)%%(2*pi) Aest<-sqrt(bb^2+gg^2) rss<- sum (( dato[,1] - (Mest+bb*cos(2*pi*time/periodo)- gg*sin(2*pi*time/periodo)))^2) #data$rrr <- cos(2 * pi * dato[,timevar] / period) #data$sss <- sin(2 * pi * dato[,timevar] / period) data$rrr <- xx data$sss <- zz spec_dex <- unlist(attr(Terms, "special")$amp.acro) - 1 55 mainpart <- c(varnames[c(-spec_dex, - (attr(Terms, "special")$time - 1))], "rrr", "sss") acpart <- paste(sort(rep(varnames[spec_dex], 2)), rep(c("rrr", "sss"), length(spec_dex)), sep = ":") newformula <- as.formula(paste(rownames(attr(Terms, "factors"))[1], paste(c(mainpart, acpart), collapse = " + "), sep = " ~ ")) fit <- lm(newformula, data, na.action = na.action) mf <- fit r.coef <- c(FALSE, as.logical(attr(mf$terms, "factors")["rrr",])) s.coef <- c(FALSE, as.logical(attr(mf$terms, "factors")["sss",])) mu.coef <- c(TRUE, ! (as.logical(attr(mf$terms, "factors")["sss",]) | as.logical(attr(mf$terms, "factors")["rrr",]))) beta.s <- mf$coefficients[s.coef] beta.r <- mf$coefficients[r.coef] groups.r <- c(beta.r["rrr"], beta.r["rrr"] + beta.r[which(names(beta.r) != "rrr")]) groups.s <- c(beta.s["sss"], beta.s["sss"] + beta.s[which(names(beta.s) != "sss")]) amp <-Aest #amp <- sqrt(groups.r^2 + groups.s^2) names(amp) <- gsub("rrr", "amp", names(beta.r)) #acr <-phiEst acr <- atan(groups.s / groups.r) # print("acrofase") # print(atan(groups.s / groups.r)) names(acr) <- gsub("sss", "acr", names(beta.s)) coef <- c(mf$coefficients[mu.coef], amp, acr) #coef <- c(Mest, amp, acr) structure(list(fit = fit, Call = match.call(), Terms = Terms, coefficients = coef, period = period, rss=rss), class = "cosinor.lm") } get_varnames <- function(Terms){ spec <- names(attr(Terms, "specials")) tname <- attr(Terms, "term.labels") dex <- unlist(sapply(spec, function(sp){ attr(Terms, "specials")[[sp]] - 1 56 })) tname2 <- tname for(jj in spec){ gbl <- grep(paste0(jj, "("), tname2, fixed = TRUE) init <- length(gbl) > 0 if( init ){ jlack <- gsub(paste0(jj, "("), "", tname2, fixed = TRUE) tname2[gbl] <- substr(jlack[gbl], 1, nchar(jlack[gbl]) - 1) } } tname2 } update_covnames <- function(names){ covnames <- grep("(amp|acr|Intercept)", names, invert = TRUE, value = TRUE) lack <- names for(n in covnames){ lack <- gsub(paste0(n, ":"), paste0("[", n, " = 1]:"), lack) lack <- gsub(paste0("^", n, "$"), paste0("[", n, " = 1]"), lack) } lack } ggplot.cosinor.lm <- function(object,l, x_str = NULL){ timeax <- seq(0, object$period, length.out = l) covars <- grep("(rrr|sss)", attr(object$fit$terms, "term.labels"), invert = TRUE, value = TRUE) newdata <- data.frame(time = timeax, rrr = cos(2 * pi * timeax / object$period), sss = sin(2 * pi * timeax / object$period)) for(j in covars){ newdata[,j] <- 0 } if(!is.null(x_str)){ for(d in x_str){ tdat <- newdata tdat[,d] <- 1 newdata <- rbind(newdata, tdat) } newdata$levels <- "" for(d in x_str){ newdata$levels <- paste(newdata$levels, paste(d, "=", newdata[,d]))