# TFM Blanca :smile:
###### tags: `Estudiantes` `WGS`
[ToC]
## Análisis metagenomas WGS.
Yo he estado trabajando con [MetaWRAP](https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md). Mi idea era seguir este pipeline completo, pero no tenemos capacidad para hacer la asignación de MAGs :sweat:. Así que decidí hacer asignación taxonómica con metaphlan... y ya puedo luego comparar mis datos con lo que esté en *curatedMetagenomicData* :smirk:
El procesamiento para *metaphlan* :custard: que he usado en mi proyecto de celiacos es el siguiente:
1. Ver QC y quitar adaptadores.... (MetaWRAP usa Read_QC y Trimgallore)
2. Filtrar seq. humanas. **Con BMTagger** Descargar: hg38 index.
3. Asignación taxonómica - [Metaphlan3](https://github.com/biobakery/MetaPhlAn/wiki/MetaPhlAn-3.0) :custard:
4. Análisis funcional - Humann3 :woman_in_lotus_position:
> **Nota sobre paired-end**
MetaPhlAn can also natively handle paired-end metagenomes (but does not use the paired-end information), and, more generally, metagenomes stored in multiple files (but you need to specify the --bowtie2out parameter):
---
<font color = 'steelblue'>
<p style="text-align:center;"><b><i>Comentarios de Blanca</i></b></p>
</font>
#### Comprobación de los ficheros
Para asegurarnos de que la descarga ha ido bien, comparamos los **códigos md5** de los ficheros descargados con los que venían en el ENA. Para ello:
- Sacamos los archivos md5 de cada fastq con este comando (por ahora **sigo con los `.gz`**):
```bash
for file in $(find -mindepth 2 -type f); do echo $file; md5sum $file > $file.md5; done # pongo el echo para poder seguir el proceso
```
- Lo concatenamos todo en un solo fichero:
```bash
for file in $(find -mindepth 2 -type f | grep md5); do cat $file >> md5_after_unzipping.txt ; done
```
- Ahora tenemos que generar el fichero `md5_paths_from_ENA.txt`, que tiene esta pinta:
```
9676cfdd5e318bcbe70219d314089925 ./ERR414230/ERR414230_1.fastq.gz
3ee413f2e6cb17dbf5deac241e19f512 ./ERR414231/ERR414231_1.fastq.gz
be2a52a759d8cd172714494d1eeb8f69 ./ERR414232/ERR414232_1.fastq.gz
df45c11ad71f7b5c1ec5d7f6eac08ff2 ./ERR414233/ERR414233_1.fastq.gz
c9d2db34f482833801f9df5535b3be86 ./ERR414234/ERR414234_1.fastq.gz
b27940088b9b8df21396f39d986f3637 ./ERR414235/ERR414235_1.fastq.gz
```
> - En este caso, yo generé el fichero así:
> ```bash
>cat subsetENA_ftp_md5.tsv | awk '{print $4 " " $5}' > tmp45.txt #este archivo tiene solamente los archivos que queremos
>wc -l tmp45.txt
>tail -n 176 tmp45.txt > tmp_nohead.txt # queremos evitar encabezado
>rm tmp45.txt
>sed 's/ftp.sra.ebi.ac.uk\/vol1\/fastq\/ERR..././g' tmp_nohead.txt > tmp_paths.txt # elimina parte del path que no queremos
>rm tmp_nohead.txt
>grep -v ';' tmp_paths > part1.txt
>grep ';' tmp_paths | tr ' ' ';' > semicolon.txt # columnas 1,3 son read1 y 2,4 son read2
>cut -d ';' -f1,3 semicolon.txt | tr ';' ' ' > read1.txt
>cut -d ';' -f2,4 semicolon.txt | tr ';' ' ' > read2.txt
>cat read1.txt read2.txt part.txt > md5_paths_from_ENA.txt
>cp md5_paths_from_ENA.txt samples
>cd samples
>```
> Para las muestras del PRJEB2054:
> ```bash
>cat filereport_read_run_PRJEB2054_tsv.txt | grep 'bgi-MH' | cut -f29 > fastq_md5.txt #obtiene columna fastq_md5
>cat filereport_read_run_PRJEB2054_tsv.txt | grep 'bgi-MH' | cut -f30 > fastq_ftp.txt #columna fastq_ftp
>paste fastq_md5.txt fastq_ftp.txt > both.txt #aqui podemos eliminar los dos .txt anteriores
>sed -i 's/ftp.sra.ebi.ac.uk\/vol1\/fastq\/ERR..././g' both.txt
>cut -d ';' -f1,4 both.txt | tr ';' ' ' > read.txt
>cut -d ';' -f2,5 both.txt | tr ';' ' ' > read1.txt
>cut -d ';' -f3,6 both.txt | tr ';' ' ' > read2.txt
>cat read.txt read1.txt read2.txt > md5_paths_from_ENA.txt
>mv md5_paths_from_ENA.txt RAW_READS
>cd RAW_READS
>```
:::success
Finalmente, el comando que corremos es:
```bash
md5sum --check md5_paths_from_ENA.txt > md5sum_check_log.txt 2>&1
```
En el fichero `md5sum_check_log.txt` podremos ver si los ficheros se descargaron bien:
```
./ERR414230/ERR414230_1.fastq.gz: OK
```
O si no fue el caso:
```
./ERR414241/ERR414241_1.fastq.gz: FAILED
```
:::
> Para las muestras del PRJEB5024 me salen varios `FAILED open or read` porque me olvidé de **excluir las muestras con `lane`** de los MH0006 y MH0012. Con una comprobación rápida podemos ver que los FAILED coinciden con estos ficheros y que no hay de qué preocuparse:
> ```bash
> cat filereport_read_run_PRJEB2054_tsv.txt | grep 'lane' | cut -f6 #vemos los ERR de esos ficheros
> cat md5sum_check_log.txt | grep 'FAILED open' | cut -d '/' -f2 | sort | uniq # OK!
> ```
> Ha habido tres que sí que han fallado:
> ```bash
> cat md5sum_check_log.txt | grep 'FAILED' | grep -v open
> ./ERR011189/ERR011189_1.fastq.gz: FAILED
> ./ERR011088/ERR011088_2.fastq.gz: FAILED
> ./ERR011229/ERR011229_2.fastq.gz: FAILED
> ```
> Saco las ftp para intentar descargarlos con `wget`:
> - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011189/ERR011189_1.fastq.gz
> - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011088/ERR011088_2.fastq.gz
> - ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011229/ERR011229_2.fastq.gz
>
> Tampoco ha funcionado, revisar qué pasa aquí. De momento las guardo en un directorio aparte y ya pensaré en ellas.
>
>Por otro lado, en estas muestras se da el problema de que tenemos un tercer fichero con el que no sé qué pasa.
> - Mirando con el run ERR011087, los \_1 y \_2 tienen 46560224 líneas, mientras que el tercero en discordia tiene 3680 (12652 veces menos líneas). Si hago grep con el id del run, son 11640056 lecturas frente a 920 (lo mismo entre 4).
> - Mirando los ID, haciendo head, tail, etc., no veo que los reads del tercer fichero aparezcan en los paired. Creo que son lecturas que no parearon con los otros dos ficheros, ahora queda **decidir qué hacemos con ellas.**
> - De momento voy a correr metaWRAP con las paired y ya veremos qué pasa con las otras.
#### Configuración del entorno e instalaciones
Para la instalación seguí las [instrucciones del github](https://github.com/bxlab/metaWRAP#installation)
Una cosa que tuve que hacer durante la instalación fue **descargar e indexar el hg38**. Seguí las [instrucciones de metaWRAP](https://github.com/bxlab/metaWRAP/blob/master/installation/database_installation.md#making-host-genome-index-for-bmtagger). El path que puse a la DB en el fichero config-metawrap es el siguiente:
```
BMTAGGER_DB=/home/lmarcos/BMTAGGER_INDEX
```
Los pasos que seguí para crear el entorno fueron:
```bash=v
conda create --name metawrap-2 python=2.7
conda activate metawrap-2
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
conda config --add channels ursky
conda install --only-deps -c ursky metawrap-mg # tarda un poco pero lo saca
```
<font color = 'steelblue'>
<p style="text-align:center;"><b><i>Fin comentarios de Blanca</i></b></p>
</font>
-------
### Mis scripts y cosas
#### MetaWRAP
Sacar data de c/carpeta y juntar en una. Luego descomprimir
Cambiar nombre/ quitar toda la :shit: del nombre
> De D1807_ESFP200009043-1a_HCYJCDSXY_L2_1.fq
a
D1807_1.fq
```bash=v
# movemos todos los fastq a carpeta reads:
fastqfiles=$(find -mindepth 2 -type f | grep -v md5)
for file in $fastqfiles; do mv $file reads; done
ls reads | wc -l # todo ok, hay 250 (264 - 14 malos)
# descomprimimos:
cd reads
gunzip *.gz
```
> - Para PRJEB5024:
```bash!
# movemos todos los fastq a carpeta RAW_READS:
fastqfiles=$(find -mindepth 2 -type f | grep -v md5 | grep '_')
for file in $fastqfiles; do mv $file ../RAW_READS; done
ls ../RAW_READS | wc -l # todo ok, hay 356: 178 runs descargados x 2 (forward + reverse) (178 porque partiamos de 186 bgi-MH - 8 con 'lane' de 06 y 12)
# descomprimimos:
cd ../RAW_READS
gunzip *.gz
```
- Limpiar y filtrar reads del host... (Paciencia :hourglass_flowing_sand:)
```bash
for F in RAW_READS/*_1.fq; do
R=${F%_*}_2.fq
BASE=${F##*/}
SAMPLE=${BASE%_*}
metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE &
done
```
:::success
```bash!
export PATH=$PATH:/home/lmarcos/metaWRAP/bin
for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_03-11-2021-10-40.txt 2>&1
for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_19-11-2021-18-23.txt 2>&1
```
:::
-------
<font color = 'steelblue'>
<p style="text-align:center;"><b><i>Comentarios de Blanca</i></b></p>
</font>
#### Ejecución del comando
Antes de correr todo esto hay que __añadir metawrap al PATH__. No me atreví a cambiar el path de forma definitiva por si acaso la liaba con algo, que alguna vez con mi ordenador me ha pasado. El comando que he usado es:
:::success
```export PATH=$PATH:/home/lmarcos/metaWRAP/bin```
:::
> ***Nota:*** al correr `which config-metawrap` en este nuevo entorno y revisar el fichero, parece que sigue apuntando correctamente al hg18.
Modifiqué un poco el comando porque los ficheros se han descomprimido como `.fastq`, no como `.fq`. Además, **antes de correrlo creé el directorio `READ_QC` porque me daba problemas si no lo hacía.** El comando que puse era:
```bash!
for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE & done > meta20211102.log
```
> <font color = 'gray'>_3 nov 10:20_</font> Esto no ha funcionado. Pongo una prueba pequeña otra vez, no vaya a ser que haya cambiado algo del environment sin darme cuenta:
> ```bash!
> for F in RAW_READS_prueba/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC_prueba/$SAMPLE & done > read_qcERROR.txt 2>&1
> ```
>Esta prueba sí funciona. Mirando en el github de metaWRAP hubo una persona a la que le pasó algo parecido y le dijeron que [a veces hay problemas con el "&" para paralelizar](https://github.com/bxlab/metaWRAP/issues/344).
:::success
Probamos sin `&` (y usamos **`2>&1`** para **redirigir el standard error además del standard output**):
```bash!
for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_03-11-2021-10-40.txt 2>&1
for F in RAW_READS/*_1.fastq; do R=${F%_*}_2.fastq; BASE=${F##*/}; SAMPLE=${BASE%_*}; metawrap read_qc -1 $F -2 $R -t 6 -o READ_QC/$SAMPLE; done > read_qcERROR_19-11-2021-18-23.txt 2>&1
```
>Este comando corrió sin problemas. Estuvo hasta las 21:40 h del 06/11 --> unas 83 horas, entre 75 muestras = como **1 h 6 min por muestra**
>
> ___Healthy MetaHit:___ corre hasta el 22/11/2021 a las 7:20 h, salen 175 carpetas. No aparecen mensajes de error. Habíamos descartado los runs ERR011189_1, ERR011088_2 y ERR011229_2, así que nos sale para las dos reverse que no existen (y por tanto no se hizo su control de calidad).
:::
:::warning
>Algunas cosas que detecté en el run del fin de semana del 30/oct-1/nov (__*nota:*__ las he anotado por si nos pudieran ayudar, aunque no sean de la ejecución buena):
- Al principio salen errores del estilo `RAW_READS/(...)_2.fastq file does not exist. Exiting...` con los siguientes archivos:
- ERR414342_2.fastq
- ERR414248_2.fastq
- ERR414346_2.fastq
- ERR414240_2.fastq
- ERR414269_2.fastq
- Estos archivos no existen en el directorio RAW_READS porque son los que tienen un código md5 erróneo, si nos fijamos los 5 aparecen en el directorio `samples-bad-md5`, que es donde almacenamos estos archivos. Hay otros archivos en este directorio, pero imagino que no nos salió el mensaje de error porque son terminados en `_1`, y el comando busca solo la pareja entre los `_1` que encuentra.
- Esto significa que hay un total de 14 ficheros que no estamos mirando (5 en los que falta el _2 + 8 en los que falta el _1 + 1 single-end)
- No tengo muy claro qué pasa con los archivos de lecturas single-end, en principio con este comando no habríamos hecho nada con ellos...
- Una cosa que he comprobado es que en todos los single-end hay pocas lecturas, así que igual no perdemos mucha información si pasamos de ellos
- **Making Pre-QC report:** Salen 6 líneas con el mensaje `Something went wrong with making pre-QC fastqc report. Exiting.`
**Estos mensajes ya no salen en el run definitivo :confetti_ball: (el primero de los \_1 y \_2 sí)**
:::
### Selección de runs
#### MetaHit T2D
Hemos visto que había muestras a las que le correspondía más de un run. Lo que hemos decidido es:
- Las runs single-end se descartan, ya que solo contienen el 4.6% de las lecturas.
- De las paired-end, juntamos las lecturas que tengan el mismo sample accession y miramos cuántas lecturas se juntan en total:
- Salen un par de muestras con muchas menos lecturas de las demás (~5e+06 lecturas frente a 2.5e+07 lecturas como mínimo en todos los demás): SAMEA2338656 y SAMEA2338670
- Los runs correspondientes a estas muestras son ERR414270 y ERR414286. Los descartamos:
```bash
mv READ_QC/ERR414270 muestras-descartadas/low-readcount
mv READ_QC/ERR414286 muestras-descartadas/low-readcount
```
- Nos quedamos entonces con **73 runs,** con un total de **65 accessions únicos.**
- En cuanto a las muestras que tienen más de un run, lo que vamos a hacer es **juntar todas las reads en uno solo antes de pasárselo a metaphlan.**
- He preparado un excel con los 73 runs que queremos para poder ver claramente cuáles son de la misma muestra y juntarlos. Lo voy a hacer a mano porque para qué pasar 10 h automatizando algo que a mano haces en 10 minutos... :smile_cat::
- Lo que voy a hacer es crear directorios nuevos para estos runs. Por ejemplo, para los runs ERR414230 y ERR414231:
- Creo directorio `READ_QC/ERR414230_1`
- Concateno:
- `cat ERR414230/final_pure_reads_1.fastq ERR414231/final_pure_reads_1.fastq > ERR414230_1/final_pure_reads_1.fastq`
- Lo mismo para el `_2`
- Comprobamos suma líneas individual y concatenada con `wc -l` --> :-1:
- Cuando termine, muevo todos los directorios antiguos a `muestras-descartadas/a-concatenadar`
- Los run_accession de estas muestras están en un fichero de texto `Preprocessing/runs_a_juntar.txt`
| Sample_accesion| Run_accession |
| -------- | -------- |
| SAMEA2338623 | **ERR414230** |
| SAMEA2338623 | **ERR414231** |
| SAMEA2338629 | ERR414237 |
| SAMEA2338629 | ERR414238 |
| SAMEA2338643 | **ERR414253** |
| SAMEA2338643 | **ERR414254** |
| SAMEA2338646 | ERR414257 |
| SAMEA2338646 | ERR414258 |
| SAMEA2338647 | **ERR414259** |
| SAMEA2338647 | **ERR414260** |
| SAMEA2338657 | ERR414271 |
| SAMEA2338657 | ERR414272 |
| SAMEA2338683 | **ERR414299** |
| SAMEA2338683 | **ERR414300** |
| SAMEA2338730 | ERR414349 |
| SAMEA2338730 | ERR414350 |
He concatenado así (toma como 2 min por pareja de muestras):
```bash
mkdir ERR414230_1
cat ERR414230/final_pure_reads_1.fastq ERR414231/final_pure_reads_1.fastq > ERR414230_1/final_pure_reads_1.fastq
cat ERR414230/final_pure_reads_2.fastq ERR414231/final_pure_reads_2.fastq > ERR414230_1/final_pure_reads_2.fastq
mkdir ERR414237_8
cat ERR414237/final_pure_reads_1.fastq ERR414238/final_pure_reads_1.fastq > ERR414237_8/final_pure_reads_1.fastq
cat ERR414237/final_pure_reads_2.fastq ERR414238/final_pure_reads_2.fastq > ERR414237_8/final_pure_reads_2.fastq
mkdir ERR414253_4
cat ERR414253/final_pure_reads_1.fastq ERR414254/final_pure_reads_1.fastq > ERR414253_4/final_pure_reads_1.fastq
cat ERR414253/final_pure_reads_2.fastq ERR414254/final_pure_reads_2.fastq > ERR414253_4/final_pure_reads_2.fastq
mkdir ERR414257_8
cat ERR414257/final_pure_reads_1.fastq ERR414258/final_pure_reads_1.fastq > ERR414257_8/final_pure_reads_1.fastq
cat ERR414257/final_pure_reads_2.fastq ERR414258/final_pure_reads_2.fastq > ERR414257_8/final_pure_reads_2.fastq
mkdir ERR414259_0
cat ERR414259/final_pure_reads_1.fastq ERR414260/final_pure_reads_1.fastq > ERR414259_0/final_pure_reads_1.fastq
cat ERR414259/final_pure_reads_2.fastq ERR414260/final_pure_reads_2.fastq > ERR414259_0/final_pure_reads_2.fastq
mkdir ERR414271_2
cat ERR414271/final_pure_reads_1.fastq ERR414272/final_pure_reads_1.fastq > ERR414271_2/final_pure_reads_1.fastq
cat ERR414271/final_pure_reads_2.fastq ERR414272/final_pure_reads_2.fastq > ERR414271_2/final_pure_reads_2.fastq
mkdir ERR414299_0
cat ERR414299/final_pure_reads_1.fastq ERR414300/final_pure_reads_1.fastq > ERR414299_0/final_pure_reads_1.fastq
cat ERR414299/final_pure_reads_2.fastq ERR414300/final_pure_reads_2.fastq > ERR414299_0/final_pure_reads_2.fastq
mkdir ERR414349_0
cat ERR414349/final_pure_reads_1.fastq ERR414350/final_pure_reads_1.fastq > ERR414349_0/final_pure_reads_1.fastq
cat ERR414349/final_pure_reads_2.fastq ERR414350/final_pure_reads_2.fastq > ERR414349_0/final_pure_reads_2.fastq
```
Comprobamos las sumas:
- 230_1 :+1:
- 237_8 :+1:
- 253_4 :+1:
- 257_8 :+1:
- 259_0 :+1:
- 271_2 :+1:
- 299_0 :+1:
- 349_0 :+1:
COMPROBAR QUE EN LOS NORMALES EL wc -l SALE IGUAL PARA \_1 Y PARA \_2 :+1: Lo he mirado en el ERR414295, por ejemplo
Movemos a `muestras-descartadas/a-concatenar`:
- Corremos este comando en el directorio costea_dnk, al que habremos movido previamente el archivo runs_a_juntar.txt con los run_acc de los ficheros que hemos concatenado
- Hemos importado el runs_a_juntar.txt de windows, con lo que el programa falla debido a los retornos de carro, que son distintos en windows y linux. Para solucionarlo:
- `sed -i 's/\r//g' runs_a_juntar.txt`
```bash!
for file in $(cat runs_a_juntar.txt); do mv READ_QC/$file muestras-descartadas/a-concatenar/ ; done
```
:::warning
Una vez terminado el pipeline, los fastq de las muestras antes de concatenar ya no sirven de nada. Los eliminamos con:
```bash=v
#directorio:blanca/costea-dnk/muestras_descartadas/a-concatenar
find . -name 'final_pure_reads_?.fastq' | xargs rm
```
:::
#### MetaHit sanos - 22 noviembre
En este caso tenemos bastantes más ficheros que concatenar. Tengo que encontrar la forma de automatizarlo --> primero vamos a estudiar bien la tabla del ENA, a ver qué podemos hacer.
```bash=v
cut -f3,6 filereport_read_run_PRJEB2054_tsv.txt | awk 'NR>1' > sample-run-acc.txt # sampleacc + runacc
```
- Si hacemos `more_runs=$(cut -f1 sample-run-acc.txt | uniq -c | sort | tail -n 9 | awk '{print $2}')` obtendremos las muestras que tienen más de 2 runs, y sacamos una expresión regular con `runs=$(echo $more_runs | tr ' ' '|'`.
- Con `grep -vE $runs sample-run-acc.txt > two_runs.txt` nos quedamos con los ficheros con solo dos runs.
- Nos quedamos con el primer run de cada pareja con `awk 'NR % 2 == 1' two-runs.txt > first-of-each.txt`, y con `cut -f2 first-of-each.txt > runacc-first-of-each.txt` sacamos el run accession.
- Al final terminamos con un fichero `to-concatenate.txt` que contiene en cada línea los runs a concatenar (ERRetc1;ERRetc2) para aquellas muestras que cuentan con dos runs. Les pasamos el script `concatenator.sh` (sesión de screen `blanca-descargas`) con `./concatenator.sh to-concatenate.txt concatenated-runs/`:
```bash=v
#!/bin/bash
file=$1
conc_dir=$2
while IFS= read -r line; do
base_name=$(echo $line | cut -d ';' -f1)
for_1=$(echo $line | cut -d ';' -f1)/final_pure_reads_1.fastq
for_2=$(echo $line | cut -d ';' -f2)/final_pure_reads_1.fastq
rev_1=$(echo $line | cut -d ';' -f1)/final_pure_reads_2.fastq
rev_2=$(echo $line | cut -d ';' -f2)/final_pure_reads_2.fastq
cat $for_1 $for_2 > $conc_dir/$base_name.conc_1.fastq
cat $rev_1 $rev_2 > $conc_dir/$base_name.conc_2.fastq
rm $for_1 $for_2 $rev_1 $rev_2
done < $file
```
- El problema es que en el fichero de entrada que le hemos pasado tiene IDs de runs que no teníamos nosotros, por eso salen bastantes errores y obtenemos ficheros vacíos. Los eliminamos:
```bash
rm READ_QC/concatenated-runs/ERR0113*
rm READ_QC/concatenated-runs/ERR01128*
rm READ_QC/concatenated-runs/ERR01129*
# lo mismo para las 73, 75 y 79
```
Además, concatenamos estos dos (no aparecen porque es una de las muestras con muchos runs): :heavy_check_mark:
| Sample_accesion| Run_accession |
| -------- | -------- |
|SAMEA728680|ERR011099|
|SAMEA728680|ERR011100|
```bash=v
cat ERR011099/final_pure_reads_1.fastq ERR011100/final_pure_reads_1.fastq > concatenated-runs/ERR011099.conc_1.fastq
cat ERR011099/final_pure_reads_2.fastq ERR011100/final_pure_reads_2.fastq > concatenated-runs/ERR011099.conc_2.fastq
rm ERR011099/final_pure_reads_1.fastq ERR011100/final_pure_reads_1.fastq ERR011099/final_pure_reads_2.fastq ERR011100/final_pure_reads_2.fastq
```
y los ficheros con tres runs:
| Sample_accesion| Run_accession |
| -------- | -------- |
|**SAMEA728635**|ERR011089|
|**SAMEA728635**|ERR011090|
|**SAMEA728635**|ERR011091|
|SAMEA728854|ERR011092|
|SAMEA728854|ERR011093|
|SAMEA728854|ERR011094|
|**SAMEA728690**|ERR011109|
|**SAMEA728690**|ERR011110|
|**SAMEA728690**|ERR011111|
|SAMEA728667|ERR011114|
|SAMEA728667|ERR011115|
|SAMEA728667|ERR011116|
|**SAMEA728672**|ERR011117|
|**SAMEA728672**|ERR011118|
|**SAMEA728672**|ERR011119|
|SAMEA728658|ERR011126|
|SAMEA728658|ERR011127|
|SAMEA728658|ERR011128|
|**SAMEA728617**|ERR011131|
|**SAMEA728617**|ERR011132|
|**SAMEA728617**|ERR011133|
|SAMEA728749|ERR011236|
|SAMEA728749|ERR011237|
|SAMEA728749|ERR011238|
- Hacemos script que tenga los tres runs a concatenar en una misma línea, separados por ;
- Script `concatenator3.sh` (lo mismo pero añadiendo un tercer forward y reverse):
```bash=v
#!/bin/bash
file=$1
conc_dir=$2
while IFS= read -r line; do
base_name=$(echo $line | cut -d ';' -f1)
for_1=$(echo $line | cut -d ';' -f1)/final_pure_reads_1.fastq
for_2=$(echo $line | cut -d ';' -f2)/final_pure_reads_1.fastq
for_3=$(echo $line | cut -d ';' -f3)/final_pure_reads_1.fastq
rev_1=$(echo $line | cut -d ';' -f1)/final_pure_reads_2.fastq
rev_2=$(echo $line | cut -d ';' -f2)/final_pure_reads_2.fastq
rev_3=$(echo $line | cut -d ';' -f3)/final_pure_reads_2.fastq
cat $for_1 $for_2 $for_3 > $conc_dir/$base_name.conc_1.fastq
cat $rev_1 $rev_2 $rev_3 > $conc_dir/$base_name.conc_2.fastq
rm $for_1 $for_2 $for_3 $rev_1 $rev_2 $rev_3
done < $file
```
:::success
- En el archivo del ENA podemos comprobar que entre las muestras de daneses hay 85 sample accessions únicos. Por tanto, debería haber 170 archivos en `concatenated-runs` (85 muestras x 2 archivos de lecturas, forward y reverse) :+1:
- En cuanto a los ficheros relacionados con los runs ERR011189_1, ERR011088_2 y ERR011229_2, no tenemos sus lecturas limpias, por lo que en `concatenated-runs` solo queda el fichero con uno de los runs de la pareja (el otro, en este caso los runs ERR011188, ERR011087 y ERR011228).
- Parece que los números de líneas cuadran :+1:
:::
<font color = 'steelblue'>
<p style="text-align:center;"><b><i>Fin comentarios de Blanca</i></b></p>
</font>
---
- Mover las secuencias limpias y sin human-host
```bash
mkdir CLEAN_READS
for i in READ_QC/*; do
b=${i#*/}
mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastq
mv ${i}/final_pure_reads_2.fastq CLEAN_READS/${b}_2.fastq
done
```
```bash
for f in READ_QC/concatenated-runs/; do mv $f CLEAN_READS/; done
cd CLEAN_READS/concatenated-runs
mv * ..
cd ..
rm -r concatenated-runs
```
:::danger
**IMPORTANTE:** no nombrar los ficheros concatenados con una \_, porque el script de Perl que usamos para metaphlan se va a hacer un lío.
:::
:::warning
Quitamos la primera "\_" de los nombres:
```bash
conc_files=$(ls CLEAN_READS/ | grep -E '_._..fastq')
for file in $conc_files; do
new_name=$(sed 's/\_//' <(echo $file))
mv CLEAN_READS/$file CLEAN_READS/$new_name
done
```
Me he equivocado y he corrido un comando en el que la línea de `mv` era mv `CLEAN_READS/$file CLEAN_READS$new_name`. Movemos los ficheros que se han generado a CLEAN_READS y solucionamos los nombres:
```bash
clean_files=$(ls CLEAN_READS/ | grep CLEAN')
for file in $clean_files; do
new_name=$(sed 's/CLEAN_READS//' <(echo $file))
mv CLEAN_READS/$file CLEAN_READS/$new_name
done
```
:::
Ahora sí que sí. Tenemos 130 ficheros en CLEAN_READS, dos por cada uno de los 65 individuos. Los ficheros cuyo nombre tiene un carácter de más se corresponden con los que hemos concatenado.
### Metaphlan :custard:
#### Instalación:
Creamos entorno de conda:
```bash
conda create --name metaphlan
conda activate metaphlan
conda install -c bioconda metaphlan
```
Ahora tenemos que instalar también la base de datos:
> - "_If you have installed MetaPhlAn using Anaconda, it is advised to install the database in a folder outside the Conda environment. To do this, run:_
> `metaphlan --install --bowtie2db <database folder>`
>- _If you install the database in a different location, remember to run MetaPhlAn using `--bowtie2db <database folder>`!_"
>`metaphlan --install --bowtie2db ~/blanca/bowtiedb`
#### Comandos:
**Nota:** he modificado la línea:
```bash!
system("metaphlan $r1,$r2 --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt")
```
Para incluir lo de `--bowtie2db`:
```bash!
system("metaphlan $r1,$r2 --bowtie2db ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt")
```
De manera que al final el script `metaphlan-script.pl` tiene esta pinta:
:::success
```perl!
#!/usr/bin/perl
#dar de argumento solo los _1
foreach $ar (@ARGV)
{
# sacamos cada archivo:
@field = split (/_/, $ar);
$name = $ar; # name = nombre completo
$name =~ s/_1.fastq//; # name = nombre antes de _
# obtenemos el _2 correspondiente:
$r1 = $ar; # nombre completo de lectura con _1
$r2 = $ar; # lo mismo que r1 por ahora
$r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2
# pasamos ambos a metaphlan:
system("metaphlan $r1,$r2 --bowtie2db ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt");
}
```
:::
Me he atrevido a correrlo con todas las muestras a la vez (este script tiene que estar en la carpeta `CLEAN_READS`):
```bash!
./metaphlan-script.pl *_1.fastq > metaphlanOUT_08-11-2021-13-09.txt 2>&1
```
>**Nota 1:** antes de correrlo me he ido a la sesión de screen blanca (`screen -r blanca`), he desactivado el entorno metawrap-2, he activado el de metaphlan y ahí ya sí he ejecutado el comando
>**Nota 2:** nada más empezar a correr, en `CLEAN_READS` había 850 GB
> **Nota 3:** A las 15:25, lleva 7 archivos: va a unos 20 min por muestra --> 20 * 65 = 1300 min (21 h 40 min)
:::success
El script para las muestras de sanos cambia un poco, añadimos el argumento `--read_min_len 40` (ver diario 22/nov):
```perl!
#!/usr/bin/perl
#dar de argumento solo los _1
foreach $ar (@ARGV)
{
# sacamos cada archivo:
@field = split (/_/, $ar);
$name = $ar; # name = nombre completo
$name =~ s/_1.fastq//; # name = nombre antes de _
# obtenemos el _2 correspondiente:
$r1 = $ar; # nombre completo de lectura con _1
$r2 = $ar; # lo mismo que r1 por ahora
$r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2
# pasamos ambos a metaphlan:
system("metaphlan $r1,$r2 --bowtie2db ~/blanca/bowtiedb --bowtie2out $name.bowtie2.bz2 --nproc 6 --input_type fastq -o profiled_$name.txt --read_min_len 40");
}
```
```bash!
#./metaphlan-script.pl *_1.fastq > metaphlanOUT_22-11-2021-13-44.txt 2>&1
./metaphlan-script.pl *_1.fastq > metaphlanOUT_22-11-2021-16-23.txt 2>&1
```
:::
>- Sesión de screen `blanca-metaphlan`
>- Como daba fallos de 'ZeroDivisionError', he cambiado la longitud mínima para que entren todas las lecturas (hay ficheros donde la longitud de la lectura original es de 44 bp, y otros donde la longitud de original es 75 bp pero se ha reducido mucho al hacer el trimming).
:::warning
_Notas:_
- No hace falta hacer nada especial ni para el script ni para los argumentos (`$./script.pl *_1.fastq`)
- Sobre el comando de metaphlan:
- La entrada `metaphlan $r1,$r2` está en el formato necesario para manejar metagenomas paired-end
- El argumento `--bowtie2out` permite guardar la salida de bowtie2, para que si tenemos que volver a correr este proceso sea más rápido al ahorrar este paso.
- Para ello usaríamos `--input_type bowtie2out` (eso sí, tiene que estar generado con MetaPhlAn 3.0 en adelante, anteriores no son compatibles)
:::
- Unir resultados:
```bash!
python utils/merge_metaphlan_tables.py profiled_*.txt > output/merged_abundance_table.txt
```
> El comando que usé es ligeramente distinto:
> ```bash!
> merge_metaphlan_tables.py profiled_*.txt > ../OUTPUT/merged_abundance_table.txt
> ```
## Siguientes pasos
Entiendo que lo próximo es **juntar todas las OTUs generadas con las disponibles en curatedMetagenomicData** y ponernos a investigar.
### Juntar todos los datos
Aquí la primera pregunta que me surge es **cómo estratificar los datos.** No me queda claro si son cuatro grupos (MUO, MHO, MUNO, MHNO) o si primero miraremos MU vs. MH y ya veremos qué pasa con los obesos.
Otra cosa importante es mirar cómo importar en R la salida de metaphlan de manera que tenga el mismo formato que la información de curatedMetagenomicData.
- [Metaphlan to phyloseq](https://gist.github.com/lwaldron/512d1925a8102e921f05c5b25de7ec94), esto es del mismo que hizo curatedMetagenomicData
- Mirar también preguntas de biostars como [esta](https://www.biostars.org/p/449688/)
- En la página de curatedMetagenomicData hay ejemplos de cómo hacer análisis de [alfa diversidad](https://waldronlab.io/curatedMetagenomicData/articles/curatedMetagenomicData.html#alpha-diversity-1) y otros, utilizando paquetes como `mia` y `vegan`
- Según explican ellos mismos, no hay ningún método en `curatedMetagenomicData` que permita convertir datos a phyloseq, pero sí lo hay en `mia`: `makePhyloseqFromTreeSummarizedExperiment`
- Entiendo entonces que podemos ahorrarnos poner la salida de MetaPhlAn en el formato de curated y ponerlo directamente en formato de phyloseq, si lo de curated también lo tenemos que convertir...
Para los datos que vienen de MetaPhlAn :custard: he usado el script de Laura modificando alguna cosa (ver más abajo).
Para los datos de curatedMetagenomicData, según su [página](https://waldronlab.io/curatedMetagenomicData/articles/curatedMetagenomicData.html):
- Primero he instalado el paquete `mia` para utilizar la función que he mencionado más arriba, haciendo una instalación normal de bioconductor.
- Para sacar los datos de curated creo que hay que hacer algo así: `curatedMetagenomicData("AsnicarF_20.+.relative_abundance", dryrun = FALSE, counts = TRUE, rownames = "short")`
- `dryrun = FALSE` devuelve un `TreeSummarizedExperiment` o una lista de `SummarizedExperiment`
- `rownames = 'short'` indica que el nombre de las filas de `relative_abundance` es el nombre de la especie. Por defecto usa `'long'`, otra opción es `NCBI` para el NCBI Taxonomy ID.
- `counts = TRUE` multiplica las abundancias relativas por la profundidad de lectura y redondea al entero más cercano.
- Luego no me queda claro si hay que pipearlo a un `mergeData()`
- Otra opción es usar `returnSamples()`, que permite seleccionar solo las muestras y los metadatos que queramos. Requiere un argumento `dataType`, que creo que tengo que marcar como `"relative_abundance"`.
- Creo que esto puede ser interesante sobre todo en datasets donde no queremos todas las muestras (por ejemplo, si hay pacientes CRC o que tomen antibióticos). Tengo que revisar si se puede elegir un estudio en concreto, porque el ejemplo que viene selecciona muestras con unas características concretas sin especificar estudio.
### Alfa y beta diversidad
En cualquier caso, entiendo que el siguiente paso es **mirar alfa y beta diversidad.**
- La página de phyloseq tiene bastantes [tutoriales](https://joey711.github.io/phyloseq/plot_richness-examples.html) que pueden venir bien
- [Producing P values for alpha diversity data?](https://github.com/joey711/phyloseq/issues/704)
- Habrá que ver qué estadísticos tengo que usar...
---
Nuevos comentarios Laura
# Análisis en R :computer:
:::info
`Phyloseq` es un paquete de R que se usa mucho para análisis de microbioma, los datos de curated se pueden pasar fácilmente. Un objeto `Phyloseq` está formado por tres elementos necesarios:
**1. La tabla de taxonomía.
2. La tabla de metadata.
3. La tabla de abundancia.** Y además...
4. La información filogenética. No es fundamental, pero sin ella no podemos calcular la distancia Unifrac.
:::
Para pasar los datos de Metaphlan a phyloseq usamos la función metaphlantophyloseq:
```r=1
metaphlanToPhyloseq <- function(
tax,
metadat=NULL,
simplenames=TRUE,
roundtointeger=FALSE,
split="|"){
## tax is a matrix or data.frame with the table of taxonomic abundances, rows are taxa, columns are samples
## metadat is an optional data.frame of specimen metadata, rows are samples, columns are variables
## if simplenames=TRUE, use only the most detailed level of taxa names in the final object
## if roundtointeger=TRUE, values will be rounded to the nearest integer
xnames = rownames(tax)
shortnames = gsub(paste0(".+\\", split), "", xnames)
if(simplenames){
rownames(tax) = shortnames
}
if(roundtointeger){
tax = round(tax * 1e4)
}
x2 = strsplit(xnames, split=split, fixed=TRUE)
taxmat = matrix(NA, ncol=max(sapply(x2, length)), nrow=length(x2))
colnames(taxmat) = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species", "Strain")[1:ncol(taxmat)]
rownames(taxmat) = rownames(tax)
for (i in 1:nrow(taxmat)){
taxmat[i, 1:length(x2[[i]])] <- x2[[i]]
}
taxmat = gsub("[a-z]__", "", taxmat)
taxmat = phyloseq::tax_table(taxmat)
otutab = phyloseq::otu_table(tax, taxa_are_rows=TRUE)
if(is.null(metadat)){
res = phyloseq::phyloseq(taxmat, otutab)
}else{
res = phyloseq::phyloseq(taxmat, otutab, phyloseq::sample_data(metadat))
}
return(res)
}
```
:::warning
Te dejo mis scripts del último proyecto que he estado analizando por si te son de utilidad... Eso sí son un poco :spaghetti:
:::
## Convertir en Objeto phyloseq
* Metadata: Añadir metadata ojo con los Sample_names :eyes: que si no coinciden peta...
```r=1
DATA_plus <- read_csv("Data/Objetos/data_plus.csv")
DATAmfa.cluster <- read_csv("Data/Objetos/DATAmfa.cluster.csv")
DATA <- DATAmfa.cluster %>% select(ID,Cluster)
names(DATA) <- c("Sample_ID", "Cluster")
DATA <- DATA %>% left_join(DATA_plus)
DATA$Cluster <- factor(DATA$Cluster, levels = c("1", "2"), labels=c("1","2"))
row.names(DATA) <- DATA$Sample_ID
sample <- sample_data(DATA)
sample_names(sample) <- DATA$Sample_ID
```
>Los archivos de metadatos que tengo son:
>- Para las muestras de MetaHit:
> - `ena-metadata-selection.csv` tiene una selección de metadata del ENA, con una fila por muestra (los run_accession de las muestras concatenadas están actualizados)
> - `metadata-costea.csv` tiene la información de la tabla EV1 de Costea, importante porque tenemos ahí BMI, edad, género y diabetes.
> - Podemos hacer el join con `ena_sample_accession` o con `internal_sample_name`
>- Para las muestras de Karlsson en principio lo tenemos todo en `curatedMetagenomicData`
> Mi script para preparar metadatos:
```r=1
data.metahit.1 <- read.csv('./Datasets/ena-metadata-selection.csv')
data.metahit.2 <- read.csv('./Datasets/metadata-costea.csv')
data.metahit.2$sample_accession <- data.metahit.2$ena_sample_accession
# une en una sola tabla:
data.metahit <- data.metahit.1 %>% left_join(data.metahit.2, by = 'sample_accession')
# elimina columnas no interesantes:
data.metahit <- data.metahit %>%
select(-c(internal_sample_name, ena_study_accession,
ena_sample_accession, m_otu_diversity, inserts_post_qc,
inserts_aligned_to_spec_i))
```
:::warning
**Blanca:** en este punto añadí un par de columnas al dataframe para que sirvieran de etiquetas una vez construyamos el clasificador (y también para poder comparar los grupos más fácilmente cuando hagamos análisis metagenómicos). Tenía dudas de cómo era mejor hacerlo, así que al final lo que decidí fue:
- Añadir una columna `mu` con valores binarios, `1` indica **metabólicamente insano** y `0` indica sano. Con esta columna haríamos la clasificación en dos grupos.
- En la columna `group` tenemos cuatro posibles valores:
- `0` para no obesos metabólicamente sanos (**MHNO**)
- `1` para obesos metabólicamente sanos (**MHO**)
- `2` para no obesos metabólicamente insanos (**MUNO**)
- `3` para obesos metabólicamente insanos (**MUO**)
En este dataset de MetaHit, tendríamos solo 1's en la columna `mu`, mientras que en la columna `group` tendremos 2's y 3's.
:::
```r=1
# crea labels:
data.metahit$mu <- rep(1, dim(data.metahit)[1])
data.metahit$group <- rep(5, dim(data.metahit)[1])
data.metahit$group[data.metahit$bmi < 30] <- 2
data.metahit$group[data.metahit$bmi >= 30] <- 3
# termina de arreglar el dataframe:
data.metahit$group <- factor(data.metahit$group)
data.metahit$mu <- factor(data.metahit$mu)
row.names(data.metahit) <- data.metahit$run_accession
# y aqui ya entra phyloseq:
sample.metahit <- sample_data(data.metahit)
sample_names(sample.metahit) <- gsub('E', 'profiled_E',
data.metahit$run_accession)
```
* **Taxonomía:** Leemos la tabla de abundancia "merged_abundance_table.txt". Quitamos la primera línea y arreglamso SampleID para que coincida.
> Revisar que no haya errores con los números.
```r=1
Abundance <- read_csv("Data/Metagenomica/merged_abundance_table.txt")
Abundance <- column_to_rownames(Abundance, var="clade_name")
Abundance$NCBI_tax_id <- NULL
```
> Blanca:
```r=1
abundances.metahit<- read.table('metaphlan_phyloseq/merged_abundance_table.txt',
skip = 1, sep = '\t', header = TRUE)
rownames(abundances.metahit) <- abundances.metahit$clade_name
abundances.metahit <- abundances.metahit %>% select(-c(clade_name))
abundances.metahit$NCBI_tax_id <- NULL
```
* Cargar script Function metaphlantophyloseq
```r=1
phyloseqin= metaphlanToPhyloseq(Abundance, metadat = sample)
phyloseqin
```
* Añadir **árbol filogenético** para poder calcular distancia Unifrac
```r=1
#Arbol
random_tree=rtree(ntaxa(phyloseqin), rooted = T, tip.label = taxa_names(phyloseqin))
#Objeto phyloseq con arbol
physeq.tree <- merge_phyloseq(phyloseqin, random_tree)
```
:::success
Lo ideal sería en este punto unir con los resultados de `curatedmetagenomicData`
:::
## Análisis de diversidades:
:::info
El análisis de diversidad puedes hacerlo con `phyloseq` o con `microbiome`
:::
> En los datos generados con metaphlan :custard:
> Necesitamos transformar a conteos absolutos para las diversidades. Porque metaphlan sólo dá abundancias relatvas por millon.
Para eso multiplicamos la abundancia de cada taxón por un millón - conteos por millón.
Debido a cómo funciona MetaPhlAn, no es posible tener "read counts". Sin embargo, se pueden tener "pseudo" conteos multiplicando las abundancias relativas por una constante y redondeando al número entero más cercano.
```r=1
GP1 <- transform_sample_counts(physeq.tree, function(x) 1E6 * x)
```
* Alfa diversidad
```r=1
plot_richness(GP1) # phyloseq
rich <- richness(GP1) # microbiome
plot_richness(GP1, color = "Cluster", x = "Cluster",
measures = c("Chao1", "Simpson", "Shannon")) +
geom_boxplot(aes(fill = Cluster), alpha=.7) +
scale_color_manual(values = c("#B2ABD2", "#b2df8a")) +
scale_fill_manual(values = c("#B2ABD2", "#b2df8a")) +
theme_bw() # phyloseq
```
> Aquí me sale un aviso:
> ```
> Warning message:
> In estimate_richness(physeq, split = TRUE, measures = measures) :
> The data you have provided does not have any singletons. This is highly suspicious. Results of richness estimates (for example) are probably unreliable, or wrong, if you have already trimmed low-abundance taxa from the data.
> We recommended that you find the un-trimmed data and retry.
> ```
> De momento no le doy mucha importancia por lo que he visto googleando un poco ([aquí](https://forum.qiime2.org/t/singletons-and-diversity-richness-indices/2971/5) y [aquí](https://github.com/benjjneb/dada2/issues/214)).
* Transformamos datos composicionales con el paquete microbiome.
```r=1
pseq.compositional <- transform(physeq.tree, "compositional")
```
* Beta diversidad
```r=1
ordinated_taxa_unifrac = ordinate(physeq.tree, method="MDS", distance="unifrac",weighted=TRUE)
plot_ordination(physeq.tree, ordinated_taxa_unifrac, color="Cluster", title = "Unifrac MDS Analysis") + theme_bw()
```
---
# LEfSe
## Instalación
La instalación me ha dado algunos problemas. Al final:
- Está instalado en `estudiante@totoro`, en el entorno de conda `biobakery`
- Antes de correr la instalación, añadí los canales que puse al instalar Humann en Gandalf:
```bash=v
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --add channels biobakery
```
- Lo he instalado con el siguiente comando:
```bash=v
conda install -c biobakery lefse --verbose > lefse_installation.txt 2>&1
```
## Correr LEfSe
Antes de correrlo hay que **formatear los datos.** Esto está hecho con el script `phyloseq-to-lefse.R`.
Para comprobar que el formateo era correcto, he probado primero a correr el LEfSe en el [servidor de Galaxy](https://huttenhower.sph.harvard.edu/galaxy/):
- No ha encontrado nada interesante al usar enfermedad y obesidad como clase y subclase: `No differentially abundant features found in export/galaxy-central/database/files/003/275/dataset_3275306.dat` (con los parámetros por defecto, pidiendo p < 0.05 para Kruskal-Wallis y Wilcoxon y LDA > 2.0).
- Sí ha encontrado alguna feature usando solamente enfermedad o solamente obesidad como clase.
Una vez instalado y comprobado que el problema no está en los ficheros de entrada, lo he ejecutado en el ordenador:
```bash
lefse_format_input.py lefse.txt lefse1.in -c 1 -s 2 -u 3 -o 1000000
lefse_run.py lefse1.in lefse1.res
Number of significantly discriminative features: 0 ( 65 ) before internal wilcoxon
No features with significant differences between the two classes
Number of discriminative features with abs LDA score > 2.0 : 0
```
:::warning
Ojo: creo que el setting por defecto `-y 0` de `lefse_run.py` no es el mismo que en la página:
```bash
-y {0,1} (for multiclass tasks) set whether the test is performed in a
one-against-one ( 1 - more strict!) or in a one-against-all
setting ( 0 - less strict) (default 0)
```
:::
- Pidiendo LDA > 1.0 sale el mismo output (`lefse_run.py -l 1 lefse1.in lefse1_LDA1.res`)
- Con `-y 1` también (`lefse_run.py -y 1 lefse1.in lefse1_y1.res`)
- Con ambas opciones, lo mismo (`lefse_run.py -y 1 -l 1 lefse1.in lefse1_y1_LDA1.res`)
:::danger
Voy a revisar lo que hice para generar el archivo de entrada, no vaya a ser que falte algo o que falle alguna cosa...
:::
Hemos decidido **probar con una sola variable categórica de 4 clases:** control-sano, control-CRC, obeso-sano, obeso-CRC.
```bash=v
lefse_format_input.py lefse-one-class.txt lefse-one-class.in -c 1 -u 2 -o 1000000
lefse_run.py -y 1 lefse-one-class.in lefse-one-class.res
```
# Humann3 :woman_in_lotus_position:
::: info
:eye-in-speech-bubble: Aqui se recomienda para paired-ends concatenar cada read y mandar sólo un archivo fasta.
:::
## Instalación (Blanca)
[Instrucciones que he seguido](https://huttenhower.sph.harvard.edu/humann/)
Lo vamos a instalar en el entorno `metaphlan`:
```bash
conda activate metaphlan
# por si acaso:
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --add channels biobakery
# INSTALAMOS:
conda install humann -c biobakery
```
Para testear la instalación:
```bash
humann_test # OK
```
:::danger
Tuve bastantes problemas porque humann, diamond y bowtie2 estaban ya instalados en el ordenador, en la carpeta `/home/lmarcos/.local/bin` (fuera de entornos de conda). Tuve que borrar esas instalaciones porque los programas se ejecutaban ahí en vez de en el entorno `metaphlan`, y las versiones no cuadraban (la de diamond).
:::
- Databases
```bash
humann_databases --download chocophlan full /home/laura/Compubio
humann2_databases --download uniref uniref90_diamond $INSTALL_LOCATION
humann_databases --download uniref uniref90_diamond /home/laura/Compubio
```
> ```bash
> humann_databases --download chocophlan full humanndb
> humann_databases --download uniref uniref90_diamond humanndb
> ```
- Poner más cores
> Estos comandos van al archivo de configuración, cambian la subsección `threads` de la sección `run_modes` y lo cambian a 6.
```bash
# humann2_config --update run_modes threads 6
humann_config --update run_modes threads 6
```
:::warning
Para instalar en la workstation de Marco:
- Creé entorno `metaphlan` a partir del archivo .yml del entorno metaphlan en Gandalf :mage:
- Salía error al invocar bowtie2, así que hice esto `conda install tbb=2020.2` tal y como vi en el [foro de biobakery](https://forum.biobakery.org/t/workflow-conda-install-bowtie-2-issue-bowtie2-align-s-error-while-loading-shared-libraries-libtbb-so-2/1831)
- Además, he ajustado el número de cores a 4 (ya que en la otra workstation usamos 6/8, aquí usaremos 4/6)
humann_config --update run_modes threads 4
:::
- Concatenar fastqs
```bash
perl cat_files.pl *1.fastq
```
> Me falta el script `cat_files.pl`. He puesto este a correr (dentro del directorio `CLEAN_READS`, que tiene las lecturas tras el trimming y eliminación del host):
```perl
#!/usr/bin/perl
#dar de argumento solo los _1
foreach $ar (@ARGV)
{
# sacamos cada archivo:
@field = split (/_/, $ar);
$name = $ar; # name = nombre completo
$name =~ s/_1.fastq//; # name = nombre antes de _
# obtenemos el _2 correspondiente:
$r1 = $ar; # nombre completo de lectura con _1
$r2 = $ar; # lo mismo que r1 por ahora
$r2 =~ s/1.fastq/2.fastq/; # cambia _1 por _2
# tenemos cada read en $r1 y $r2 y el nombre base en $name
system("cat $r1 $r2 > $name.fastq; rm $r1 $r2; sed -i 's/ /_/' $name.fastq");
}
```
> Creo que ha funcionado bien. El tamaño total de la carpeta sigue siendo 850 GB, tenemos el número de archivos que deberíamos tener y todos tienen el formato de nombre que esperamos. Haciendo un `head` o un `tail` no parece que haya nada raro. Los tamaños de archivo individuales también cuadran (antes estábamos entre los 5,5 y los 7 GB y ahora entre los 11-15 GB).
:::success
- Con las muestras de T2D tardó unos 7 min por pareja
- Con las muestras de sanos tarda unos 4,5 min por pareja
:::
:::warning
- He probado el script con las primeras 16 líneas de dos parejas de archivos y parece que funciona :crossed_fingers:
- De todas formas, no me fío todavía, así que voy a copiar dos parejas de archivos (enteros) a `CLEAN_READS/prueba` por si acaso. He buscado dos que ocuparan "poco":
- `cp ../ERR414296_1.fastq .`
- `cp ../ERR414296_2.fastq .`
- `cp ../ERR4142578_1.fastq .`
- `cp ../ERR4142578_2.fastq .`
Miro antes el número de líneas en cada archivo por si acaso:
```bash
wc -l *fastq
114561776 ERR4142578_1.fastq
114561776 ERR4142578_2.fastq
108002052 ERR414296_1.fastq
108002052 ERR414296_2.fastq
445127656 total
```
- Pongo la prueba a correr en `screen -S blanca-cat` (**16:24 h, 11 nov**)
```bash
time perl cat_files.pl *1.fastq
real 12m56,155s
user 3m33,183s
sys 7m34,580s
```
Tarda unos **7 min por pareja** --> 7\*65 = 455 min = **7 h 35 min**
Chequeamos:
- Los números de líneas cuadran :+1:
```bash
wc -l *fastq
229123552 ERR4142578.fastq
216004104 ERR414296.fastq
445127656 total
```
- Las cabeceras cuadran :+1:
```bash=v
head ERR4142578.fastq
head ../ERR4142578_1.fastq
tail ERR4142578
tail ../ERR4142578_2.fastq
head ERR414296.fastq
head ../ERR414296_1.fastq
tail ERR414296.fastq
tail ../ERR414296_2.fastq
```
:::
## Ejecutar Humann
```bash
for f in *.fastq; do humann -i $f -o hmp_subset; done
```
>Probamos el siguiente comando en un directorio de prueba `CLEAN_READS/prueba` (con los ficheros ERR414295 y ERR414278, que ocupan "poco") primero:
>```bash!
>for f in *.fastq; do humann -i $f -o HUMANN --verbose; done > humannOUT-12-11-2021-1010.txt 2>&1
>```
:::warning
**Blanca**
En el github de Humann explican que se puede [saltar metaphlan](https://github.com/biobakery/humann#custom-taxonomic-profile) si pasamos un perfil taxonómico:
_\[...] To run HUMAnN 3.0 with the custom taxonomic profile ($FILE), use the option "`--taxonomic-profile $FILE`". This will bypass running MetaPhlAn, which creates a taxonomic profile, and instead will use the custom taxonomic profile provided \[...]_
:::
El tema es que los archivos taxonómicos son de la forma:
/home/lmarcos/blanca/costea_dnk/OUTPUT/profiles/profiled_{RUNACC}.txt
```bash!
for F in CLEAN_READS/*.fastq; do BASE=${F##*/}; TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt; echo $BASE $TAX; done
```
- `$BASE` tiene el nombre sin CLEAN_READS y `$TAX` tiene la dirección completa de cada perfil taxonómico:
```bash!
for F in CLEAN_READS/prueba/*.fastq; do BASE=${F##*/}; TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt; humann -i $F -o CLEAN_READS/prueba/HUMANN --verbose --taxonomic-profile $TAX; done > humannOUT-12-11-2021-11-30.txt 2>&1
```
:::danger
Ahora empecé a tener fallos de DIAMOND:
`CRITICAL ERROR: Please update diamond from version 0.8.22 to version 0.9.36`
No conseguí resolverlos ni instalando ni actualizando Diamond, ni tampoco probando en nuevos environments. **El problema era que usaba el diamond de `home/lmarcos/.local/bin` y no el del paquete.** La solución fue **eliminar esa instalación,** así como las de bowtie2 y humann, para que solo funcionaran las del environment.
Después probé con un par de muestras:
```bash=v
for F in CLEAN_READS/prueba/*.fastq; do BASE=${F##*/}; TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt; humann -i $F -o CLEAN_READS/prueba/HUMANN --verbose --taxonomic-profile $TAX; done > humannOUT-12-11-2021-15-30.txt 2>&1
```
Corrió hasta el día siguiente a las 5:21, **unas 7 h por muestra.**
**Los archivos temporales ocupan mucho espacio,** sobre todo por los índices de bowtie2. La opción `--remove-temp-output` permite borrar archivos intermedios, y según leo [aquí](https://github.com/biobakery/humann#4-intermediate-temp-output-files), los que tanto ocupan se borrarían :+1:.
- Los ficheros que se generan son:
- $SAMPLENAME_genefamilies.tsv, $SAMPLENAME_pathabundance.tsv y $SAMPLENAME_pathcoverage.tsv. Estos tres son la salida de HUMAnN.
- Diez archivos intermedios: $SN_bowtie2_aligned.sam, $SN_bowtie2_aligned.tsv, $SN_bowtie2_index*, $SN_bowtie2_unaligned.fa, $SN_custom_chocophlan_database.ffn, $SN_metaphlan_bowtie2.txt, $SN_metaphlan_bugs_list.tsv, \$SN_$TRANSLATEDALIGN_aligned.tsv, \$SN_\$TRANSLATEDALIGN_unaligned.fa y el log ($SN.log). De estos solo nos interesaría el log.
:::
:::success
Borro las dos muestras del directorio CLEAN_READS (quedan dentro de "prueba"), muevo HUMANN para arriba y borro las carpetas temporales. Pongo el comando completo en la sesión blanca-viernes:
```bash!
for F in CLEAN_READS/*.fastq; do BASE=${F##*/}; TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt; humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX --remove-temp-output; done > humannOUT-14-11-2021-15-06.txt 2>&1
```
**Lleva unas 7 horas por muestra, por lo que hemos interrumpido la ejecución por ahora para seguir una vez tengamos todas las tablas de taxonomía.**
:::
```bash!
for F in CLEAN_READS/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt
humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX
rm HUMANN/*bowtie2*
rm HUMANN/*aligned*
rm HUMANN/*phlan*
done > humannOUT-14-11-2021-15-06.txt 2>&1
```
:::warning
Lo voy a probar con unas pocas muestras primero:
```bash
mkdir prueba
cp ERR01110* prueba #ERR01105, ERR01107 y ERR01109.conc.fastq
```
Corremos en `blanca-metaphlan`:
```bash
for F in CLEAN_READS/prueba/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT_MPA/profiles/profiled_${BASE%.*}.txt
humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX
rm HUMANN/*bowtie2*
rm HUMANN/*aligned*
rm HUMANN/*phlan*
done > humannOUT-23-11-2021-17-03.txt 2>&1
```
Ha corrido en 7 h 45 min --> unas 2 h 35 por muestra! :confetti_ball:
- Borro las tres muestras de prueba del directorio `CLEAN_READS`, quedan solo en `CLEAN_READS/prueba` para ahorrar un poco de tiempo.
- Modifico el comando para que borre bien los archivos temporales:
```bash
for F in ../CLEAN_READS/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT_MPA/profiles/profiled_${BASE%.*}.txt;
echo HUMANN/${BASE%.*}_humann_temp
done
```
:::
Funciona! Con esto **mantenemos el log de cada muestra,** pero no los demás archivos temporales (son muy pesados). Vamos allá:
:::success
```bash
for F in CLEAN_READS/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT_MPA/profiles/profiled_${BASE%.*}.txt
humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX
rm HUMANN/${BASE%.*}_humann_temp/*bowtie2*
rm HUMANN/${BASE%.*}_humann_temp/*aligned*
rm HUMANN/${BASE%.*}_humann_temp/*phlan*
done > humannOUT-24-11-2021-09-36.txt 2>&1
```
:::
### Segunda parte Humann
- Hay que localizar los que **ya corrieron y tenemos los resultados:**
- `ls ../HUMANN/ | cut -d '_' -f1 | grep -v log | sort | uniq `
- Son estos:
- ERR4142301|ERR414232|ERR414233|ERR414234|ERR414235|ERR414236|ERR4142378|ERR414239|ERR414242|ERR414244|ERR414278|ERR414295
- En Gandalf: ERR4142301 y ERR4142378
- Creamos directorio old-humann y metemos ahí esas dos lecturas y los resultados
- ERR414295 y ERR414278 están missing (1 de 2)
- Resto:
- ERR414232|ERR414233|ERR414234|ERR414235|ERR414236|ERR414239|ERR414242|ERR414244 (**8** archivos de 32 --> procesaremos **24**)
- Todas en Marco, tenemos que organizarlo
**En worsktation Marco (sesión screen `blanca`)**
Tenemos 32 archivos, de los cuales procesaremos 24.
> El fichero `humannOUT-Marco-02-12-2021-15-05.txt` recoge el log de un primer intento que ha dado error. He solucionado (espero) el problema, que creo que es con bowtie2, y lo vuelvo a poner. (Ver cuadro amarillo en Instalación)
1. Seleccionamos los archivos:
```bash
already=$(ls CLEAN_READS/ | grep -E 'ERR414232|ERR414233|ERR414234|ERR414235|ERR414236|ERR414239|ERR414242|ERR414244')
cd CLEAN_READS
for file in $already; do mv $file ../already-humann; done
```
2. Ahora podemos poner esto a correr sin problema (espero), en screen `blanca`:
```bash
for F in CLEAN_READS/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT_MPA/profiles/profiled_${BASE%.*}.txt
humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX
rm HUMANN/${BASE%.*}_humann_temp/*bowtie2*
rm HUMANN/${BASE%.*}_humann_temp/*aligned*
rm HUMANN/${BASE%.*}_humann_temp/*phlan*
done > humannOUT-Marco-02-12-2021-17-15.txt 2>&1
```
**En Gandalf :mage:**
Aquí he movido 31 ficheros fastq, de los cuales correremos 29.
En screen `blanca`:
```bash
for F in CLEAN_READS/*.fastq; do
BASE=${F##*/}
TAX=OUTPUT/profiles/profiled_${BASE%.*}.txt
humann -i $F -o HUMANN --verbose --taxonomic-profile $TAX
rm HUMANN/${BASE%.*}_humann_temp/*bowtie2*
rm HUMANN/${BASE%.*}_humann_temp/*aligned*
rm HUMANN/${BASE%.*}_humann_temp/*phlan*
done > humannOUT-Gandalf-02-12-2021-16-55.txt 2>&1
```
:::warning
Faltan 2 ficheros fastq en CLEAN_READS. Los he identificado y he movido las lecturas crudas a RAW_READS, de momento no hacemos más con ellos porque sus resultados están en HUMANN, pero **hay que volver a obtener las lecturas limpias.**
:::
## Procesado post-Humann
- Movemos todos los archivos a un directorio común:
```bash=
# directorio: ~/blanca/costea-dnk/HUMANN
for file in $(find -type f | grep tsv); do cp $file ../../humann_output_files/; done
# directorio: ~/blanca/healthy-metahit/HUMANN
for file in $(find -type f | grep tsv); do cp $file ../../humann_output_files/; done
```
- Activamos screen `blanca-stamp`
- Activamos entorno `metaphlan`
- Unir las tablas
```bash=v
humann_join_tables -i hmp_subset -o hmp_subset_genefamilies.tsv --file_name genefamilies
humann_join_tables --input hmp_subset --output hmp_subset_pathcoverage.tsv --file_name pathcoverage
humann_join_tables --input hmp_subset --output hmp_subset_pathabundance.tsv --file_name pathabundance
```
```bash=+
humann_join_tables -i humann_output_files/ -o humann_merge/metahit_genefamilies.tsv --file_name genefamilies
humann_join_tables -i humann_output_files/ -o humann_merge/metahit_pathcoverage.tsv --file_name pathcoverage
humann_join_tables -i humann_output_files/ -o humann_merge/metahit_pathabundance.tsv --file_name pathabundance
```
- Transformar contajes
```bash
humann_renorm_table -i hmp_subset_genefamilies.tsv -o hmp_subset_genefamilies-cpm.tsv --units cpm --update-snames
```
```bash=+
cd humann_merge
humann_renorm_table -i metahit_genefamilies.tsv -o metahit_genefamilies-cpm.tsv --units cpm --update-snames
```
- Abundancia relativa
```bash
humann_renorm_table -i hmp_subset_genefamilies.tsv -o hmp_subset_genefamilies-relab.tsv --units relab
humann_renorm_table -i hmp_subset_pathabundance.tsv -o hmp_subset_pathabundance-relab.tsv --units relab
humann_renorm_table -i hmp_subset_pathcoverage.tsv -o hmp_subset_pathcoverage-relab.tsv --units relab
```
```bash=+
humann_renorm_table -i metahit_genefamilies.tsv -o metahit_genefamilies-relab.tsv --units relab
humann_renorm_table -i metahit_pathcoverage.tsv -o metahit_pathcoverage-relab.tsv --units relab
humann_renorm_table -i metahit_pathabundance.tsv -o metahit_pathabundance-relab.tsv --units relab
```
:::info
Todos estos pasos los he metido en un script, `humann_postprocess.sh`
:::
- Comprobar que sea 1 millón
```
grep -v "#" hmp_subset_genefamilies-cpm.tsv | grep -v "|" | cut -f2 | paste -s -d+ - | bc
```
```
grep -v "#" metahit_genefamilies-cpm.tsv | grep -v "|" | cut -f2 | paste -s -d+ - | bc
```
## Unir los datos de `curatedMetagenomicData`
- Lo que estoy haciendo está en `write_humann_from_curated.R`.
- No estoy consiguiendo recuperar los objetos `gene_family`, así que por lo pronto voy a trabajar con el resto.
- Lo que me está pareciendo es que la salida de Humann que pasamos a STAMP son las **cuentas,** por lo que habría que poner el argumento `counts = TRUE` al hacer `returnSamples` con los datos de curated.
- Luego bastaría con unir los datasets: las nuevas muestras se unen como **columnas**, y si aparecen rutas nuevas se añadirán como nuevas filas --> quizá lo más fácil sea unirlos con `dplyr::full_join`
:::warning
- [ ] Comprobar si podemos pasar abundancias relativas --> es necesario esto realmente??
- [ ] Cargar todo en R y hacer el join
- [ ] Corregir los nombres de archivo en el encabezado del fichero final (quitar terminación `_Abundance` o `.conc_Abundance` con un `sed`)
:::
He hecho el join con este script:
```r=
library(tidyverse)
library(janitor)
setwd('C:/Users/usuario/Desktop/TFM')
## Abundances:
feng.abundance <- read_tsv('./Humann_to_stamp/pathway_abundance_feng.tsv')
karl.abundance <- read_tsv('./Humann_to_stamp/pathway_abundance_karlsson.tsv')
mh.abundance <- read_tsv('./Humann_to_stamp/metahit_pathabundance-relab.tsv') %>%
rename(pathway = `# Pathway`)
all.abundances <- full_join(feng.abundance, karl.abundance,
by = 'pathway') %>%
full_join(mh.abundance)
write.table(all.abundances,
'./Humann_to_stamp/all_pathway_abundances-relab.tsv',
quote = FALSE, row.names = FALSE, sep = '\t')
## Coverage:
feng.coverage <- read_tsv('./Humann_to_stamp/pathway_coverage_feng.tsv')
karl.coverage <- read_tsv('./Humann_to_stamp/pathway_coverage_karlsson.tsv')
mh.coverage <- read_tsv('./Humann_to_stamp/metahit_pathcoverage.tsv') %>%
rename(pathway = `# Pathway`)
all.coverages <- full_join(feng.coverage, karl.coverage,
by = 'pathway') %>%
full_join(mh.coverage)
write.table(all.coverages,
'./Humann_to_stamp/all_pathway_coverages.tsv',
quote = FALSE, row.names = FALSE, sep = '\t')
```
Por si acaso, luego he editado el fichero en `vim` como explican [aquí](https://stackoverflow.com/questions/10208411/convert-windows-file-to-linux-file) por si hubiera problemas con los retornos de carro (son distintos en Windows, que es donde tengo R, y en Linux)
:::danger
Me confunde un poco lo de estar usando abundancias relativas o cuentas. Según [esto](https://www.biorxiv.org/content/10.1101/103085v2.full#:~:text=The%20two%20abundance%20output%20files%20were%20normalized%20in%20terms%20of%20relative%20abundance%20through%20the%20%E2%80%9Chumann2_renorm_table%E2%80%9D) los de curated pasaron todo a abundancias relativas, así que he usado eso para las `pathway_abundances`. Para los ficheros de `coverage` no hicieron la normalización, según eso, así que he usado las cuentas tal cual.
:::
## Humann to STAMP
```
humann_to_stamp.pl hmp_subset_pathabundance.tsv > hummann_path.spf
humann2_split_stratified_table -i hmp_subset_pathabundance.tsv -o split_folder
sed -s hmp_subset_pathabundance_unstratified.tsv > pathabundance_unstratified_head-fix.tsv
sed -s/# Pathway/MetaCyc_pathway/‘split_folder/pathabundance_example_unstratified.tsv > pathabundance_example_unstratified_head-fix.tsv
```
```bash=
humann_split_stratified_table -i all_pathway_abundances-relab.tsv -o split_folder
humann_split_stratified_table -i all_pathway_coverages.tsv -o split_folder_coverage
sed 's/pathway/MetaCyc_pathway/' split_folder/all_pathway_abundances-relab_unstratified.tsv > pathabundance_unstratified_head-fix.tsv
sed 's/pathway/MetaCyc_pathway/' split_folder_coverage/all_pathway_coverages_unstratified.tsv > pathcoverage_unstratified_head-fix.tsv
# fix headers:
sed -i 's/.conc_Abundance//g' pathabundance_unstratified_head-fix.tsv
sed -i 's/_Abundance//g' pathabundance_unstratified_head-fix.tsv
# change NAs to 0:
sed -i 's/NA/0/g' pathabundance_unstratified_head-fix.tsv
sed -i 's/.conc_Coverage//g' pathcoverage_unstratified_head-fix.tsv
sed -i 's/_Coverage//g' pathcoverage_unstratified_head-fix.tsv
sed -i 's/NA/0/g' pathcoverage_unstratified_head-fix.tsv
```
## STAMP
:::spoiler
Si al cargar los datos en STAMP sale un error del tipo "Key error: **''** at line 44 of file **Metadata.pyo**", cerrar y volver a abrir sin más! Por algún motivo se raya al cargar cosas en una sesión en la que ya hemos cargado datos.
:::
Cargando esto en STAMP, vemos que en la PCA las muestras se separan por estudio:
_Usando de entrada fichero de pathway coverage:_
Por lo que he probado a pasárselo a `MMUPHin::adjust_batch`, y parece que se soluciona:
Los ficheros que estoy usando por ahora son los de **pathway abundance** y **pathway coverage**, ya que los de **gene families** son muy pesados y no consigo recuperarlos de `curatedMetagenomicData` (para las muestras de Karlsson y Feng).
Como filtrando únicamente por p-valor ajustado salen gráficas con muchísimas rutas, estoy jugando también con los filtros de effect size:
### Pathway coverage
- Features con q-valor < 0.05 y effect size > 1 ("difference between proportions"):
- 
- Si filtro solo por q-valor, tengo que bajar mucho para reducir el número de features que aparecen. Lo mínimo que permite es $q < 10^{-5}$, con lo que obtenemos 49 features:
- 
:::warning
He probado con opciones en "Unclassified > Remove unclassified reads" y siguen apareciendo las UNMAPPED (secuencias no asignadas a ningún gen) y UNINTEGRATED (genes no asignados a ningún pathway).
:::
### Pathway abundance
- En este fichero no obtenemos features significativas si filtramos por q < 0.05 y effect size > 1 como en el anterior. Ni siquiera si reducimos a > 0.5 obtenemos nada.
- Filtrando solo por q-valor < 0.05 salen 206 features significativas, y bajando hasta $q < 10^{-5}$ quedan 83:
Sign in with Wallet
Connect another wallet