# Base de datos FooDrugs :pill: :pizza:
###### tags: `Marco`, `FNS-Cloud`, `FooDrugs`, `MySQL`, `R`, `database`
[ToC]
## Preparando la base de datos
El archivo está en la partición de Ubuntu de Marco. Lo he encontrado con estos comandos desde mi directorio home:
```bash
ls -lht ../mgarranzo | grep 'sql'
```
```
-rw-rw-r-- 1 1000 1000 1,3G ene 31 07:06 FinalFooDrugs_v2.sql
```
Vemos que el archivo es del **31 de enero de 2022**, por lo que sí parece ser la última versión de la base de datos.
### Configuración de MySQL
Ahora quedaría instalarla para tenerla en MySQL. He intentado instalar Workbench siguiendo [este tutorial](https://platzi.com/tutoriales/1566-bd/8226-como-instalar-mysql-y-workbench-en-ubuntu-sin-morir-en-el-intento/), pero me da un error porque no encuentra el cliente de línea de comandos.
Así pues, he seguido con la línea de comandos. Para **generar la base de datos:**
```bash
mysql -u root -p FooDrugs_v2 < FinalFooDrugs_v2.sql
```
:::warning
Es importante que pongamos `FooDrugs_v2` como nombre de la base de datos porque es ese el que viene en el script `FinalFooDrugs_v2.sql`:
```sql=!
CREATE DATABASE /*!32312 IF NOT EXISTS*/ `FooDrugs_v2` /*!40100 DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_0900_ai_ci */ /*!80016 DEFAULT ENCRYPTION='N' */;
```
:::
Esto para el portátil. En Totoro he tenido que crear manualmente la base de datos primero:
```bash
blanca@totoro:~$ mysql -u root -p
```
```sql
mysql> CREATE DATABASE FooDrugs_v2;
mysql> exit
```
```bash
blanca@totoro:~$ mysql -u root -p FooDrugs_v2 < FinalFooDrugs_v2.sql
```
## Buscando la información relevante
### Número de textos de cada base de datos
En la tabla `texts` tenmos la URL y el texto de cada texto. Podemos ver las URLs con:
```sql
SELECT link FROM texts;
```
Así veremos todas las URLs. Al final vemos que la tabla tiene **2849 filas**. Observaremos algunos links repetidos. Por ejemplo:
```sql=!
SELECT texts_ID FROM texts WHERE link = 'https://www.drugs.com/food-interactions/fentanyl,sublimaze.html?professional=1';
```
```
+----------+
| texts_ID |
+----------+
| 2 |
| 2931 |
| 2934 |
| 2937 |
| 2940 |
| 2943 |
+----------+
```
Si nos fijamos en la columna `document`, veremos que los textos también son idénticos:
```sql=!
SELECT document FROM texts WHERE link = 'https://www.drugs.com/food-interactions/fentanyl,sublimaze.html?professional=1';
```
Asumo que esto es porque en un mismo texto se puede mencionar varias veces cada interacción. Por ejemplo:
```sql=!
SELECT * FROM TM_interactions WHERE texts_id = '2' OR texts_id = '2931' OR texts_ID = '2934' OR texts_id = '2937' OR texts_ID = '2940' OR texts_ID = '2943';
```
```
+--------------------+----------+-------------+-----------+------------+----------+
| TM_interactions_ID | texts_ID | start_index | end_index | food | drug |
+--------------------+----------+-------------+-----------+------------+----------+
| 3 | 2 | 385 | 641 | grapefruit | fentanyl |
| 4 | 2 | 643 | 846 | grapefruit | fentanyl |
| 5 | 2 | 1953 | 2155 | grapefruit | fentanyl |
| 6407 | 2931 | 385 | 641 | grapefruit | fentanyl |
| 6408 | 2931 | 643 | 846 | grapefruit | fentanyl |
| 6409 | 2931 | 1953 | 2155 | grapefruit | fentanyl |
| 6413 | 2934 | 385 | 641 | grapefruit | fentanyl |
| 6414 | 2934 | 643 | 846 | grapefruit | fentanyl |
| 6415 | 2934 | 1953 | 2155 | grapefruit | fentanyl |
| 6419 | 2943 | 385 | 641 | grapefruit | fentanyl |
| 6420 | 2943 | 643 | 846 | grapefruit | fentanyl |
| 6421 | 2943 | 1953 | 2155 | grapefruit | fentanyl |
+--------------------+----------+-------------+-----------+------------+----------+
```
Si nos fijamos en el texto 2, veremos que en efecto hay tres menciones de la interacción entre el zumo de uva y el fentanil:
>| Major Food Interaction||GENERALLY AVOID: Alcohol may potentiate the central nervous system (CNS) depressant effects of opioid analgesics including fentanyl. Concomitant use may result in additive CNS depression and impairment of judgment, thinking, and psychomotor skills. In more severe cases, hypotension, respiratory depression, profound sedation, coma, or even death may occur. **GENERALLY AVOID: Consumption of grapefruit juice during treatment with oral transmucosal formulations of fentanyl may result in increased plasma concentrations of fentanyl, which is primarily metabolized by CYP450 3A4 isoenzyme in the liver and intestine.** ***Certain compounds present in grapefruit are known to inhibit CYP450 3A4 and may increase the bioavailability of swallowed fentanyl (reportedly up to 75% of a dose) and/or decrease its systemic clearance.*** The clinical significance is unknown. In 12 healthy volunteers, consumption of 250 mL regular-strength grapefruit juice the night before and 100 mL double-strength grapefruit juice one hour before administration of oral transmucosal fentanyl citrate (600 or 800 mcg lozenge) did not significantly affect fentanyl pharmacokinetics, overall extent of fentanyl-induced miosis (miosis AUC), or subjective self-assessment of various clinical effects compared to control. However, pharmacokinetic alterations associated with interactions involving grapefruit juice are often subject to a high degree of interpatient variability. The possibility of significant interaction in some patients should be considered. MANAGEMENT: Patients should not consume alcoholic beverages or use drug products that contain alcohol during treatment with fentanyl. Any history of alcohol or illicit drug use should be considered when prescribing fentanyl, and therapy initiated at a lower dosage if necessary. Patients should be closely monitored for signs and symptoms of sedation, respiratory depression, and hypotension. **Due to a high degree of interpatient variability with respect to grapefruit juice interactions, patients treated with fentanyl should preferably avoid the consumption of grapefruit and grapefruit juice.** In addition, patients receiving transdermal formulations of fentanyl should be cautioned that drug interactions and drug effects may be observed for a prolonged period beyond removal of the patch, as significant amounts of fentanyl are absorbed from the skin for 17 hours or more after the patch is removed. |
Entonces, lo que tenemos que hacer es **quedarnos con los `texts_ID` únicos y ver su `link`**. El problema es que **todas las entradas** hacen referencia a las mismas frases del mismo texto, para una misma interacción...
:::info
Por ahora voy a quedarme con las **URL únicas** y a contar cuántas hay. Luego veremos. Miraré también las **combinaciones URL-texto únicas**.
:::
:::warning
**Guardar la salida en fichero de texto**
Para manejarme mejor he intentado guardar la salida en un fichero de texto, y me he topado con un error:
```sql=!
SELECT * FROM texts INTO OUTFILE '~/texts.txt';
```
```!
ERROR 1290 (HY000): The MySQL server is running with the --secure-file-priv option so it cannot execute this statement
```
Si miramos la variable `secure_file_priv`, veremos que los ficheros se guardarían en un directorio para el cual no tenemos acceso:
```sql=!
SHOW VARIABLES LIKE "secure_file_priv";
```
```
+------------------+-----------------------+
| Variable_name | Value |
+------------------+-----------------------+
| secure_file_priv | /var/lib/mysql-files/ |
+------------------+-----------------------+
1 row in set (0,01 sec)
```
Y para el cual no podemos cambiar los permisos:
```bash!
chmod u+rwx /var/lib/mysql-files/
```
```
chmod: changing permissions of '/var/lib/mysql-files/': Operation not permitted
```
Vamos a intentar solucionalo con [este tutorial](https://computingforgeeks.com/how-to-solve-mysql-server-is-running-with-the-secure-file-priv-error/). Voy a seguir la segunda opción:
```bash=!
sudo gedit /etc/mysql/my.cnf # no tengo vim instalado aun
```
>**Nota:** por si acaso he guardado una copia del fichero `my.cnf` en el home.
Añadimos las líneas:
```
[mysqld]
secure-file-priv = ""
```
Reiniciamos mysqld:
```bash!
sudo systemctl restart mysql
```
***
Solucionado! Guardamos la salida:
```sql
SELECT * FROM texts INTO OUTFILE '/tmp/texts.txt';
```
Le cambiamos los permisos:
```bash=!
sudo chmod ugo+rwx /tmp/texts.txt
```
:::
Y listo! Lo he copiado al home para manejarlo más fácilmente. s un fichero delimitado por tabuladores y con tres columnas: el ID del texto, el texto como tal y la URL. Vamos a quedarnos con las URL únicas:
```bash
cut -f3 texts.txt | sort | uniq | wc -l
```
Vemos que tenemos 2839 líneas --> solo 10 se repiten, y son:
```bash
cut -f3 texts.txt | sort | uniq -d
```
```
https://www.drugs.com/food-interactions/abemaciclib,verzenio.html?professional=1
https://www.drugs.com/food-interactions/fentanyl,sublimaze.html?professional=1
```
Y ya por fin:
```bash
cut -f3 texts.txt | sort | uniq | cut -d '/' -f3 | uniq -c
```
```
285 go.drugbank.com
2312 pubmed.ncbi.nlm.nih.gov
242 www.drugs.com
```
Si contamos con los que se repiten varias veces:
```bash
cut -f3 texts.txt | sort | cut -d '/' -f3 | uniq -c
```
```
285 go.drugbank.com
2312 pubmed.ncbi.nlm.nih.gov
252 www.drugs.com
```
### Número de perfiles GEO
Contando accession únicos:
:::spoiler
```
mysql> SELECT COUNT(DISTINCT accession) FROM study;
+---------------------------+
| COUNT(DISTINCT accession) |
+---------------------------+
| 161 |
+---------------------------+
1 row in set (0,00 sec)
```
La tabla tiene 161 filas, así que no se repite ningún estudio.
En cuanto al **"número total de experimentos"**, contamos con 2996 nodos. Cada nodo es una combinación **estudio-tratamiento-muestra-concentración**. Creo que Enrique quiere también ver cuántas muestras de distintos tratamientos tenemos. No sé si esto podría referirse a combinaciones **tratamiento-estudio** únicas.
Mirando en la tabla **`sample`**, tenemos 4873 combinaciones únicas estudio-tratamiento-muestra-concentración, y 3376 si excluimos los controles:
```sql!
SELECT COUNT(DISTINCT study_id,treatment,concentration,sample_id) FROM sample;
SELECT COUNT(DISTINCT study_id,treatment,concentration,sample_id) FROM sample WHERE treatment != 'control';
```
Si ahora nos vamos a la tabla **`node`**, veremos que tenemos 510 `node_id` distintos y 2996 `sample_id` distintos (la tabla tiene 2996 filas).
```sql
SELECT COUNT(DISTINCT node_id) FROM nodes;
SELECT COUNT(DISTINCT sample_id) FROM nodes;
```
:::
:::info
El número que le he dado a Enrique (**426 estudios**) sale de haber replicado la consulta Food-drug interactions con "resveratrol" usando MySQL y después haciendo una "búsqueda vacía":
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration, origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes NATURAL JOIN cmap_foodrugs
WHERE nodes.node_id = cmap_foodrugs.node_id; -- ambigüedad de atributos
-- 426 rows in set (0,02 sec)
```
Código completo:
:::spoiler
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration, origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes NATURAL JOIN cmap_foodrugs
WHERE nodes.node_id = cmap_foodrugs.node_id
INTO OUTFILE '/tmp/FinalFooDrugs_v2_experiments_no_header.txt';
-- Query OK, 426 rows affected (0,02 sec)
```
La movemos a un sitio donde pueda modificarla para incluir los headers:
```bash
# ejecutados en home
sudo cp /tmp/FinalFooDrugs_v2_experiments_no_header.txt ./FooDrugs/
sudo chmod ugo+rw FooDrugs/FinalFooDrugs_v2_experiments_no_header.txt
```
Y se los añado manualmente, separados por tabuladores
```
accession node_id treatment time_point concentration origin_name
```
```bash
# no me dejaba enviar el original por slack :/
cp experiments_FinalFooDrugs_v2.txt all_experiments_FinalFooDrugs_v2.txt
rm experiments_FinalFooDrugs_v2.txt
```
:::
Me ha costado un poco llegar hasta aquí. Algunas **notas:**
* El set de 401 entradas que obtuve al principio venía de que no incluía el `origin_name` en la selección, de forma que en estudios donde se ha usado el mismo tratamiento, durante el mismo tiempo y a la misma concentracón en líneas celulares diferentes solo veíamos una entrada.
* Siguiendo ese razonamiento, tendremos números diferentes quitando cualquiera de las otras categorías del `SELECT`.
* Si hacemos la query eliminando el `NATURAL JOIN cmap_foodrugs` salen 513 entradas, tenemos entonces 86 estudios para los que no hay nodos en CMAP:
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration
FROM study NATURAL JOIN sample NATURAL JOIN nodes;
-- 513 rows in set (0,02 sec)
```
#### Replicando una búsqueda en FooDrugs
Cuando entras en FooDrugs y buscas "resveratrol", sale:
La query que se hace por debajo es:
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration, origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes NATURAL JOIN cmap_foodrugs
WHERE nodes.node_id = cmap_foodrugs.node_id -- resuelve ambigüedad de atributos
AND treatment = 'resveratrol' AND (tau >= 90 OR tau <= -90);
```
```
+-----------+---------+-------------+------------+---------------+
| accession | node_id | treatment | time_point | concentration |
+-----------+---------+-------------+------------+---------------+
| 32357 | 350 | resveratrol | 30 days | 150 mg/day |
| 25412 | 351 | resveratrol | 48 hours | 150 mM |
| 25412 | 352 | resveratrol | 48 hours | 250 mM |
| 85871 | 353 | resveratrol | 12 hours | 10 µM |
| 42432 | 354 | resveratrol | 30 days | 150 mg/day |
| 78188 | 355 | resveratrol | 24 hours | - |
+-----------+---------+-------------+------------+---------------+
6 rows in set (0,03 sec)
```
:::success
Creo que lo que quiere Enrique es la salida de:
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration, origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes NATURAL JOIN cmap_foodrugs
WHERE nodes.node_id = cmap_foodrugs.node_id;
-- el WHERE es recomendable al tener ambigüedad de atributos
```
```
426 rows in set (0,02 sec)
```
Es decir, la misma query sin especificar ningún tratamiento ni restringir tau.
:::
:::info
Según el archivo de Marco:
```sql
SELECT DISTINCT study.accession, nodes.node_id, sample.treatment, sample.time_point, sample.concentration, sample.origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes
RIGHT JOIN cmap_foodrugs ON nodes.node_id = cmap_foodrugs.node_id
WHERE sample.treatment = ?; -- AND cell_line = ?
```
_"Cuando el usuario busca un fármaco con interacciones con compuestos alimenticios basados en evidencia molecular. Puedes añadir “`AND cell_line = ?`” si el usuario añade cell line como párametro de búsqueda también."_
:::
:::spoiler
Me queda la duda, sin embargo, porque si quitamos el join con `cmap_foodrugs` salen más entradas:
```sql=!
SELECT DISTINCT accession, node_id, treatment, time_point, concentration
FROM study NATURAL JOIN sample NATURAL JOIN nodes;
```
```
513 rows in set (0,02 sec)
```
Asumo que esto puede ser porque tengamos nodos para los cuales no haya información en el connectivity map.
Me guardo ambas para trastear:
```sql
SELECT DISTINCT accession, node_id, treatment, time_point, concentration
FROM study NATURAL JOIN sample NATURAL JOIN nodes
ORDER BY node_id
INTO OUTFILE '/tmp/njoin_513.txt';
SELECT DISTINCT accession, node_id, treatment, time_point, concentration
-- meter study_id en el SELECT no cambia el n de entradas recuperadas
FROM study NATURAL JOIN sample NATURAL JOIN nodes NATURAL JOIN cmap_foodrugs
ORDER BY node_id
INTO OUTFILE '/tmp/njoin_401.txt';
```
Si miramos la cabecera de ambos ficheros:
Vemos que en la query sin `cmap_foodrugs` aparecen los nodos 2, 4 y 5; y que en la otra no. De hecho, es un fichero de 513 líneas en el que el último `node_id` es el 519: solo faltan 6 nodos. En el otro, llegamos igualmente hasta el 519, pero solo tenemos 401 líneas: faltan 118 nodos. **Si buscamos estos nodos que faltan en `cmap_foodrugs`, vemos que no hay entradas con ese `node_id`:**
Mientras que, por ejemplo, para `node_id = 1` recuperamos 70895 filas.
#### Otra aproximación
Otra cosa que se me ha ocurrido es ver cuántos tratamientos diferentes tenemos por estudio:
```sql=!
SELECT study_id,treatment,COUNT(sample_id)
FROM sample
GROUP BY study_id,treatment;
-- sale lo mismo con SELECT DISTINCT
```
Con esto tenemos una tabla de 576 filas. Echemos un vistazo a las primeras líneas:
```
+----------+---------------------------+------------------+
| study_id | treatment | count(sample_id) |
+----------+---------------------------+------------------+
| 0 | cranberry extract | 3 |
| 0 | grapeseed extract | 3 |
| 0 | control | 3 |
| 0 | red wine extract | 3 |
| 1 | control | 18 |
| 1 | sulforaphane | 18 |
| 2 | NULL | 28 |
| 3 | control | 3 |
| 3 | orange extract | 3 |
| 4 | eusynstyelamide B | 3 |
| 4 | control | 6 |
```
Si quitamos los controles, nos quedan 421 filas; 420 si quitamos un estudio en el que `treatment = '-'`. El NULL que sale arriba no sé de dónde sale porque al buscar `WHERE treatment != 'control'` ya deja de salir. De hecho para este estudio solo tenemos las columnas `sample_id`, `study_id` y `GSM`.
No me termina de quedar claro qué diferencia hay entre esto y la "búsqueda vacía". Probando esto:
```sql=!
SELECT DISTINCT study_id,accession,treatment,COUNT(sample_id),COUNT(node_id)
FROM sample NATURAL JOIN study NATURAL JOIN nodes
GROUP BY study_id,treatment
ORDER BY accession;
```
Vemos que las columnas `COUNT(sample_id)` y `COUNT(node_id)` son iguales, eso sí, tenemos **418 filas**.
```
----------+-----------+-------------------------------------------------------------------------------------------+------------------+----------------+
| study_id | accession | treatment | COUNT(sample_id) | COUNT(node_id) |
+----------+-----------+-------------------------------------------------------------------------------------------+------------------+----------------+
| 69 | 5145 | vitamin D | 3 | 3 |
| 0 | 5556 | cranberry extract | 3 | 3 |
| 0 | 5556 | grapeseed extract | 3 | 3 |
| 0 | 5556 | red wine extract | 3 | 3 |
| 28 | 9647 | apple procyanidins extract | 3 | 3 |
| 28 | 9647 | apple procyanidins extract + tnf-α | 3 | 3 |
| 28 | 9647 | tnf-α | 3 | 3 |
| 76 | 11355 | Lactobacillus plantarum dead | 8 | 8 |
| 76 | 11355 | Lactobacillus plantarum midlog | 8 | 8 |
| 76 | 11355 | Lactobacillus plantarum stationary | 8 | 8 |
```
Comparo las dos queries:
```sql=!
SELECT DISTINCT study_id,accession,treatment,COUNT(sample_id),COUNT(node_id)
FROM sample NATURAL JOIN study NATURAL JOIN nodes
GROUP BY study_id,treatment
ORDER BY accession; -- 418 filas
SELECT DISTINCT accession, node_id, treatment, time_point, concentration
FROM study NATURAL JOIN sample NATURAL JOIN nodes
ORDER BY accession; -- 513 filas
```
Aquí veremos que hay diferencias porque en la segunda búsqueda diferenciamos también muestras con distintos tiempos o concentraciones, mientras que la primera lo agrupa todo junto. Podríamos sacar lo mismo si además agrupásemos por `time_point` y `concentration`:
```sql
SELECT DISTINCT study_id,accession,treatment,time_point,concentration,COUNT(sample_id)
FROM sample NATURAL JOIN study NATURAL JOIN nodes
GROUP BY study_id,treatment,time_point,concentration
ORDER BY accession; -- 513 filas
```
Y para replicar la búsqueda de FooDrugs:
```sql
SELECT study_id,accession,treatment,time_point,concentration,COUNT(sample_id),COUNT(node_id)
FROM sample NATURAL JOIN study NATURAL JOIN nodes
WHERE treatment = 'resveratrol'
GROUP BY study_id,treatment,time_point,concentration
ORDER BY accession;
```
:::success
En esta búsqueda **no tenemos controles** al estar haciendo el join con la tabla `nodes`.
:::
:::info
- [x] Enseñarle esto a Enrique y ver si es por aquí por donde van los tiros.
- [x] Preguntarle si puedo ir mirando algo más de esto.
:::
## Entendiendo las consultas
### Fármaco que interacciona con un compuesto alimenticio
Acabamos de verlo:
```sql
SELECT DISTINCT study.accession, nodes.node_id, sample.treatment, sample.time_point, sample.concentration, sample.origin_name
FROM study NATURAL JOIN sample NATURAL JOIN nodes
RIGHT JOIN cmap_foodrugs ON nodes.node_id = cmap_foodrugs.node_id
WHERE sample.treatment = ? -- AND cell_line = ?
```
* El símbolo `?` en `WHERE sample.treatment = ?` sirve como _placeholder_: ahí irá lo que introduzca el usuaio en la barra de búsqueda.
:::warning
- [ ] Queda resolver **por qué hace un `RIGHT JOIN`** en lugar de otro `NATURAL JOIN`
:::
### Elección de un alimento específico
Si tras la búsqueda que acabamos de hacer hacemos clic en la primera fila:
Estamos encontrando información para el nodo 350. Para replicar estos resultados necesitamos usar este comando:
```sql
SELECT node_id, cmap_node_id, tau, compound, cell_line
FROM cmap_foodrugs NATURAL JOIN cmap
WHERE node_id = 350 AND ABS(tau) >= 90 AND pert_type
= 'trt_cp'
ORDER BY ABS(tau) DESC
LIMIT 10;
```
```
+---------+--------------+----------+----------------+-----------+
| node_id | cmap_node_id | tau | compound | cell_line |
+---------+--------------+----------+----------------+-----------+
| 350 | 66114 | -99.9636 | goserelin | HT29 |
| 350 | 10272 | -99.9562 | proadifen | VCAP |
| 350 | 51780 | 99.9263 | eliprodil | HCC515 |
| 350 | 11398 | 99.8821 | iproniazid | HCC515 |
| 350 | 56906 | 99.882 | SB-216763 | A375 |
| 350 | 7296 | -99.8701 | BRD-K00256256 | A375 |
| 350 | 55115 | 99.8111 | SU-1498 | A375 |
| 350 | 29143 | 99.7939 | hydroxyfasudil | HT29 |
| 350 | 61490 | -99.7789 | CGP-54626 | HCC515 |
| 350 | 12281 | -99.6964 | esomeprazole | HT29 |
+---------+--------------+----------+----------------+-----------+
10 rows in set (0,20 sec)
```
:::success
Vemos que esto coincide con la query que venía en el documento:
```sql
SELECT DISTINCT cmap_foodrugs.*, cmap.compound, cmap.cell_line,
REPLACE(REPLACE(cmap.pert_type, 'trt_sh.cgs', 'Consensus signature FROM shRNAs targeting the same gene'), 'trt_oe', 'cDNA for overexpression of wild-type gene') as pert_type
FROM cmap_foodrugs NATURAL JOIN cmap
WHERE cmap_foodrugs.node_id = ? AND abs(tau) >= ?;
```
> "Cuando el usuario escoge un alimento especifico. A continuacion se añade `AND ( cmap.pert_type = 'trt_cp' or cmap.pert_type = 'trt_lig' )` si busca **ali-fármaco** o `AND (cmap.pert_type = 'trt_oe' or cmap.pert_type = 'trt_oe.mut' or cmap.pert_type = 'trt_sh' or cmap.pert_type = 'trt_sh.cgs' or cmap.pert_type = 'trt_sh.css' )` si busca **ali-gen**. Al final se añade `order by abs(cmap_foodrugs.tau) desc`"
:::spoiler
Para un ejemplo **ali-gen:**
```sql
SELECT DISTINCT cmap_foodrugs.*, cmap.compound, cmap.cell_line,
REPLACE(REPLACE(cmap.pert_type, 'trt_sh.cgs', 'Consensus signature FROM shRNAs targeting the same gene'), 'trt_oe', 'cDNA for overexpression of wild-type gene') as pert_type
FROM cmap_foodrugs NATURAL JOIN cmap
WHERE cmap_foodrugs.node_id = 350 AND abs(tau) >= 90 AND (cmap.pert_type = 'trt_oe' or cmap.pert_type = 'trt_oe.mut' or cmap.pert_type = 'trt_sh' or cmap.pert_type = 'trt_sh.cgs' or cmap.pert_type = 'trt_sh.css' )
ORDER BY ABS(tau) DESC
LIMIT 10;
```
```
+---------+--------------+----------+----------+-----------+---------------------------------------------------------+
| node_id | cmap_node_id | tau | compound | cell_line | pert_type |
+---------+--------------+----------+----------+-----------+---------------------------------------------------------+
| 350 | 46711 | 99.9728 | DCK | HEPG2 | cDNA for overexpression of wild-type gene |
| 350 | 61603 | 99.9562 | DDX10 | VCAP | Consensus signature FROM shRNAs targeting the same gene |
| 350 | 45768 | 99.9562 | ZNF385B | VCAP | Consensus signature FROM shRNAs targeting the same gene |
| 350 | 18934 | 99.9519 | IL18RAP | PC3 | cDNA for overexpression of wild-type gene |
| 350 | 2711 | 99.9514 | SLC7A1 | MCF7 | cDNA for overexpression of wild-type gene |
| 350 | 16116 | 99.9479 | FOXP3 | HT29 | cDNA for overexpression of wild-type gene |
| 350 | 37783 | 99.9415 | NGRN | HA1E | Consensus signature FROM shRNAs targeting the same gene |
| 350 | 43972 | 99.941 | MAP3K8 | HCC515 | Consensus signature FROM shRNAs targeting the same gene |
| 350 | 36718 | -99.941 | NRP1 | HCC515 | cDNA for overexpression of wild-type gene |
| 350 | 5673 | -99.9344 | BDKRB2 | VCAP | Consensus signature FROM shRNAs targeting the same gene |
+---------+--------------+----------+----------+-----------+---------------------------------------------------------+
10 rows in set (0,21 sec)
```
:::
### Compuesto alimenticio que interacciona con fármacos
Probamos con:
```sql
SELECT nodes.node_id, cmap.compound, cmap.cell_line, study.accession
FROM cmap NATURAL JOIN cmap_foodrugs NATURAL JOIN nodes NATURAL JOIN sample NATURAL JOIN study
WHERE cmap.cmap_node_id = cmap_foodrugs.cmap_node_id
AND cmap_foodrugs.node_id = nodes.node_id
AND nodes.sample_id = sample.sample_id
AND sample.study_id = study.study_id
AND cmap.compound = 'ibuprofen'
AND (tau >= 90 OR tau >= -90);
```
Esto tarda mucho. Si miramos el documento vemos que la solución es mucho más sencilla:
```sql
SELECT DISTINCT cmap_node_id AS node_id, compound, cell_line
FROM cmap
WHERE cmap.compound = ? AND (cmap.pert_type = 'trt_cp')
```
:::warning
* Aquí lo que no me queda claro es que este resultado realmente no te enseña con qué alimentos interacciona el fármaco...
* Además, en `node_id` en realidad vemos el `cmap_node_id`, no es un poco confuso? (para mí lo es)
:::
Para intentar resolver estas dudas voy a obtener los `node_id` correspondientes a los 8 `cmap_node_id` que tenemos en los resultados:
```sql
SELECT DISTINCT cmap_node_id, node_id, compound, cell_line
FROM cmap NATURAL JOIN cmap_foodrugs
WHERE cmap.compound = 'ibuprofen' AND cmap.pert_type = 'trt_cp';
```
:::warning
Y vemos que ahora nos salen **3184 filas** en vez de 8, porque cada nodo del CMAP se conecta con varios nodos nuestros (es decir, un solo compuesto del CMAP muestra interacciones con numerosos de los tratamientos de los estudios recogidos). Así, esta query nos dice **que este fármaco tiene interacciones con alimentos de los estudios recogidos**, pero no nos dice con cuáles (podrían ser muchos).
:::
:::info
**Columna `pert_type`**
Esta columna tiene 5 valores posibles:
```
mysql> select distinct pert_type from cmap;
+-------------+
| pert_type |
+-------------+
| trt_oe |
| trt_cp |
| trt_sh.cgs |
| ctl_vector |
| ctl_vehicle |
+-------------+
5 rows in set (0,02 sec)
```
En la [documentación de CMAP](https://clue.io/connectopedia/perturbagen_types_and_controls) podemos ver que `pert_type` viene de _**perturbagen type**_:
>"A perturbagen is a (chemical or genetic) reagent used in the laboratory to treat cells and measure the resulting biological response. In the case of CMap, changes in gene expression are measured following treatment with a perturbagen."
Hay más valores posibles en esta columna más allá de los 5 que tenemos nosotros. Centrándonos en estos 5 (el resto están en el link):
* **`trt_oe`**: cDNA for overexpression of wild-type gene
* **`trt_cp`**: Compound
* **`trt_sh.cgs`**: Consensus signature from shRNAs targeting the same gene
* **`ctl_vector`**: Controls - vector for genetic perturbation (e.g empty vector, GFP)
* **`ctl_vehicle`**: Controls - vehicle for compound treatment (e.g DMSO)
:::
### Elección de un nodo concreto
Si tras la búsqueda anterior seleccionamos la quinta fila:
Es lo mismo que esta query:
```sql
SELECT DISTINCT treatment, concentration, time_point, origin_name, cmap_foodrugs.*
FROM cmap_foodrugs NATURAL JOIN nodes NATURAL JOIN sample
WHERE cmap_node_id = '59203' AND ABS(tau) >= 90
ORDER BY ABS(tau) DESC;
```
```
+---------------------------------------+---------------+------------+-------------------------------+---------+--------------+----------+
| treatment | concentration | time_point | origin_name | node_id | cmap_node_id | tau |
+---------------------------------------+---------------+------------+-------------------------------+---------+--------------+----------+
| wortmannin | - | 24 hours | CUTLL1 | 512 | 59203 | 98.7589 |
| 3-o-acetyl-β-boswellic_acid | 46 μg/ml | 6-8 hours | MDA-MB-231 | 12 | 59203 | -97.4848 |
| tcdc | - | 11 days | Epithelium | 412 | 59203 | 97.3318 |
| phenoxibenzamine | - | 24 hours | CUTLL1 | 321 | 59203 | 96.2468 |
| after 2 weeks milk | - | - | Whole blood | 25 | 59203 | 94.7371 |
| H202 | 300µM | 24 hours | SK-N-MC | 208 | 59203 | -94.528 |
| boswellia serrata extract | 128 µg/mL | 6-8 hours | MDA-MB-231 | 76 | 59203 | -94.0978 |
| vitamin D | 100 nM | 24 hours | RWPE-1 | 480 | 59203 | 93.3946 |
| fruit and vegetable juice concentrate | - | 8 weeks | PBMC | 164 | 59203 | -93.126 |
| ferulic acid | - | 6 hours | breast cancer cell line MCF-7 | 147 | 59203 | 92.9195 |
+---------------------------------------+---------------+------------+-------------------------------+---------+--------------+----------+
10 rows in set (0,01 sec)
```
:::success
Lo que tenemos en el excel es:
```sql
SELECT DISTINCT sample.treatment, sample.concentration, sample.time_point, sample.origin_name, cmap_foodrugs.*
FROM study NATURAL JOIN sample NATURAL JOIN nodes RIGHT JOIN cmap_foodrugs ON nodes.node_id = cmap_foodrugs.node_id NATURAL JOIN cmap
WHERE cmap_foodrugs.cmap_node_id = ? AND abs(tau) >= ?
```
> Cuando el usuario escoge un farmaco especifico. A continuacion se añade `AND ( cmap.pert_type = 'trt_cp' or cmap.pert_type = 'trt_lig' )` y `order by abs(cmap_foodrugs.tau) desc`
:::warning
- [ ] Queda ver por qué se está hacendo el `RIGHT JOIN`
:::
### Text mining
Buscamos "resveratrol". Nos salen entradas con resveratrol tanto en la columna `food` como en `drug`, y con distintas mayúsculas / minúsculas:
:::warning
creo que aquí me quedé a medias
:::
## Entendiendo la tabla `nodes`
Vamos a empezar por las dos muestras a las que se les ha asignado el `node_id = 1`:
```sql
select * from nodes limit 3;
```
```
+---------+-----------+
| node_id | sample_id |
+---------+-----------+
| 1 | 4487 |
| 1 | 4503 |
| 2 | 1136 |
+---------+-----------+
3 rows in set (0,00 sec)
```
```sql
select * from sample where sample_id = 4487;
```
```
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
| sample_id | study_id | treatment | origin_type | origin_name | GSM | time_point | concentration |
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
| 4487 | 133 | prodrug g668684 | cell line | HCT-116 | 767904 | 72 hours | 150 nm |
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
1 row in set (0,01 sec)
```
```sql
select * from sample where sample_id = 4503;
```
```
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
| sample_id | study_id | treatment | origin_type | origin_name | GSM | time_point | concentration |
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
| 4503 | 133 | prodrug g668684 | cell line | HCT-116 | 767903 | 72 hours | 150 nm |
+-----------+----------+------------------+-------------+-------------+--------+------------+---------------+
1 row in set (0,00 sec)
```
Vemos que solo cambian las columnas `GSM` y `sample_id`: los **identificadores de las muestras**.
Probando con `node_id = 2` sacamos 160 entradas, que se quedan en 80 si hacemos:
```sql
select distinct study_id, treatment, origin_type, origin_name, GSM, time_point, concentration from sample natural join nodes where node_id =2;
```
Si quitamos `GSM` de la query anterior, obtenemos una sola fila.
Si ahora hacemos esto:
```sql
select distinct node_id, study_id, treatment, origin_type, origin_name, time_point, concentration from sample natural join nodes;
```
Recuperamos **571** entradas, de manera que hay `node_id` repetidos. Me llevo esto a un fichero para echar un vistazo.
:::warning
El `node_id = 74` no existe, el `465` tampoco... y así con más
:::
### A partir de excel
Me he quedado con la tabla que sale de la última query (combinaciones únicas de ID de nodo, ID de estudio, tratamiento, tipo de origen y nombre, tiempo y concentración) y me he dedicado a buscar **muestras distintas con el mismo ID de nodo**. Tenemos 510 `node_id` distintos y 2996 `sample_id` distintos en la tabla `nodes`, así que algo no cuadra.
:::info
Aquí ya estamos filtrando mucho, porque puede haber duplicados por ejemplo que no vamos a ver al tener todo igual salvo por los identificadores de muestra (`sample_id` y `GSM`)
:::
Aquí vemos que hay filas repetidas, por ejemplo, para muestras donde lo único que cambia es el `origin_name`, a veces junto con `origin_type`; y unas pocas donde lo único que cambia es el tiempo.
:::danger
De esto lo que me raya es que se supone que los nodos son combinaciones **tiempo, concentración y tratamiento** únicas, y estamos viendo que hay veces donde dos muestras con distintos tiempos se han asignado al mismo nodo... Nodos 265, 302 y 427.
Esto **no es un problema de las muestras con `time_point = 0`.** Por ejemplo, los nodos 357 y 203 (estudio 52) únicamente tienen muestras con ese tiempo, sin mezclar con otras muestras del estudio. El nodo 357 solo tiene muestras de t = 0, y en ese mismo estudio se miran también 6 meses y 12 meses. El nodo 203 es lo mismo, porque es el mismo estudio solo que con otro compuesto.
Ejemplo:
:::spoiler
```
mysql> select * from nodes natural join sample where node_id = 265;
+-----------+---------+----------+--------------+-------------+-------------+--------+------------+---------------+
| sample_id | node_id | study_id | treatment | origin_type | origin_name | GSM | time_point | concentration |
+-----------+---------+----------+--------------+-------------+-------------+--------+------------+---------------+
| 354 | 265 | 15 | low fat diet | tissue | PBMC | 701130 | 3 months | - |
| 357 | 265 | 15 | low fat diet | tissue | PBMC | 701133 | 0 | - |
| 358 | 265 | 15 | low fat diet | tissue | PBMC | 701136 | 3 months | - |
| 360 | 265 | 15 | low fat diet | tissue | PBMC | 701139 | 0 | - |
| 363 | 265 | 15 | low fat diet | tissue | PBMC | 701142 | 0 | - |
| 364 | 265 | 15 | low fat diet | tissue | PBMC | 701145 | 3 months | - |
| 369 | 265 | 15 | low fat diet | tissue | PBMC | 701135 | 0 | - |
| 374 | 265 | 15 | low fat diet | tissue | PBMC | 701138 | 3 months | - |
| 378 | 265 | 15 | low fat diet | tissue | PBMC | 701141 | 0 | - |
| 381 | 265 | 15 | low fat diet | tissue | PBMC | 701144 | 0 | - |
| 382 | 265 | 15 | low fat diet | tissue | PBMC | 701147 | 3 months | - |
| 385 | 265 | 15 | low fat diet | tissue | PBMC | 701126 | 0 | - |
| 386 | 265 | 15 | low fat diet | tissue | PBMC | 701127 | 3 months | - |
| 391 | 265 | 15 | low fat diet | tissue | PBMC | 701143 | 3 months | - |
| 393 | 265 | 15 | low fat diet | tissue | PBMC | 701129 | 3 months | - |
| 395 | 265 | 15 | low fat diet | tissue | PBMC | 701146 | 0 | - |
| 399 | 265 | 15 | low fat diet | tissue | PBMC | 701152 | 3 months | - |
| 410 | 265 | 15 | low fat diet | tissue | PBMC | 701140 | 3 months | - |
| 413 | 265 | 15 | low fat diet | tissue | PBMC | 701128 | 3 months | - |
| 414 | 265 | 15 | low fat diet | tissue | PBMC | 701134 | 3 months | - |
| 415 | 265 | 15 | low fat diet | tissue | PBMC | 701137 | 0 | - |
+-----------+---------+----------+--------------+-------------+-------------+--------+------------+---------------+
21 rows in set (0,00 sec)
```
:::
También sobra una fila correspondiente al estudio 27, donde en lugar de control pone `'ref'`. Lo mismo con el 40.
Los nodos se han sacado de una forma muy parecida a esto:
```sql
select distinct study_id, time_point, concentration, treatment
from sample
where treatment != 'control';
-- 530 entradas, sobran 20
select distinct study_id, time_point, concentration, treatment
from sample
where treatment != 'control'
into outfile '/tmp/distinct_study_tp_conc_ttmt.txt';
select distinct node_id, study_id, treatment, time_point, concentration
from sample natural join nodes
into outfile '/tmp/nodes_513.txt';
-- Query OK, 513 rows affected (0,01 sec)
```
:::success
Sacamos estas dos tablas, miramos cuántas veces se repite cada `study_id` y miramos aquellas filas donde no coincida. Si sumamos las diferencias da 17. La cosa está en que en la tabla que hemos hecho nosotros manualmente sin mirar `nodes` hay 17 filas que en la otra no están. Vamos a ver qué pasa... :male-detective:
:::
```r
setwd('~/FooDrugs')
library(dplyr)
# nodes.tab contains the result of a natural join operation between the
# sample and the nodes tables in FooDrugs database:
nodes.tab <- read.table('nodes_513.txt', header = FALSE, sep = '\t',
na.strings = '-')
colnames(nodes.tab) <- c('node_id', 'study_id', 'treatment',
'time_point', 'concentration')
# distinct.tab contains the result of a select operation within the
# sample table selecting unique study-time point-concentration-treatment
# combinations, **excluding controls**
distinct.tab <- read.table('distinct_study_tp_conc_ttmt.txt', header = FALSE,
sep = '\t', na.strings = '-')
colnames(distinct.tab) <- c('study_id', 'time_point', 'concentration',
'treatment')
# now we will compare the number of times that each study_id appears
# in both tables:
nodes.times <- nodes.tab %>% count(study_id)
distinct.times <- distinct.tab %>% count(study_id)
times <- nodes.times %>% left_join(distinct.times, by = 'study_id')
colnames(times) <- c('study_id', 'nodes', 'distinct')
times$equal <- times['nodes'] == times['distinct']
# see difference:
tmp <- times %>% filter(equal == FALSE)
tmp[distinct] - tmp[nodes]
tmp$distinct - tmp$nodes
sum(tmp$distinct - tmp$nodes)
# we can see that the difference equals 17, mistery solved!
```
```
study_id nodes distinct nodes
1 17 4 7 FALSE
2 27 1 2 FALSE
3 40 1 2 FALSE
4 80 3 4 FALSE
5 96 1 2 FALSE
6 103 2 3 FALSE
7 110 1 2 FALSE
8 113 1 2 FALSE
9 115 1 2 FALSE
10 117 6 9 FALSE
11 141 3 6 FALSE
```
- **17**: faltan 3 con `time_point = 0`
- **27**: la que falta yo creo que es un control, solo que pone `treatment = ref` en lugar de control
- **40**: `treatment = ref`
- **80**: falta fila com `time_point = 0`
- **96**: falta fila com `time_point = 0`
- **103**: `treatment = ref`
- **110**: falta fila com `time_point = 0`
- **113**: falta fila com `time_point = 0`
- **115**: falta fila com `time_point = 0`
- **117**: faltan 3 con `time_point = 0`
- **141**: faltan 3 con `time_point = 0`
:::danger
Al buscar:
```r
nodes.tab %>% filter(time_point == 0)
```
Salen 6 filas:
```
node_id study_id treatment time_point concentration
1 203 52 grape extract 0 8 mg
2 265 15 low fat diet 0 <NA>
3 302 15 nuts 0 <NA>
4 357 52 resveratrol-enriched grape extract 0 8 mg
5 427 15 Traditional Mediterranean Diet enriched with virgin olive oil 0 <NA>
```
No sé si las que fallan son estas o las que hemos visto antes: **¿los nodos con `time_point = 0` deben existir o no?**
Además, los nodos 265. 302 y 427 tienen dos combinaciones distintas asignadas, porque se les ha asignado la muestra com `time_point = 3 months` y también la muestra con `time_point = 0`
:::
## Tabla Enrique
### Duda condiciones
Enrique me ha pedido que considere que cada condición incluya también los `origin_name` distintos, excluyendo controles. Entonces ahora nos surge una nueva duda: al hacer los "nodos" manualmente nos salen más filas que al hacer el natural join con nodes:
```sql
-- hacemos los nodos manualmente: 585 filas
SELECT DISTINCT study_id, treatment, time_point, concentration, origin_name, origin_type
FROM study NATURAL JOIN sample
WHERE treatment != 'control' and treatment != 'ref'
INTO OUTFILE '/tmp/origin_585.txt';
-- 585 filas
```
```sql
-- natural join con nodes: 571 filas
SELECT DISTINCT study_id, node_id, treatment, time_point, concentration, origin_name, origin_type
FROM study NATURAL JOIN sample NATURAL JOIN nodes
INTO OUTFILE '/tmp/origin_571.txt';
```
Vamos a hacer lo mismo de mirar en R los sitios donde haya diferencias.
Primero sacamos las tablas:
```
mysql> SELECT DISTINCT study_id, node_id, treatment, time_point, concentration, origin_name, origin_type
-> FROM study NATURAL JOIN sample NATURAL JOIN nodes
-> INTO OUTFILE '/tmp/origin_571.txt';
Query OK, 571 rows affected (0,01 sec)
mysql> SELECT DISTINCT study_id, treatment, time_point, concentration, origin_name, origin_type
-> FROM study NATURAL JOIN sample
-> WHERE treatment != 'control' and treatment != 'ref'
-> INTO OUTFILE '/tmp/origin_585.txt';
Query OK, 585 rows affected (0,04 sec)
```
Y corremos el script exactamente igual. Sale esta tabla:
```
study_id nodes distinct nodes
1 17 4 7 FALSE
2 80 3 4 FALSE
3 96 1 2 FALSE
4 110 1 2 FALSE
5 113 1 2 FALSE
6 115 1 2 FALSE
7 117 6 9 FALSE
8 141 3 6 FALSE
```
:::success
Vemos que son los mismos estudios que antes, quitando los de `treatment = 'ref'`. Es decir, **nos faltan nodos con `time_point = 0`**, que según lo hablado con Enrique sí deben existir.
**Por ello para el report voy a trabajar con 585 condiciones.**
:::
### Texto
A total of ~~119~~ **161** GEO series with ~~334~~ **585** conditions for food or bioactive compounds were retrieved following the strategy described above. **Each condition is defined as a food/biocomponent per time point per concentration per cell line (or tissue, or primary cell) per study.** Characteristics of these series are shown in Table XXX.
:::info
El **426** sale de una búsqueda vacía hecha según se veía [aquí](https://hackmd.io/zw-eHRfgT-GWhNebtAh7AA?view#N%C3%BAmero-de-perfiles-GEO).
El **510** sale del número de nodos (`node_id` únicos).
:::
One-color arrays
~~92~~ **131**
Two-color arrays
<font color = 'red'> ~~15~~ **10** </font>
RNA seq
~~12~~ **20**
:::info
Esto sale de:
```sql
select study_type,count(distinct study_id) from study group by study_type;
```
```
+-------------------+--------------------------+
| study_type | count(distinct study_id) |
+-------------------+--------------------------+
| HTS | 20 |
| One color array | 131 |
| Two channel array | 10 |
+-------------------+--------------------------+
3 rows in set (0,00 sec)
```
:::danger
Esto es un poco raro. Perdemos 5 estudios de arrays de dos colores...
- [ ] Buscar GEO accession en TFM / Preguntar a Enrique
- [ ] Si no... Se me ocurre mirar si en la versión antigua que me pasó David pasa lo mismo
:::
Studies conducted on cell lines
~~63~~ **96**
Studies conducted on primary cells
~~66~~ **27**
Studies conducted on biopsies
~~10~~ **41**
:::info
Si hacemos lo mismo que acabamos de ver pero con `origin_type` vemos que vamos a tener que mirar alguna cosa manualmente:
```sql
select origin_type, count(distinct study_id)
from sample
group by origin_type;
```
```
+-----------------------------------------------+--------------------------+
| origin_type | count(distinct study_id) |
+-----------------------------------------------+--------------------------+
| NULL | 1 |
| generated from H9 human embryonic stem cells | 1 |
| vehicle treated | 1 |
| vorinostat treated | 1 |
| - | 3 |
| CCD-18Co | 1 |
| cel type | 1 |
| cell line | 81 |
| cell type | 34 |
| regulatory T cell | 1 |
| tissue | 40 |
+-----------------------------------------------+--------------------------+
11 rows in set (0,02 sec)
```
Hemos ido mirando los estudios que no coincidían con `cell_line`, `cell_type` ni con `tissue`. Ahora queda mirar estas tres categorías:
- en `cell_line` solo hay líneas celulares;
- en `cell_type` hay línas celulares y células primarias;
- en `tissue` hay biopsias y células primarias.
Voy a sacar una tablita para ir haciendo cálculos:
```sql
select distinct origin_name, accession
from sample natural join study
where origin_type = 'tissue' or origin_type = 'cell type'
order by origin_name
into outfile '/tmp/cell_type_or_tissue.txt';
-- Query OK, 84 rows affected (0,01 sec)
```
Y ahora vamos a ir eliminando cosas.
1. Sabemos que todo lo que tenga `primary` son **células primarias:**
```bash
grep -i 'primary' cell_type_or_tissue.txt # 5 todos distintos
grep -iv 'primary' cell_type_or_tissue.txt > ctt2.txt
```
Vemos que hay mezcla de `cell type` y `tissue` en **`origin_type`**:
```
(base) blanca@totoro:~/FooDrugs$ grep 'primary' ctt_3cols.txt
cell type primary chondrocytes 75181
tissue primary human bronchial smooth muscle cells 5145
cell type primary myotubes 18589
tissue primary tumor sample 110248
```
2. Todo lo que contenga 'blood' o 'PBMC' viene de una **biopsia líquida**:
```bash
grep -iE 'blood|PBMC' ctt2.txt # funciona bien, todos unicos
grep -iEc 'blood|PBMC' ctt2.txt # 20 lineas
grep -iEv 'blood|PBMC' ctt2.txt > ctt3.txt
```
Aquí también tenemos mezcla.
3. Hay algunas cosas que sabemos que son **cell lines**:
- ARPE-19 Cells
- BEAS-2B immortalized human bronchial epithelium
- CACO-2
- HUVEC
- LNCaP
- MCF7
- PC-3
- Raji cells
- SUDHL10 Cells
- SUDHL2 Cells
```bash
grep -iEc 'arpe|beas-2b|caco-2|huvec|lncap|mcf7|pc-3|raji|sudhl' ctt3.txt # 10 lineas, 7 unicos
grep -iEv 'arpe|beas-2b|caco-2|huvec|lncap|mcf7|pc-3|raji|sudhl' ctt3.txt > ctt4.txt
```
Ya solo nos queda discriminar **biopsia _vs._ primarias**.
4. En `classification_cl_prim_biopsy.ods` pude verse el resto de la clasificación (hecha manualmente)
Otro:
- PrEC **creo** que son **primarias** cotilleando el [fabricante](https://bioscience.lonza.com/lonza_bs/CH/en/Primary-and-Stem-Cells/p/000000000000185069/PrEC-%E2%80%93-Human-Prostate-Epithelial-Cells)
- RCC PDX tumor son **primarias** pero derivadas de una **biopsia**... ([link](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1973613)) --> para mí esto cuenta como primarias
Para confimar he pasado todo a ficheros independientes (uno para cada cosa clasificada) y hacemos el `sort` + `uniq`.
```bash=
grep -i 'primary' ctt_3cols.txt > primary.txt
grep -iE 'arpe|beas-2b|caco-2|huvec|lncap|mcf7|pc-3|raji|sudhl' ctt_3cols.txt > cell_line.txt
grep -iE 'blood|PBMC' ctt_3cols.txt > biopsy.txt
# pego manualmente las tablas del excel
# pego manualmente las del cuaderno
```
:::danger
- He decidido que si en algún sitio pone que se aisló algún tipo celular concreto es **primary** y si no, es **biopsy**
- Muestras de leche y heces son **biopsy**
- tejido 3D y células reprogramadas son **cell line**
:::
Number of different GEO platforms
~~48~~ **69**
:::info
Aquí hacemos:
```sql
select distinct gpl from study;
```
Y vemos que con esto no vale, porque hay estudios con más de una plataforma. Pongo las primeras 5 entradas que se recuperan:
```
+-------------------------------+
| gpl |
+-------------------------------+
| 570 |
| 11154 |
| 6104 |
| 16604 |
| 1;7;6;9;2 |
```
Toca entonces volver a sacar un fichero de texto.
```sql
select distinct gpl from study into outfile '/tmp/gpl.txt';
```
Ahora es muy sencillo calcular el número de valores diferentes en bash:
```bash
cat gpl.txt | tr ';' '\n' | sort | uniq | wc -l
# salida: 69
```
:::
Series which study more than one compound
~~34~~ **66**
:::info
Aquí me salen bastantes más, lo he hecho de esta forma:
```sql
select study_id, count(distinct treatment)
from sample
where treatment != 'control' and treatment != 'ref'
group by study_id
having count(distinct treatment) > 1;
-- 66 rows in set
```
Probamos otro approach por confirmar:
```sql
select distinct study_id, treatment
from sample
where treatment != 'control' and treatment != 'ref'
into outfile '/tmp/tmp.txt';
```
```bash
cut -f1 tmp.txt | uniq -d | wc -l
# 66
```
:::
Series which study more than one concentration
~~5~~ **39**
:::info
```sql
select study_id, count(distinct concentration)
from sample
where treatment != 'control' and treatment != 'ref'
group by study_id
having count(distinct concentration) > 1;
-- 39 rows
```
:::
Series which study more than one time point
~~14~~ **33**
:::info
```sql
select study_id, count(distinct time_point)
from sample
where treatment != 'control' and treatment != 'ref'
group by study_id
having count(distinct time_point) > 1;
```
:::
:::success
Todo esto lo he hecho también en `compare_tables.R` y los números coinciden :confetti_ball:
```r=
# more than 1 treatment
by.treatment <- cl.distinct.tab %>%
select(study_id, treatment) %>%
distinct() %>%
count(study_id)
by.treatment %>% filter(n > 1) # 66 rows --> OK
# more than 1 concentration
by.concentration <- cl.distinct.tab %>%
select(study_id, concentration) %>%
distinct() %>%
count(study_id)
by.concentration %>% filter(n > 1) # 39 rows --> OK
# more than 1 time point
by.time_point <- cl.distinct.tab %>%
select(study_id, time_point) %>%
distinct() %>%
count(study_id)
by.time_point %>% filter(n > 1) # 33 rows --> OK
```
:::
Table XXX.- Characteristics of GEO series retrieved
## Creación de tablas - trigger en tabla `sample`
:::spoiler
```sql
--
-- Table structure for table `sample`
--
DROP TABLE IF EXISTS `sample`;
/*!40101 SET @saved_cs_client = @@character_set_client */;
/*!50503 SET character_set_client = utf8mb4 */;
CREATE TABLE `sample` (
`sample_id` int NOT NULL AUTO_INCREMENT,
`study_id` int NOT NULL,
`treatment` varchar(163) DEFAULT 'control',
`origin_type` varchar(50) DEFAULT 'cell line',
`origin_name` varchar(150) DEFAULT NULL,
`GSM` int NOT NULL,
`time_point` varchar(50) DEFAULT NULL,
`concentration` varchar(150) DEFAULT '0',
PRIMARY KEY (`sample_id`),
KEY `study_id` (`study_id`),
CONSTRAINT `sample_ibfk_1` FOREIGN KEY (`study_id`) REFERENCES `study` (`study_id`) ON DELETE CASCADE
) ENGINE=InnoDB AUTO_INCREMENT=5069 DEFAULT CHARSET=utf8mb4 COLLATE=utf8mb4_0900_ai_ci;
/*!40101 SET character_set_client = @saved_cs_client */;
--
-- Dumping data for table `sample`
--
LOCK TABLES `sample` WRITE;
/*!40000 ALTER TABLE `sample` DISABLE KEYS */;
INSERT INTO `sample` VALUES (0,0,'cranberry extract','cell type','endothelial cells',129143,'6 hours','-'), --...
/*!40000 ALTER TABLE `sample` ENABLE KEYS */;
UNLOCK TABLES;
/*!50003 SET @saved_cs_client = @@character_set_client */ ;
/*!50003 SET @saved_cs_results = @@character_set_results */ ;
/*!50003 SET @saved_col_connection = @@collation_connection */ ;
/*!50003 SET character_set_client = utf8mb4 */ ;
/*!50003 SET character_set_results = utf8mb4 */ ;
/*!50003 SET collation_connection = utf8mb4_0900_ai_ci */ ;
/*!50003 SET @saved_sql_mode = @@sql_mode */ ;
/*!50003 SET sql_mode = 'ONLY_FULL_GROUP_BY,STRICT_TRANS_TABLES,NO_ZERO_IN_DATE,NO_ZERO_DATE,ERROR_FOR_DIVISION_BY_ZERO,NO_ENGINE_SUBSTITUTION' */ ;
DELIMITER ;;
/*!50003 CREATE*/ /*!50017 DEFINER=`marco2`@`localhost`*/ /*!50003 TRIGGER `create_node` AFTER INSERT ON `sample` FOR EACH ROW begin
declare id_to_insert integer;
declare is_occupied integer;
case
when (NEW.concentration != "0" and NEW.time_point != "0") and (NEW.treatment != "control")
then
set is_occupied :=(SELECT EXISTS (SELECT 1 FROM `nodes`));
case
when is_occupied = 1
then
set id_to_insert := (select COALESCE( ( select distinct node_id from sample natural join nodes
where treatment = NEW.treatment and time_point = NEW.time_point and concentration = NEW.concentration LIMIT 1),
(select max(node_id) from nodes) + 1));
else
set id_to_insert=0;
end case;
insert into nodes(node_id, sample_id)
values( id_to_insert, NEW.sample_id);
else begin end;
end case;
end */;;
DELIMITER ;
/*!50003 SET sql_mode = @saved_sql_mode */ ;
/*!50003 SET character_set_client = @saved_cs_client */ ;
/*!50003 SET character_set_results = @saved_cs_results */ ;
/*!50003 SET collation_connection = @saved_col_connection */ ;
```
:::
# Pipeline de transcriptómica
## Lista de scripts
1. A partir de una lista de GSEs, se corre el pipeline de análisis diferencial desde el script **`runLimmaFromGSEList.sh`** para estudios de microarrays o desde el script **`runLimmaHTSFromGSEList.sh`** para estudios de transcriptómica.
2. Estos scripts a su vez llaman a otros scripts:
- `create_DM.py` para generar la matriz de diseño y matriz de nodos
- `DownloadDataMatrix.sh` para descargar datos CEL en el caso de los estudios de microarrays (en transciptómica asume `HTS_gene_counts/numerogse_counts.tsv`)
- **`runLimma.R`** para correr el análisis y obtener las topTables y los ficheros GMT
:::info
- Los ficheros **GMT (*Gene Matrix Transposed*)** recogen los nombres de los genes asociados a cada gene set. En este caso lo que entiendo que tenemos es un fichero por gene set (UP/DN), con una sola fila que dice los nombres de los genes
- Las **topTables** son tablas rankeadas con el top de genes que salen al aplicar el modelo lineal a un estudio.
:::
3. El script **`runLimma.R`** requiere:
- **Tres argumentos:** (1) una matriz de diseño, (2) un archivo de expresión y (3) un archivo de nodos;
- **Dos directorios:** (1) gmtFiles, donde escribir los gene sets para CMAP, y (2) topTable, para escribir los topTables.
4. Además, **`runLimma.R`** recoge funciones de otros tres scripts:
- `detectDesign.R`: funciones para obtener información sobre cada estudio y generar la matriz de contraste y de diseño que hay que pasarle a `limma`
- `HTS_functions.R`: funciones para analizar RNASeq
- `normalizedArrayFunctions.R`: funciones para el análisis de microarrays. Tiene unas cuantas funciones:
- varias para pasar identificadores de genes a ENSEMBL;
- **correr el modelo lineal, hacer la corrección de eBayes y sacar los contrastes;**
- **escribir los topTables;**
- **detectar top 150 y bottom 150 genes del BINGSpace para la carpeta gmtFiles;**
- ver qué hacer si un probeset vincula a varios genes o viceversa.
## Proceso `runLimma.R`
Es un script bastante corto. Podemos quedarnos con que:
1. Incorpora las funciones necesarias:
```r=
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 8)
source("detectDesign.R")
source("normalizedArraysFunctions.R")
source("HTS_functions.R")
```
2. Parsea el input y prepara el output:
```r=+
args = commandArgs(trailingOnly=TRUE)
DM <- args[1]
data_matrix <- args[2]
nodes <- args[3]
#filename suffix for outputs
f <- strsplit(data_matrix, "/")[[1]]
f <- f[length(f)]
f <- strsplit(f, "_")[[1]][1]
```
3. Prepara la matriz de diseño y detecta si el estudio es HTS o microarrays:
```r=+
con <- file(DM,"r")
first_line <- readLines(con,n=1)
ret <- detect_design(DM, log = "", sep = "\t")
#print(ret$designMatrix)
print(ret$contMatrix)
```
4. En función de si es HTS o microarrays, prepara los datos necesarios para el resto del análisis:
```r=+
#checker wheter HTS or microarrays
if(first_line == "#HTS"){
counts <- prepare_matrix(data_matrix, sep = "\t", header = 1)
annotatedEset <- proccessHTS(counts, ret$designMatrix)
}else{
annotatedEset <- annotateGSE(data_matrix)
}
```
5. Prepara la matriz de nodos:
```r=+
node_ids <- read.table(nodes, sep = "\t", header = T)
```
6. Corre el análisis diferencial (modelos lineales, etc.):
```r=+
generalizedTtest(annotatedEset, ret$designMatrix, ret$contMatrix,node_ids,"",".", f)
```
:::info
Todos los pasos de correr el análisis, quedarse con el top de genes y preparar las cosas para llamar a CMAP están dentro de la función `generalizedTtest`, que está definida dentro del script `normalizedArrayFunctions.R`
:::
### Función `generalizedTtest`
La función toma como argumentos:
- `annotatedGseset`: datos de expresión a utilizar
- `designMatrix`: matriz de diseño
- `contMatrix`: matriz de contrastes
- `id_matrix`: matriz de nodos
- `log = ""`: para fichero log cuando se busquen genes del BINGspace (ahí veremos si se encontraon menos de 0 significativos, qué ficheros UP/DN se han grabado, etc.)
- `direc`: directorio donde se guardará la salida (archivos gmt y topTables)
- `file`: nombre de la GSE para generar después los nombres de los archivos de salida
```r=
generalizedTtest<- function(annotatedGseset, designMatrix, contMatrix,
id_matrix, log = "", direc, file ){
designMatrix <- designMatrix[order(match(paste("GSM",
rownames(designMatrix),
sep = ""),
colnames(annotatedGseset))), ,
drop = F]
```
Transforma los datos de expresión logarítmicamente:
```r=+
#reorder the design matrix
#so row ordering matches column ordering of annotatedGseset
#limma will not look for rownames
exprs(annotatedGseset) <- log2(exprs(annotatedGseset))
```
Ajusta el modelo lineal ( `lmFit`), lleva a cabo los contrastes (`contrasts.Fit`) y aplica la corrección de eBayes (`eBayes`):
:::info
Manual de `eBayes`:
> Given a linear model fit from `lmFit`, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a global value.
> **Aquí no hace ninguna corrección de p-valores**.
:::
```r=+
fit <- lmFit(annotatedGseset, designMatrix)
fit <- contrasts.fit(fit, contMatrix)
fit <- eBayes(fit)
```
Para cada contraste de interés (cada `coeff` o columna de la matriz de contraste):
```r=+
for(coeff in colnames(contMatrix)){
GSMs <- names(which(designMatrix[,coeff] != 0))
node <- id_matrix[ id_matrix$GSM %in% strtoi(GSMs), "node_id" ][1]
```
- Obtiene el top de genes con `topTable`, aportándole:
- el modelo lineal ajustado (`fit`);
- el contraste que nos interesa (`coeff`);
- el máximo de genes de la lista (`nrow(annotatedGseset)`);
- el método para ajustar el p-valor (`"BH"`);
- el estadístico usado para hacer el ranking de genes (`"logFC"`)
:::danger
Mi problema aquí es que **no se ha dado un valor al parámetro `p.value`**, que indica el **punto de corte de los p-valores tras el ajuste**.
**El valor por defecto de este parámetro es `p.value = 1`.**
:::
```r=+
tt <- topTable(fit,
coef = coeff,
number= nrow(annotatedGseset),
adjust.method="BH",
sort.by="logFC")
```
Aplica la función `getMaxMean`, definida también en el script `normalizedArrayFunctions.R`, para quedarse con una sola fila si hay varias que mapeen al mismo gen:
```r=+
summarizedGenes <- getMaxMean(tt[ , c("entrezIDs", "logFC"), drop = F])
```
Escribe la topTable en el fichero de salida:
```r=+
write.table(tt,
file = paste( direc,"/topTable/", node, "___",
file, "_top_table_", coeff ,
sep = "") )
```
Llama a la función `getBINGSpaceGenes`, definida en el propio script, para escribir finalmente los ficheros gmt que se enviarán a CMAP:
:::success
Esta es la función que obtiene el top 150 y botttom 150 de genes y que escribe los ficheros a CMAP.
:::
```r=+
getBINGSpaceGenes(signGenes = summarizedGenes,
probeToEntrezIDs = F,
name = paste(node, "_",file,coeff, sep = "_"),
log = log,
direc = direc )
}
}
```
## Función `getBINGSpaceGenes`
Lo que interpreto que hace esta función es:
1. Quedarse con todo el conjunto de genes up / down de la topTable (que está en `signGenes`, generando los dataframes `upgenes` y `dngenes`:
```r=
getBINGSpaceGenes <- function(signGenes, col = "entrezIDs",probeToEntrezIDs = F,topGenes = 150, name, log = "", direc = "."){
#probeToEntrezIDs: dataframe of rows:probeset IDs, one column: entrez gene IDs
#signGenes: dataframe of rows: probeset IDs, columns logFC
#generate gmt fileto use in clue io
#clue io uses only a subset of genespace, as well as a limit of 150 genes
#subset those allowed, and if there are more than 150, sort based on logFC
if(probeToEntrezIDs){ ## OJO: esto por default es FALSE
upgenes <- merge(subset(signGenes, signGenes$logFC > 0), probeToEntrezIDs , 0)
rownames(upgenes) <- upgenes$Row.names
upgenes$Row.names <- NULL
dngenes <- merge(subset(signGenes, signGenes$logFC < 0), gseset@featureData@data, 0)
rownames(dngenes) <- dngenes$Row.names
dngenes$Row.names <- NULL
}
else{ ## lo que corre por defecto
## col es el entrezID
signGenes[, col] <- as.integer(signGenes[, col])
upgenes <- subset(signGenes, signGenes$logFC > 0)
dngenes <- subset(signGenes, signGenes$logFC < 0)
}
```
2. Selecciona los genes up que están en el BINGSpace (el espacio de genes recogidos en CMAP), si son más de 150, los ordena por logFC y toma el top 150. Si son menos de 10, manda un aviso y **para el análisis**, de manera que **no se hace ninguna query a CMAP**:
:::danger
Como en `topTable` no se aporta el parámetro `p.value`, yo esperaba aquí encontrar que se filtrara la tabla en función de la columna correspondiente al p-valor ajustado, pero parece que ordena por logFC sin hacer ningún filtrado antes.
:::
```r=+
upgenes <- upgenes[!is.na(upgenes[, col, drop = F]), ]
upgenes <- upgenes[upgenes[, col] %in% BINGspace$Entrez.ID, ]
if (nrow(upgenes) > 150){
upgenes <- upgenes[order(-1*upgenes$logFC), col][1:150]
}else if(nrow(upgenes) < 10){
cat(paste("Less than 10 significantly differentially expressed upgenes found in BINGspace for ", name, ". Analysis stopped", sep = ""), sep = "\n"
, file = log, append = T)
return(NULL)
}else{
upgenes <- upgenes[order(-1*upgenes$logFC), col]
}
```
3. Hace lo mismo con los genes down. Sin embargo, en este caso, si hay menos de 10 genes la query se hace igualmente, dando un aviso al usuario de que no se van a enviar los genes downregulados:
```r=+
}else if(nrow(dngenes) < 10){
cat(paste("Less than 10 significantly differentially expressed dngenes found in BINGspace for ", name,
". Skipping dngenes list", sep = ""), sep = "\n", file = log, append = T)
writeDngenes <- FALSE
}
```
5. Escribe los ficheros gmt:
```r=+
upgenes <- c(paste(name, "_up", sep = ""), "TAG", upgenes)
dngenes <- c(paste(name, "_dn", sep = ""), "TAG", dngenes)
fwrite(as.list(upgenes), file = paste(direc,"/gmtFiles/",name, "UP.gmt", sep = ""), sep = "\t")
if(writeDngenes){
fwrite(as.list(dngenes), file = paste(direc,"/gmtFiles/",name, "DN.gmt", sep = ""), sep = "\t")
cat(paste("_up and _dn GMT files for ",direc,"/gmtFiles/", name, " created.", sep = ""), sep = "\n"
, file = log, append = T)
}
else{
cat(paste("_up GMT files for ", name, " created.", sep = ""), sep = "\n"
, file = log, append = T)
}
}
```
## Distribución de genes up y down
Guardo la tabla con el siguiente comando:
```sql
select node_id, geneSet, count(entrez_id)
from topTable
group by node_id, geneSet
order by node_id
into outfile '/tmp/topTable_geneSet_numbers.txt';
```
Y ahora vamos a mirar un poco. Vemos que para casi todos los nodos tenemos que el geneSet 'NULL' tiene muchísimos genes:
Mientras tanto, los geneSets UP y DN tienen menos:
> La línea horizontal marca dónde está el 150.
Podemos hacer un histograma para ver cuántos nodos tienen < 150 genes en alguno de estos geneSets. Este lo he hecho restringiendo para ver únicamente los geneSets con un tamaño < 200, porque si no no se ve nada:
> Las líneas verticales marcan el 10 (mínimo para enviar una query) y el 150 (máximo de genes por query).
En la gráfica anterior podemos comprobar que __no hay geneSets con un tamaño < 10__. De hecho, si hacemos un `summary`:
```
Min. 1st Qu. Median Mean 3rd Qu. Max.
36.0 527.5 665.0 10497.1 7276.0 239830.0
```
Por otro lado, tenemos 46 geneSets UP y 47 DN con < 150 genes (44 solapan entre ambos).
:::warning
No termino de entender ahora por qué esto sale así, porque en la tabla `cmap_foodrugs` solo tenemos 398 node_id y aquí hay 410. Si todos estos gene sets tienen más de 10 genes, no entiendo por qué no tienen conexión con cmap...
:::
En la intersección de los `node_id` de ambas tablas tenemos 398 nodos, que son los que aparecen en `cmap_foodrugs`. Los que están solo en `topTable` son:
```r
unique(genesets$node_id)[!unique(genesets$node_id) %in% unique(cmap_foodrugs$node_id)]
[1] 111 203 209 283 318 326 327 328 329 357 395 456
```
Los nodos de `topTable` con menos de 3 geneSets son:
```r
names(which(table(genesets$node_id) < 3))
[1] "111" "203" "209" "283" "357" "359" "395"
```
> Todos estos tienen en común que no tienen DN. Algunos tienen solo NULL (203, 357, 395) y otros tienen NULL y UP (111, 209, 283, 359).
* Revisando estos nodos sí es verdad que la mitad no tienen todos los geneSets, bien porque solo tengan NULL o NULL y UP.
* Ninguno de ellos tiene geneSet DN.
* Sin embargo, hay nodos a los que les falta algún geneSet que no están en esta lista (nodo 359), y hay nodos de la lista con 3 geneSets (318, 326, 327, 328, 329, 456)
* Por tanto, el número de geneSets no es el motivo por el cual faltan `node_id` en `cmap_foodrugs`.
Lo que se me ocurre ahora es mirar los p-values... No vaya a ser que sí haya un filtro que no hayamos visto, y en realidad sean los que no tenan como mínimo x valores con p < 0.05, o lo que sea.
:::danger
Pausa aquí porque no funcionaba la VPN :) voy a intentar configurar la base de datos en la partición de Windows del portátil
:::
Mientras tanto otra opción es buscar estos nodos en la página web y ver qué onda con ellos. Para saber las condiciones de cada nodo he ido a una tabla que tenía guardada en el directorio '/fnsCloud-FooDrugs/misc-stuff' llamada `nodes_sample_NJ`:
> **Están en `topTable` pero no en `cmap_foodrugs`:** 111 203 209 283 318 326 327 328 329 357 395 456
> **Tienen menos de 3 geneSets**: "111" "203" "209" "283" "357" "359" "395"
| node_id | treatment | concentration | time_point |
| ------- | ---------------------------------- | ------------- | ---------- |
| 111 | dasatinib 12 weeks | - | 13 weeks |
| 203 | grape extract | 8 mg | 0 |
| 209 | heavy drinking | - | - |
| 283 | moderate drinking | - | - |
| 318 | perifosine | 10 uM | 24 hours |
| 326 | prochlorperazine | - | 24 hours |
| 327 | prodrug g668684 | 7 uM | 72 hours |
| 328 | prodrug g669043 | 1500 nm | 72 hours |
| 329 | prodrug g669043 | 150 nm | 72 hours |
| 357 | resveratrol-enriched grape extract | 8 mg | 0 |
| 359 | resveratrol-enriched grape extract | 16 mg | 12 months |
| 395 | sorafenib | 4 uM | 72 hours |
| 456 | tvb3567 | | 48 hours |
| 456 | tvb3567 | | 48 hours |
| 456 | tvb3567 | | 48 hours |
> en el nodo 456 tenemos tres filas "iguales" porque vienen de líneas celulares diferentes :)
# Red
## Intentando entender tabla `cmap_foodrugs`
No entiendo por qué hay filas duplicadas... Vamos a ver cuántas hay y a ver si viéndolas puedo sacar algo en claro :male-detective:
```bash
## directorio: ~/fnscloud-FooDrugs/misc-stuff
# cmap_foodrugs.txt es simplemente la tabla cmap_foodrugs de mySQL exportada
cat cmap_foodrugs.txt | sort | uniq | wc -l
# 28,216,210
wc -l cmap_foodrugs.txt
# 28,683,990 cmap_foodrugs.txt
cat cmap_foodrugs.txt | sort | uniq -d | wc -l
# 396,885
cat cmap_foodrugs.txt | sort | uniq -d | head
# node_id cmap_node_id tau
# 110 0 0
# 110 10 0
# 110 100 0
# 110 1000 0
# 110 10000 0
# 110 10001 -71.5964
# 110 10002 -16.3902
# 110 10003 88.0374
# 110 10004 0
# 110 10005 0
# separamos lineas duplicadas a un archivo aparte:
cat cmap_foodrugs.txt | sort | uniq -d > cmap_foodrugs_duplines.txt
```
Si miramos el archivo que acabamos de generar para almacenar todas las líneas duplicadas (`cmap_foodrugs_duplines.txt`), veremos que hay **6 nodos de GEO**, con valores de tau para **todos los nodos de CMAP**:
```bash
# vamos a ver que node_id unicos salen:
cut -f1 cmap_foodrugs_duplines.txt | uniq
# 110
# 112
# 113
# 114
# 115
# 118
# y que cmap_node_id unicos tenemos...
# en este caso si hace falta sort al ser la 2a columna
cut -f2 cmap_foodrugs_duplines.txt | sort | uniq | wc -l
# 70,895 --> TODOS los nodos de CMAP
# las taus:
cut -f3 cmap_foodrugs_duplines.txt | sort | uniq | wc -l
# 104,505
```
De aquí puedo entender lo de CMAP, porque al final tenemos que probar todas las posibles interacciones. Lo que no entiendo es que solo haya 6 nodos de los obtenidos por nosotros para los cuales tengamos filas repetidas...
> Si miramos nodo a nodo, todos tienen 70895 repeticiones menos el 118, que se repite solo 42410 veces. Así que deberíamos tener $70895 · 5 + 42410 = 354475 + 42410 = 396885$ líneas repetidas, que es exactamente lo que encontramos! :confetti_ball:
Vamos a mirar, para el `node_id = 110`, cuántas repeticiones tenemos:
```sql
select node_id, cmap_node_id, tau
from cmap_foodrugs
where node_id = 110
order by cmap_node_id;
-- 141,790 filas
select distinct node_id, cmap_node_id, tau
from cmap_foodrugs
where node_id = 110
order by cmap_node_id;
-- 70,895 filas
```
Es decir... Este nodo tiene 70,895 conexiones posibles con nodos de CMAP (tantos como nodos de CMAP hay). Cada una de estas tiene un tau asociado. Encontramos que cada una de estas conexiones, definida por una combinación `node_id`-`cmap_node_id`-`tau`, se repite dos veces en la tabla `cmap_foodrugs`.
Ahora podemos repetir el proceso para todos los nodos con repeticiones modificando un poco el comando anterior, y guardando un nuevo fichero de texto de cuatro columnas: número de repeticiones de la fila, `node_id`, `cmap_node_id` y `tau`. Vamos a ver:
```bash
cat cmap_foodrugs.txt | sort | uniq -dc > cmap_foodrugs_duplines_count.txt
cut -f1 cmap_foodrugs_duplines_count.txt | sort | uniq
# 2 110
# 2 112
# 2 114
# 2 115
# 2 118
# 3 113
```
Sigo sin tener ni idea de qué pasa :male-detective:
Simplemente por tener la información de los estudios, concentraciones, etc. correspondientes a cada uno de los nodos de la discordia:
```sql
select distinct node_id, study_id, treatment, origin_name, time_point, concentration
from nodes natural join sample
where node_id = 110 or node_id = 112 or node_id = 113 or node_id = 114 or node_id = 115 or node_id = 118;
```
```
+---------+----------+--------------------+-------------+------------+---------------+
| node_id | study_id | treatment | origin_name | time_point | concentration |
+---------+----------+--------------------+-------------+------------+---------------+
| 110 | 140 | dasatinib 12 weeks | Kasumi-1 | 12 weeks | - |
| 112 | 133 | decitabine | HCT-116 | 72 hours | 5 µM |
| 113 | 133 | decitabine | HCT-116 | 72 hours | 120 nm |
| 114 | 89 | Delta tocotrienol | HeLa | 24 hours | 2 µg/ml |
| 115 | 60 | deoxycholic acid | MCF7 | 12 hours | 10 µM |
| 118 | 157 | diclofenac | VMM39 | 8 hours | - |
| 118 | 157 | diclofenac | SLM2 | 8 hours | - |
| 118 | 157 | diclofenac | DM331 | 8 hours | - |
+---------+----------+--------------------+-------------+------------+---------------+
8 rows in set (0,00 sec)
```
# Interacciones en ambas estrategias
## Resveratrol
## Todas las interacciones posibles
**Todas las interacciones obtenidas por GEO + CMAP.** Estos comandos tardan muchísimo, recordemos que hay ~ 1 millón de edges con tau > 90. Pruebo también con > 95 para restringir un poco más:
```sql
select distinct compound,treatment
from cmap natural join cmap_foodrugs natural join nodes natural join sample
where abs(tau) > 90
into outfile '/tmp/all_CMAP_interactions_90.txt';
-- Query OK, 524081 rows affected (10 min 0,79 sec)
select distinct compound,treatment
from cmap natural join cmap_foodrugs natural join nodes natural join sample
where abs(tau) > 95
into outfile '/tmp/all_CMAP_interactions_95.txt';
-- Query OK, 299177 rows affected (9 min 52,66 sec)
```
**Todas las interacciones obtenidas por text mining.** Esto es mucho más sencillo:
```sql
select distinct food, drug
from TM_interactions
into outfile '/tmp/all_TM_interactions.txt';
-- Query OK, 5487 rows affected (0,01 sec)
```
:::spoiler
```bash
mkdir fnscloud-FooDrugs/misc-stuff/comparing-interactions
sudo mv '/tmp/all_CMAP_interactions_90.txt' fnscloud-FooDrugs/misc-stuff/comparing-interactions
sudo mv '/tmp/all_CMAP_interactions_95.txt' fnscloud-FooDrugs/misc-stuff/comparing-interactions
sudo chmod ugo+rw fnscloud-FooDrugs/misc-stuff/comparing-interactions/all_CMAP_interactions_9*
sudo mv '/tmp/all_TM_interactions.txt' fnscloud-FooDrugs/misc-stuff/comparing-interactions/
sudo chmod ugo+rw fnscloud-FooDrugs/misc-stuff/comparing-interactions/all_TM_interactions.txt
```
:::
Primero vamos a ordenar todos los archivos:
```bash
sort all_CMAP_interactions_95.txt > sorted_95.txt
sort all_TM_interactions.txt > sorted_TM.txt
sort all_CMAP_interactions_90.txt > sorted_90.txt
```
Lo más sencillo ahora es buscar las líneas en común entre cualquiera de los ficheros de la red y el de TM:
```bash
comm -12 sorted_TM.txt sorted_90.txt
```
```
coumestrol genistein
daidzein genistein
genistein estradiol
myricetin quercetin
naringenin doxorubicin
nobiletin atorvastatin
phloretin quercetin
pinocembrin quercetin
piperine resveratrol
quercetin atorvastatin
```
```bash
comm -12 sorted_TM.txt sorted_95.txt
```
```
coumestrol genistein
myricetin quercetin
nobiletin atorvastatin
phloretin quercetin
piperine resveratrol
quercetin atorvastatin
```
```bash
comm -12 sorted_TM.txt sorted_90.txt > comm_TM_90.txt
```
Pero... También puede ser que haya compuestos que CMap haya detectado como 'compound' ( = fármaco) y el TM como 'food', y eso no lo veremos aquí. Con permutar las columnas de uno de los ficheros nos vale. Permutaremos el de text mining, para no tener que permutar ambos archivos de CMAP:
```bash
paste <(cut -f2 sorted_TM.txt) <(cut -f1 sorted_TM.txt) > permuted_TM.txt
sort permuted_TM.txt > sorted_permuted_TM.txt
rm permuted_TM.txt
```
Y volvemos a comparar:
```bash
comm -12 sorted_permuted_TM.txt sorted_90.txt
```
```
acitretin vitamin A
atorvastatin quercetin
dasatinib quercetin
erlotinib genistein
estradiol genistein
midazolam quercetin
mitoxantrone quercetin
nocodazole quercetin
pantoprazole quercetin
prazosin quercetin
pterostilbene quercetin
retinol vitamin D
rosiglitazone quercetin
simvastatin quercetin
sulfasalazine quercetin
tamoxifen estradiol
tamoxifen estrogen
tamoxifen genistein
valsartan resveratrol
```
Lo guardamos en otro fichero:
```bash
comm -12 sorted_permuted_TM.txt sorted_90.txt > comm_pTM_90.txt
```
Y los concatenamos:
```bash
cat comm_TM_90.txt comm_pTM_90.txt | sort > final_TMvs90.txt
```
:::info
El fichero final tiene **29 líneas**...
| Compound (drug) | Treatment (food) |
| ------------- | ------------- |
acitretin |vitamin A
atorvastatin| quercetin
coumestrol |genistein
daidzein |genistein
dasatinib |quercetin
erlotinib |genistein
estradiol |genistein
genistein |estradiol
midazolam |quercetin
mitoxantrone| quercetin
myricetin | quercetin
naringenin |doxorubicin
nobiletin |atorvastatin
nocodazole |quercetin
pantoprazole| quercetin
phloretin |quercetin
pinocembrin |quercetin
piperine |resveratrol
prazosin |quercetin
pterostilbene| quercetin
quercetin |atorvastatin
retinol |vitamin D
rosiglitazone| quercetin
simvastatin |quercetin
sulfasalazine| quercetin
tamoxifen |estradiol
tamoxifen |estrogen
tamoxifen |genistein
valsartan |resveratrol
:::
# Tablas para el paper
## GEO
A partir de lo que había hecho Óscar en su HackMD he sacado una tabla más o menos parecida, pero sin lo de 'other':
```sql
SELECT treatment, COUNT(DISTINCT concat(lower(treatment), concentration,time_point)) AS n_treatments,
COUNT(DISTINCT study_id) AS n_series
FROM sample
GROUP BY treatment
ORDER BY n_treatments DESC
INTO OUTFILE '/tmp/most_common_treatments.txt';
```
```bash
sudo mv /tmp/most_common_treatments.txt fnscloud-FooDrugs/misc-stuff/
sudo chmod ugo+rw fnscloud-FooDrugs/misc-stuff/most_common_treatments.txt
```
Las 10 primeras filas:
| treatment | experiments | series |
|------------------ | ----------- | ---------- |
| control |58 |153 |
| vitamin D |17 |17 |
| genistein |13 |2 |
| quercetin |11 |4 |
| frankincense |9 |1 |
| etoposide |7 |1 |
| methanesulfonate |7 |1 |
| resveratrol |5 |5 |
| vorinostat |5 |3 |
| doxorubicin |4 |2 |
Top 10 de tratamientos más frecuentes en la tabla “sample”, ordenados por el número de estudios en el que aparecen. Entre ellos aparecen control y ref para referirse a las muestras sin tratar. En total tenemos 370 tratamientos diferentes, de los cuales 346 aparecen en un solo estudio.
> Ojo: esto es un poco tramposo, porque aquí aparecen como tratamientos distintos gallic_acid y Gallic Acid, por ejemplo. Óscar lo ha contemplado en su HackMD: https://hackmd.io/@OscarPiette/ryTvmiCTc#Resultado3
## Text mining
Para las tablas de esta parte veo que Óscar se ha peleado bastante y ha tenido que hacer muchs cambios manuales. Me fío de él para esto: https://hackmd.io/@OscarPiette/ryTvmiCTc#Resultado1
> Aquí ojo porque en `n_texts` no sé si tenemos referencias a un mismo texto como si fueran de textos distintos... Y lo de poner "dopamine" como food no sé yo si es muy adecuado.
## Mirando categorías de todos
:::info
Creo que aquí nos puede ser útil hacer algo como la **tabla 2** del TFM de Marco:
- Entiendo que para los compuestos de **GEO / CMap** esto ya está hecho? Al menos se puede ver en la página web... La cosa es que no sé de primeras dónde podemos encontrarlo.
- me extraña que no haya una categoría para "polyphenols" (sé que los flavonoids lo son, pero ¿?)
- ¿Sería más o menos sencillo hacer esto también para los términos de la parte de text mining?
:::
Vamos a empezar por obtener los foods que tenemos en CMAP y en Text mining por separado, y a ver cuál es su solapamiento. Como tengo las dos tablas (`sample` y `TM_interactions`) en `~/fnscloud-FooDrugs/misc-stuff`, las importo a R para que me sea más fácil y sigo por ahí. Scrit: `overlap_food_drugs.R`
:::danger
Aquí se me complica un poco la cosa, porque no siempre coincide la consideracón de algo como **food** o como **drug**. Pero claro, si consideramos todos los tratamiento de GEO y compuestos de CMap por un lado y todos los foo y drug de text mining por otro, acabamos con muchísimos compuestos, y muchos de los cuales no nos interesan.
- Posible solución: usar solo los `pert_type == trt_cp`. El resto no nos interesan y así trabajaríamos con 2545 compuestos (de CMap) en lugar de 6395.
- **¿será suficiente con hacer esto?**
> Por ahora voy a seguir **solo con GEO** y con todos los de text mining (tanto food como drugs), porque realmente lo que puede "fallar" es eso (en CMap consideramo que todo son drugs)
:::
Así que lo que voy a mandar a KEGG phytochemicals es:
- La lista de foods y drugs de text mining (`text_mining_both.txt`);
- también solo la lista de foods (`text_mining_foods.txt`);
- la lista de treatments de GEO (`geo_treatments.txt`);
- la unión entre las dos listas anteriores (_teniendo en cuenta foods **y** drugs_) (`text_mining_geo_union.txt`);
- la intersección de las dos listas (`text_mining_geo_intersection.txt`);
- **vamos a incorporar también los tratamientos tipo _compound_ del CMap (`all_foods_and_drugs.txt`).**
:::danger
El razonamiento detrás de esto es que si incorporo los drugs del text mining, tendré que incorporar también los de CMap...
:::
Con esto podemos ver:
- si el text mining está muy sesgado con respecto a lo que hay en GEO (comparamos ambas listas de TM con GEO, si salen cosas parecidas guay y si no también será interesante porque nos dirá si la detección de foods o de drugs funciona mejor que la otra);
- si la lista de interacciones que coinciden entre ambos recursos está sesgada con respecto a lo que tenemos en la base de datos (ej. si en la base de datos hubiera un 20% de flavonoides y en la lista de interacciones en común fueran el 100%).
:::warning
Pasos para hacer esto:
1. Mandar todo a KEGG --> podemos hacer por separado GEO, TM y CMap y ya luego unir los text files resultantes
2. A ver cómo se hacía la estadística del enriquecimiento este que pretendemos hacer
3. Hacer la estadítica :-)
:::
Sign in with Wallet
Connect another wallet