Metodología
2023-06-15
1 MÉTODOS
1.1 Datos genómicos
Se descargaron 2 conjuntos de genomas de cianobacterias del servidor del NCBI (https://www.ncbi.nlm.nih.gov/genome/browse/#!/prokaryotes).
Estos conjuntos corresponden a:
- 269 genomas completos y aquellos que solo contenian el cromosoma (complete_chr)
- 165 genomas nuevos usados en Cabello-Yeves et al. (2022) (pico)
Dichos genomas fueron descargados en formato Genebank (.gbk o gbff).
1.2 Abundancia de palíndromos.
Una vez descargados los genomas, el siguiente paso fue calcular el valor observado y esperado de repeticiones de todos los posibles octámeros palindrómicos de 8 nucleótidos.
El valor observado es el número de veces que cada octámero palindrómico se repite a lo largo de cada genoma. El valor esperado se calculó mediante un modelo de markov de 3er orden.
1.2.1 Modelos de Markov
En una cadena de Markov, el valor tomado por una variable aleatoria depende de los valores tomados por la variable aleatoria en un estado anterior. El número de estados históricos que influyen en el valor de la variable aleatoria en un lugar dado a lo largo de la secuencia también se conoce como el grado del proceso de Markov. El modelo de cadena de Markov de primer grado tiene parámetros \(|\Sigma| + |\Sigma|^2\), correspondientes a las frecuencias de nucleótidos individuales así como a las frecuencias de dinucleótidos. De esta manera, este modelo permite que una posición sea dependiente de la posición anterior. Sin embargo, las frecuencias se modelan de manera invariable en la posición y, por lo tanto, pueden no ser adecuadas para modelar señales. Este modelo de secuencia \(M\) se define sobre el espacio muestral \(\Sigma^{*}\) y asigna una probabilidad a cada secuencia \(x\) de longitud \(n(x)\) sobre \(\Sigma^{*}\):
\[\begin{equation} P(x|M) = P_1(x_1) \prod_{i=2,...,n(x)} P_2(x_i|x_{i-1},...,x_{i-n}) \tag{1.1} \end{equation}\]
donde \(P_1\) es una función de probabilidad en \(\Sigma\) que modela la distribución de \(\alpha\)’s en la primera posición de la secuencia y \(P_2\) es la función de probabilidad condicional en \(\Sigma\times\Sigma\) que modela la distribución de \(\beta\)’s en la posición \(i>1\) en el símbolo alfabético \(\alpha\) en la posición \(i-1\). La estimación de parámetros se hace utilizando el estimador de probabilidad máxima. Las probabilidades de transición se estiman utilizando el teorema de Bayes, como se muestra a continuación:
\[\begin{equation} P_2(\beta|\alpha)={P(\alpha\beta) \over P(\alpha)} \tag{1.2} \end{equation}\]
De esta manera, las probabilidades transicionales condicionales de encontrar una base \(\beta\) en la posición (\(i\)) dado que la base \(\alpha\) se encontró en la posición (\(i-1\)) se calculan encontrando la abundancia del dinucleótido \(\alpha\beta\) como una fracción de la abundancia del nucleótido \(\alpha\).
Ejemplo:
Considerando una la secuencia de 25 nucleótidos. \[Seq = AACGT\space CTCTA\space TCATG\space CCAGG\space ATCTG\]
Al considerar los modelos de cadena de Markov de primer grado, es necesario calcular los \(4-parámetros\) correspondientes a las frecuencias de nucleótidos individuales y los \(4^2\) parámetros correspondientes a las frecuencias de dinucleótidos. Los parámetros de \(\Sigma\) son:
\[\begin{equation} \begin{split} \Sigma & = \{frec(A),frec(C),frec(G),frec(T),\}\\ & =\{\frac{6}{25},\frac{7}{25},\frac{7}{25},\frac{5}{25}\} \end{split} \tag{1.3} \end{equation}\]
Para calcular \(P_2\), los valores de probabilidad condicional \(\Sigma \times \Sigma\), las frecuencias de dinucleótidos y las probabilidades se calculan a partir de los datos de secuencia. Las frecuencias de los dinucleótidos y las probabilidades se muestran a continuación (con los números entre paréntesis que representan las probabilidades):
\[\begin{equation} \Sigma \times \Sigma = \begin{Bmatrix} frec(AA)=\frac{1}{24} & frec(AC)=\frac{1}{24} & frec(AT)=\frac{3}{24} & frec(AG)=\frac{1}{24} \\ frec(CA)=\frac{2}{24} & frec(CC)=\frac{1}{24} & frec(CT)=\frac{3}{24} &frec(CG)=\frac{1}{24} \\ frec(TA)=\frac{1}{24} & frec(TC)=\frac{4}{24} & frec(TT)=\frac{0}{24} & frec(TG)=\frac{1}{24} \\ frec(GA)=\frac{1}{24} & frec(GC)=\frac{1}{24} & frec(GT)=\frac{1}{24} & frec(GG)=\frac{1}{24} \end{Bmatrix} \tag{1.4} \end{equation}\]
A continuación, las probabilidades condicionales se calculan utilizando el teorema de Bayes (consulte la Ecuación (1.2) ). Por ejemplo, la probabilidad de encontrar \(C\) en la posición \(i+1\) dado que se ha encontrado una \(A\) en la posición (\(i\)) es:
\[\begin{equation} P(C|A)= \frac{P_{AC}}{P_A}=\frac{\frac{1}{24}}{\frac{6}{25}} \tag{1.5} \end{equation}\]
Para secuencias grandes, la probabilidad condicional \(P(S_i|S_{i-1})\) se aproxima a:
\[\begin{equation} P(S_i|S_{i-1}) =\frac{frec(S_iS_{i-1})}{frec(S_{i-1})} \tag{1.6} \end{equation}\]
Las probabilidades condicionales para la secuencia de ejemplo se muestran en (1.4). Usando estos parámetros del modelo, la probabilidad de encontrar el patrón \(CAAT\) en esta secuencia usando el modelo de Markov de primer orden de la secuencia subyacente sería igual a:
\[\begin{equation} \begin{split} P(C)P(A|C)P(A|A)P(T|A) & = P(C)\cdot\frac{P(CA)}{P(C)}\cdot\frac{P(AA)}{P(A)}\cdot\frac{P(AT)}{P(A)} \\ & = (\frac{7}{25}) \cdot (\frac{50}{168}) \cdot (\frac{25}{144}) \cdot (\frac{75}{144})\\ & = 0.0075 \end{split} \tag{1.7} \end{equation}\]
1.2.2 Modelos de Markov de orden superior
Una cadena de Markov de orden \(n\) sobre algún alfabeto \(A\) es equivalente a una cadena de Markov de primer orden sobre el alfabeto \(A^n\) de \(n\)-tuplas. Esto se sigue del cálculo de la probabilidad de \(A\) y \(B\), dado \(B\) es la probabilidad de \(A\) dado \(B\), es decir: \[\begin{equation} P(x_k|x_{k-1},...,x_{k-n}) = P(x_k,x_{k-1},...,x_{k-n+1}|x_{k-1},...,x_{k-n}) \tag{1.8} \end{equation}\] es decir, la probabilidad de \(x_k\) dada la \(n\)-tupla terminada en \(x_{k-1}\) es igual a la probabilidad de que la \(n\)-tupla termine en \(x_k\) dado que la \(n\)-tupla termina en \(x_{k-1}\).
1.2.3 Modelo de Markov de orden 1 para hallar octanucleótidos
Por ejemplo, para una octanucleótido de 8 letras, digamos HIP1: \[W=GCGATCGC\] Los parametros de \(\Sigma\) corresponden a: \[\begin{equation} \Sigma= \{frec(A),frec(C),frec(G),frec(T)\} \end{equation}\]
Los valores de probabilidad condicional de \(\Sigma \times \Sigma\) son: \[\begin{equation} \Sigma \times \Sigma = \begin{Bmatrix} frec(AA) & frec(AC) & frec(AT) & frec(AG) \\ frec(CA) & frec(CC) & frec(CT) & frec(CG) \\ frec(TA) & frec(TC) & frec(TT) & frec(TG) \\ frec(GA) & frec(GC) & frec(GT) & frec(GG) \end{Bmatrix} \end{equation}\]
Si queremos usar un modelo de orden 1, la probabilidad de hallar \(W\) segun las ecuaciones (1.1) y(1.2) es: \[\begin{equation} \begin{split} P(W) & = P(G) \cdot P(C|G) \cdot P(G|C) \cdot P(A|G) \cdot P(T|A) \cdot P(C|T) \cdot P(G|C) \cdot P(C|G)\\ & = {P(G)} \cdot \frac{P(GC)}{P(G)} \cdot \frac{P(CG)}{P(C)} \cdot \frac{P(GA)}{P(G)} \cdot \frac{P(AT)}{P(A)} \cdot \frac{P(TC)}{P(T)} \cdot \frac{P(CG)}{P(C)} \cdot \frac{P(GC)}{P(G)}\\ & = P(GC) \cdot \frac{P(CG)}{P(C)} \cdot \frac{P(GA)}{P(G)} \cdot \frac{P(AT)}{P(A)} \cdot \frac{P(TC)}{P(T)} \cdot \frac{P(CG)}{P(C)} \cdot \frac{P(GC)}{P(G)} \end{split} \end{equation}\]
finalmente: \[\begin{equation} P(W) = \frac{P(GC) \cdot P(CG) \cdot P(GA) \cdot P(AT) \cdot P(TC) \cdot P(CG) \cdot P(GC)}{P(C) \cdot P(G) \cdot P(A) \cdot P(T) \cdot P(C) \cdot P(G)} \tag{1.9} \end{equation}\]
1.3 Significancia de los conteos observados
Para darle una significancia estadística al conteo se usó una prueba binomial y un test FDR.
1.3.1 Prueba binomial.
Para calcular la probabilidad de que el conteo esperado, el cual sigue una distribución binomial, tome valores MAYORES O IGUALES al conteo observado, usamos la función pbinom
Donde:
- q: Cuantil o vector de cuantiles
- size: Numero de experimentos (n>=0)
- prob: Probabilidad de éxito en cada experimento
- lower.tail: si es TRUE, las probabilidades son P(X<=x), o P(X>x) en otro caso.
Tomemos un caso particular del conteo:
| Spp | Palindrome | Observed | Markov (Expected) | GenomeSize |
|---|---|---|---|---|
| 336-3 | GCGATCGC | 6202 | 65.396286071305 | 6420126 |
La probabilidad de que se observen 6202 sitios \(GCGATCGC\), O MAS, si el número de sitios posibles en el genoma es 6420119 (\(6420126-8+1\), es decir \(GenomeSize-k+1\)) y la probabilidad de observar dicho sitio es de: 1.018615e-05 (\(65.3962860713054 \over 6420126-8+1\), es decir \(Expected \over GenomeSize-k+1\)), es casi 0.
En otras palabras, la probabilidad de que suceda lo que estoy observando es muy baja.
1.3.2 FDR
Para estudios en los que se realizan miles de test de forma simultánea, el resultado de estos métodos es demasiado conservativo e impide que se detecten diferencias reales. Una alternativa es controlar el false discovery rate o FDR.
Para nuestros datos el FDR se calculó en R de usando los valores obtenidos de la prueba binomial:
Donde pval es la probabilidad obtenida de la prueba binomial.
1.3.3 Conjuntos de conteos de acuerdo a la significancia
Se crearon 4 conjuntos de resultados de acuerdo a 4 valores mínimos de significancia de acuerdo al FDR:
- sel32 (\(1 \times 10^{-32}\))
- sel64 (\(1 \times 10^{-64}\))
- sel128 (\(1 \times 10^{-128}\))
- sel256 (\(1 \times 10^{-256}\))
El conjunto más laxo corresponde a sel32 ya que su valor de corte de FDR es \(1 \times 10^{-32}\), debido a esto, es el conjunto con más palíndromos (Figure 1.1). Por otro lado, el conjunto sel256 es el conjunto más restrictivo ya que su valor de corte de FDR es de \(1 \times 10^{-256}\), y por lo tanto tiene menos palíndromos (Figure 1.2).
1.4 Visualización de la abundancia: OE vs Frecuencia Observada cada 1000nt
Para visualizar la abundancia creamos un gráfico que muestra el enriquecimiento OE vs la abundancia por cada 1000 nucleótidos. Esto se hizo para cada conjunto de significancia y para cada conjunto de genomas.
Figure 1.1: Enriquecimiento versus abundancia de palíndromos octámeros en el conjunto de genomas complete_chr con un \(FDR \leq 1 \times 10^{-32}\). Enriquecimiento (O/E) en función de la frecuencia del motivo cada 1000 nt (FrecObs). Cada punto representa un palíndromo octámero de un genoma.
Figure 1.2: Enriquecimiento versus abundancia de palíndromos octámeros en el conjunto de genomas complete_chr con un \(FDR \leq 1 \times 10^{-256}\). Enriquecimiento (O/E) en función de la frecuencia del motivo cada 1000 nt (FrecObs). Cada punto representa un palíndromo octámero de un genoma.
1.5 Filogenia
Se infirieron filogenias para los dos conjuntos de genomas. Para esto usamos el software Orthofinder (Emms and Kelly (2019)), el cual utiliza FastME para inferir la filogenia (Lefort, Desper, and Gascuel (2015)). FastME proporciona algoritmos de distancia para inferir filogenias. FastME se basa en una evolución mínima equilibrada, que es el principio mismo de Neighbor Joining (NJ).
El software se corrió en la línea de comandos de la siguiente manera:
1.5.1 Anotación de la filogenia
Para tener una forma de más visual de entender la distribución de los palíndromos en los genomas, anotamos las filogenias de acuerdo a su abundancia. Se anotaron 4 filogenias según la significancia (sel32, sel64, sel128 y sel256) para los 2 conjuntos de genomas. Además, esta anotación se hizo para la abundancia de acuerdo a la Frecuencia Observada por cada 1000 nucleotidos (\(FrecObs\)) (Figure 1.3) y a la tasa de Observados sobre esperados (\(OE\)) (Figure 1.4).
La anotación de las filogenias consistió en agregarles un heatmap que mostrara la abundancia de cada palíndromo y un diagrama de barras que indicara aquel palíndromo con mayor abundancia.
Figure 1.3: Filogenia del conjunto de genomas complete_chr anotada de acuerdo a la Frecuencia observada cada 1000 nt (FrecObs). La abundancia visualizada en esta filogenia es de acuerdo al conjunto sel256, es decir conteos con un \(FDR \leq 1 \times 10^{-256}\). La filogenia muestra 269 especies, frente a la filogenia se muestra un heatmap que indica la abundancia de cada palíndromo. Frente al Heatmap se muestra un Diagrama de barras el cual indíca el palindromo mas abundante de entre todos.
Figure 1.4: Filogenia del conjunto de genomas complete_chr anotada de acuerdo a la tasa de observados sobre esperados (OE). La abundancia visualizada en esta filogenia es de acuerdo al conjunto sel256, es decir conteos con un \(FDR \leq 1 \times 10^{-256}\). La filogenia muestra 269 especies, frente a la filogenia se muestra un heatmap que indica la abundancia de cada palíndromo. Frente al Heatmap se muestra un Diagrama de barras el cual indíca el palindromo mas abundante de entre todos.
1.6 Identificación de casos relevantes
De acuerdo a las filogenias anotadas, se buscaron aquellos casos en los que HIP1 o algún otro palíndromo se hubiera ganado o perdido abruptamente y en su lugar hubiese otro palíndromo abundante. Además, se buscó que en aquellos casos, las ramas en la filogenia no fueran tan largas. Esto se hizo de manera visual revisando el diagrama de barras que mostraba el palíndromo más abundante para cada especie. En total hubo 6 subclados que mostraban cambios abruptos en la abundancia de sus palíndromos (Figure 1.5).
Figure 1.5: Casos de interés. En la figura se muestran remarcados los casos interesantes: clado calothrix (rojo), clado cyanobacterium (naranja), clado geminocystis (amarillo), clado thermosynechococcus (azul), clado pseudoanabaena (verde).
También se hallo un caso interesante en el conjunto pico (clado A18-40) el cual sirvió como punto de partida para analisis posteriores. En este caso se muestra que la especie Synechococcus A18-40 muestra una tasa OE mucho mayor comparada con las demás especies del clado (Figure 1.6).
Figure 1.6: Casos de interés. En la figura se muestra remarcado el clado A18-40 (azul).
1.7 Reconstrucción Ancestral de sitios palindrómicos en ortólogos
Para tratar de entender como es que los sitios HIP1 han ido evolucionando, hicimos una reconstrucción de sitios ancestrales y posteriormente construimos varios conjuntos de redes para visualizar dicha evolución.
1.7.1 Ortólogos
Para simplificar la reconstrucción de secuencias ancestrales usamos unicamente los ortólogos. Para obtener esto usamos el pipeline get_homologues:
Después de obtener los ortólgos filtramos:
- aquellos que no estuvieran en las 6 especies del clado
- aquellos que tuvieran mas de una copia (parálogos)
- aquellos sin sitios HIP1
1.7.2 Reconstrucción
Para hacer la reconstrucción usamos la pagueteria de R phangorn, la cual proporciona varios métodos para estimar estados de caracteres ancestrales con Máxima Parsimonia (MP) o Máxima Verosimilitud (ML). En este caso usamos ML. Adicionalmente podemos asignar los estados ancestrales según la máxima verosimilitud (“ml”):
\[P(x_r=A) = {{L(x_r=A)} \over {\sum\limits_{k\in \{A,C,G,T\}}} L(x_r= k) }\]
y el criterio de mayor probabilidad posterior (“bayes”): \[P(x_r=A) = {{\pi_A L(x_r=A)} \over {\sum\limits_{k\in \{A,C,G,T\}}} \pi_k L(x_r= k) }\]
dónde \(L(x_r)\) es la probabilidad conjunta de los estados en las puntas y el estado en la raíz \(x_r\) y \(\pi_i\) son las frecuencias base estimadas del estado \(i\).
Toda la información de la reconstrucción fue guardada en dos tablas las cuales contienen listas de cada transicion entre cada estado. Estas tablas fueron creadas con la siguiente función: