1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 1/21 Cuaderno de Google Colaboratory para análisis de secuencias de virus utilizando la tecnología de nanoporos Elaborado por: Alejandra Gil-Ordóñez. a.gil@cgiar.org Aprobado por: Wilmer Cuéllar. w.cuellar@cgiar.org Laboratorio de Virología y Protección de Cultivos. Programa de Yuca. Centro Internacional de Agricultura Tropical (CIAT). Palmira, Valle del Cauca. Colombia. Parte 1 Guppy basecaller En primer lugar, se debe realizar el llamado de bases. Este proceso consiste en asignar una base nucleotídica de acuerdo con el cambio de corriente registrado. El MinION tiene esta función incluída (fast basecalling), sin embargo es preferible realizar el llamado de bases en un cluster de cómputo, ya que el comando a ejecutar permite aumentar la calidad del llamado (algoritmo de mayor precisión) y realizar algunos ajustes adicionales como retirar barcodes y adaptadores. NOTA: Recordemos que esta tecnología se basa en nanoporos a través de los cuales circula una corriente iónica constante. Con ayuda de una proteína motora, la doble cadena de ADN (o híbrido ARN-ADN) se desenrrolla, y cada porción monocatenaria pasa a través del nanoporo impulsado por el gradiente energético. En este proceso, cada nucleótido produce un cambio de voltaje característico que es censado por el equipo y empleado para determinar la identidad nucleotídica a ~450 bases por segundo. Ejemplo de señal eléctrica de buena (A) y baja (B) consistencia mailto:a.gil@cgiar.org mailto:w.cuellar@cgiar.org 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 2/21 Este cuaderno describe el procesamiento de datos de secuenciación de Nanopore (archivos fast5) en un entorno de cuaderno interactivo de Google Colab. Esto es posible gracias al tiempo de ejecución habilitado para GPU que está disponible a través de Colab. Algunas cosas a tener en cuenta antes de continuar: Necesitará una cuenta en el foro de la comunidad ONT para descargar Guppy, así que asegúrese de tener una y poder acceder a la sección de descargas. Este es un enfoque basado en la nube, lo que signi�ca que todos los datos se ubicarán en alguna instancia de la nube en algún lugar (uso Google Drive en este ejemplo). Esto puede no ser apropiado para los datos que tiene Considere esto cuidadosamente antes de cargar cualquiera de sus datos. Este es actualmente un servicio gratuito, es posible que se elimine en cualquier momento; esto queda completamente a discreción de Google. Como parte de esto, Google tiene el derecho de monitorear el uso y puede acelerar o denegar la asignación de recursos a los usuarios que se ejecutan constantemente. La cantidad y el tipo de recursos asignados pueden y probablemente cambiarán. La instancia de GPU actual usa GPU que funcionan con Guppy, y el disco disponible es de aproximadamente 64 Gb, pero esto puede cambiar. La desconexión del tiempo de ejecución es una cosa, si la computadora portátil está inactiva demasiado tiempo, se desconectará. Es posible quedarse sin memoria/RAM (Se le recomienda inactivar la suspención de su equipo mientras realiza el procedimiento). Actualmente colab dispone de tres modelos principales de GPU, uno de ellos gratuito (T4) y otros dos disponibles con cuenta premium (A100 y T100). Guppy basecaller por GPU 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 3/21 Iniciar tiempo de ejecución de GPU Antes de comenzar es necesario con�gurar con GPU. Lo primero es asegurarse de que el tiempo de ejecución esté con�gurado para usar una GPU. Hacer esto es bastante simple: Ir al menú a Entorno de ejecución o Runtime Seleccione la opción Cambiar tipo de entorno de ejecución o Change runtime type Asegúrese de que Hardware accelerator esté con�gurado en GPU Una vez que lo anterior esté con�gurado, debería poder ejecutar el siguiente bloque de código. Si tiene éxito, debería ver algo /device:GPU:0 como la salida. Esto signi�ca que la GPU está disponible para su uso. Comprobar la presencia de una GPU /device:GPU:0' ' import tensorflow as tf tf.test.gpu_device_name() # Esto le dirá el número de dispositivo (debe ser 0 con una sola GPU) El siguiente código indicará el modelo de GPU asignado. Deberá tener disponible la GPU "Tesla T4". Tesla T4' ' import torch torch.cuda.get_device_name(0) # Esto le dirá el nombre/modelo de la GPU 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 4/21 También podemos ver algunas características adicionales de la GPU como también la versión del CUDA. !nvidia-smi !/usr/local/cuda/bin/nvcc --version # Esto le mostrará las características de la GPU y el Cuda Mon Nov 27 08:46:59 2023 +-----------------------------------------------------------------------------+ | NVIDIA-SMI 525.105.17 Driver Version: 525.105.17 CUDA Version: 12.0 | |-------------------------------+----------------------+----------------------+ | GPU Name Persistence-M| Bus-Id Disp.A | Volatile Uncorr. ECC | | Fan Temp Perf Pwr:Usage/Cap| Memory-Usage | GPU-Util Compute M. | | | | MIG M. | |===============================+======================+======================| | 0 Tesla T4 Off | 00000000:00:04.0 Off | 0 | | N/A 51C P0 28W / 70W | 839MiB / 15360MiB | 0% Default | | | | N/A | +-------------------------------+----------------------+----------------------+ +-----------------------------------------------------------------------------+ | Processes: | | GPU GI CI PID Type Process name GPU Memory | | ID ID Usage | |=============================================================================| +-----------------------------------------------------------------------------+ nvcc: NVIDIA (R) Cuda compiler driver Copyright (c) 2005-2022 NVIDIA Corporation Built on Wed_Sep_21_10:33:58_PDT_2022 Cuda compilation tools, release 11.8, V11.8.89 Build cuda_11.8.r11.8/compiler.31833905_0 Deberá tener acceso al foro de la comunidad ONT aquí para poder acceder a la sección de descargas y obtener una copia de Guppy. Una vez que tenga acceso y pueda navegar a la sección "Descargas de software" del foro de la comunidad ONT, verá una lista de Guppy. Recomiendo obtener los archivos binarios precompilados, es decir, la versión que aparece como GPU Linux x64-bit, debe tener un nombre de archivo similar a ont-guppy_X.X.X_linux64.tar.gz , donde las X indican el número de versión. Descargar Guppy https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fnanoporetech.com%2Fcommunity 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 5/21 Puede copiar el enlace a esta descarga y pegarlo en el bloque de código a continuación, es decir, reemplace la sección GuppyBinary . Ejecute el bloque de código y Guppy se descargará. %%shell GuppyBinary=https://cdn.oxfordnanoportal.com/software/analysis/ont-guppy_6.5.7_linux64.tar.gz wget $GuppyBinary --2023-11-27 08:47:06-- https://cdn.oxfordnanoportal.com/software/analysis/ont-guppy_6.5.7_linux64.tar Resolving cdn.oxfordnanoportal.com (cdn.oxfordnanoportal.com)... 52.85.247.23, 52.85.247.114, 52.85.247 Connecting to cdn.oxfordnanoportal.com (cdn.oxfordnanoportal.com)|52.85.247.23|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1865218294 (1.7G) [application/x-tar] Saving to: ‘ont-guppy_6.5.7_linux64.tar.gz’ ont-guppy_6.5.7_lin 100%[===================>] 1.74G 43.6MB/s in 42s 2023-11-27 08:47:49 (42.5 MB/s) - ‘ont-guppy_6.5.7_linux64.tar.gz’ saved [1865218294/1865218294] Antes de que podamos usar los binarios de Guppy, necesitamos extraer el archivo que descargamos. Modi�que en el bloque de código a continuación la versión que descargó y luego ejecute el bloque de código. Si usamos la Extraiga los binarios Guppy comprimidos https://cdn.oxfordnanoportal.com/software/analysis/ont-guppy_6.5.7_linux64.tar.gz 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 6/21 versión 6.4.8 como ejemplo: %%shell tar -xzvf ont-guppy_6.5.7_linux64.tar.gz ont-guppy/bin/ ont-guppy/bin/autoconfigure_guppy_server.py ont-guppy/bin/guppy_basecall_server ont-guppy/bin/guppy_basecaller ont-guppy/bin/guppy_basecaller_duplex ont-guppy/data/ ont-guppy/data/YHR174W.fasta ont-guppy/data/adapter_scaling_dna_r10.3_min.jsn ont-guppy/data/adapter_scaling_dna_r10.3_prom.jsn ont-guppy/data/adapter_scaling_dna_r10.4_e8.1.jsn ont-guppy/data/adapter_scaling_dna_r9.4.1_e8.1.jsn ont-guppy/data/adapter_scaling_dna_r9.4.1_min.jsn ont-guppy/data/adapter_scaling_dna_r9.4.1_prom.jsn ont-guppy/data/certs-bundle.crt ont-guppy/data/dna_r10.3_450bps_fast.cfg ont-guppy/data/dna_r10.3_450bps_fast_mk1c.cfg ont-guppy/data/dna_r10.3_450bps_fast_prom.cfg ont-guppy/data/dna_r10.3_450bps_hac.cfg ont-guppy/data/dna_r10.3_450bps_hac_mk1c.cfg ont-guppy/data/dna_r10.3_450bps_hac_prom.cfg ont-guppy/data/dna_r10.3_450bps_sup.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_fast.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_fast_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_fast_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_hac.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_hac_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_hac_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_fast.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_fast.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_fast_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_fast_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_hac.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_hac.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_hac_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_hac_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_sup.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_sup.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5hmc_5mc_cg_sup_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_fast.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_fast.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_fast_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_fast_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_hac.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_hac.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_hac_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_hac_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_sup.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_sup.jsn ont-guppy/data/dna_r10.4.1_e8.2_260bps_modbases_5mc_cg_sup_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_260bps_sup.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_fast.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_fast_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_fast_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_hac.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_hac_mk1c.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_hac_prom.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_modbases_5hmc_5mc_cg_fast.cfg ont-guppy/data/dna_r10.4.1_e8.2_400bps_5khz_modbases_5hmc_5mc_cg_fast.jsn t /d t /d 0 8 2 00b kh db h f t k f 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 7/21 Podremos veri�car el correcto funcionamiento del programa guppy_basecaller ejecutando alguna función directamente desde la carpeta con los binarios que hemos descargado. %%shell /content/ont-guppy/bin/guppy_basecaller --help : Guppy Basecalling Software, (C) Oxford Nanopore Technologies plc. Version 6.5.7+ca6d6af, minimap2 Use of this software is permitted solely under the terms of the end user license agreement (EULA). By running, copying or accessing this software, you are demonstrating your acceptance of the EULA. The EULA may be found in /content/ont-guppy/bin Usage: With config file: guppy_basecaller -i -s -c [options] With flowcell and kit name: guppy_basecaller -i -s --flowcell --kit List supported flowcells and kits: guppy_basecaller --print_workflows Use GPU for basecalling: guppy_basecaller -i -s -c --device [options] Command line parameters: --adapter_pt_range_scale Set polyT/adapter range scale for setting read scaling median absolute deviation. --as_cpu_threads_per_scaler Number of CPU worker threads per adapter scaler. --dmean_threshold Threshold for coarse stall event detection --dmean_win_size Window size for coarse stall event detection --as_gpu_runners_per_device Number of runners per GPU device for adapter scaling. --jump_threshold Threshold level for rna stall detection --max_search_len Maximum number of samples to search through for the stall --as_model_file Path to JSON model file for adapter scaling. --noisiest_section_scaling_max_size Threshold read size in samples under which nosiest-section scaling will be performed. --as_num_scalers Number of parallel scalers for adapter scaling. --override_scaling Manually provide scaling parameters rather than estimating them from each read. --pt_median_offset Set polyT median offset for setting read scaling median. --pt_minimum_read_start_index Set minimum index for read start sample required to attempt polyT scaling. --pt_required_adapter_drop Set minimum required current drop from adapter max to polyT detection. --pt_scaling Enable polyT/adapter max detection for read scaling. --as_reads_per_runner Maximum reads per runner for adapter scaling. --scaling_mad Median absolute deviation to use for manual scaling. --scaling_med Median current value to use for manual scaling. 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 8/21 %%shell /content/ont-guppy/bin/guppy_aligner --help guppy_aligner, part of Guppy basecalling suite, (C) Oxford Nanopore Technologies plc. Version 6.5.7+ Use of this software is permitted solely under the terms of the end user license agreement (EULA). By running, copying or accessing this software, you are demonstrating your acceptance of the EULA. The EULA may be found in /content/ont-guppy/bin Usage: guppy_aligner -i -s --align_ref --align_type Command line parameters: --bam_out Emit BAM files instead of SAM. --index Output BAM index file. -q [ --records_per_file ] Maximum number of records per output file, 0 means no maximum. --alignment_filtering Specify whether to filter reads based on their alignment status --align_type Specify whether you want full or coarse alignment. Valid values are (auto/full/coarse). --bed_file Path to .bed file containing areas of interest in reference genome. -a [ --align_ref ] Reference FASTA or index file. --minimap_opt_string Specify minimap2 options. See `guppy_aligner --minimap_opt_string --help` for details). -c [ --config ] Configuration file for application. -d [ --data_path ] Path to use for loading any data files the application requires. -i [ --input_path ] Path to input files. --progress_stats_frequency Frequency in seconds in which to report progress statistics, if supplied will replace the de -z [ --quiet ] Quiet mode. Nothing will be output to STDOUT if this option is set. -r [ --recursive ] Search for input file recursively. -s [ --save_path ] Path to save output files. -h [ --help ] Display the application usage help. -v [ --version ] Display the application version information. --trace_category_logs Enable trace logs - list of strings with the desired names. --trace_domains_config Configuration file containing list of trace domains to include in verbose logging (if enable --verbose_logs Enable verbose logs. -t [ --worker_threads ] Number of worker threads. --disable_pings Disable the transmission of telemetry pings. --ping segment duration 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 9/21 Al montar su Google Drive, podrá cargar archivos fast5 que se pueden procesar y la salida se puede volver a escribir en la misma ubicación dentro de Drive. El siguiente fragmento realiza el montaje. Se le pedirá que se autentique, solo siga las instrucciones y todo debería ir bien. Vincular con tu Google Drive from google.colab import drive drive.mount('/content/drive') Mounted at /content/drive Para este caso, vamos tener nuestros archivos .fast5 en una carpeta dentro de nuestro Drive. Es necesario que los archivos esten dentro de su unidad de Drive personal (no en compartidos), por lo tanto, podemos generar un acceso directo de nuestra carpeta del curso "Curso_Virus_Nanoporos" como se indica en las siguientes imagenes. En este paso, debemos ir al Drive a la carpeta compartida y crear un acceso directo al Drive, en la carpeta de Mi Unidad. Seleccionamos "Mi unidad" y damos click en "Añadir". 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 10/21 Podemos veri�car que la carpeta y los archivos montados estén vinculados a nuestra cuenta personal. !ls -lah /content/drive/MyDrive/curso_virus_2023 lrw------- 1 root root 0 Nov 27 08:52 /content/drive/MyDrive/curso_virus_2023 -> /content/drive/MyDrive ¡Se ve bien! Podemos ver una lista de archivos fast5. Con todo lo anterior funcionando, ahora podemos hacer el basecalling a nuestros datos. Primero estableceremos algunas variables. El siguiente bloque de código crea variables de shell para las ubicaciones de entrada y salida, el binario guppy (basecaller) y varios archivos modelo para la llamada base (Ej. fast, hac y bases modi�cadas). Basecalling con Guppy A continuación se explica en detalle cada argumento: -i . Ruta del input o archivos de entrada. Considera todos archivos en formato fast5 presentes en la ruta de trabajo que de�namos. Este argumento nos va a ser útil en las siguientes situaciones: Si se ha decidido realizar la secuenciación con fast base calling (llamado de bases de baja precisión), es recomendable realizar el llamado de bases de precisión sobre todas las carpetas fast5 generadas (pass, skip y fail) (la carpeta fast5_skip se genera cuando se detiene la secuenciación sin �nalizar el fast base calling). Si se debe realizar el llamado de bases sobre archivos recuperados de la carpeta temporales (tmp), debido a la interrupción abrupta de la secuenciación por pérdida de la fuente de energía (desconexión accidental del equipo o bajón de luz). -s guppy_out El save o carpeta de salida. Crea la carpeta donde se almacenan los resultados del llamado de bases. Con�guración 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 11/21 --recursive El argumento recursive o recursivo indica que se van a incluir todos los archivos fast5 presentes en carpetas dentro de carpetas de la ruta de trabajo en la que nos encontremos. Este argumento es útil cuando tenemos una estructura de trabajo similar al siguiente esquema: ├── fast5_pass │ ├── barcode01 │ │ ├── file1.fast5 │ │ ├── file2.fast5 │ │ └── file3.fast5 │ ├── barcode02 │ │ ├── file1.fast5 │ │ └── file2.fast5 │ └── barcode03 │ ├── file1.fast5 │ ├── file2.fast5 │ ├── file2.fast5 │ └── file4.fast5 └── En este ejemplo, se ha realizado la secuenciación con fast base calling y los archivos fast5 se encuentran demultiplexados dentro de la carpeta fast5_pass. NOTA: La demultiplexación es el proceso en el que se reconoce el barcode de cada secuencia y se agrupan las lecturas en carpetas por dicha etiqueta de identidad. Se lleva a cabo durante el llamado de bases: fast base calling y/o llamado de bases preciso por comando. En nuestro caso, mantenemos el argumento recursivo por defecto. --config dna_r9.4.1_450bps_hac.cfg El argumento con�g o con�guración indica el modelo de llamado de bases que queremos realizar. Es bastante recomendable mantener el modelo de alta precisión (dna_r9.4.1_450bps_hac.cfg). Alternativamente, se puede emplear el modelo fast base calling (dna_r9.4.1_450bps_fast.cfg) disponible en el equipo. Debemos considerar que aunque el modelo preciso es aproximadamente >10X (GPU) y >20X (CPU) más lento que el modelo rápido (fast base calling), el ajuste realizado con el modelo alta de precisión facilita los análisis corriente abajo. --device cuda:0 El argumento device o dispositivo indica que el proceso se va a realizar en núcleos de GPU paralelizados. Estos núcleos están ajustados para realizar una gran cantidad de cálculos simultáneamente, lo que optimiza el proceso reduciendo su tiempo de ejecución. --gpu_runners_per_device 15 Número de subprocesos GPU de trabajo a emplear. --num_callers 1 Número de basecallers paralelos a crear. --barcode_kits La referencia del kit de barcode usado, este va estar relacionado con el número de muestras. EXP- NBD104 con barcodes del 1-12, EXP-NBD114 contiene los barcodes del 13-24, mientras que EXP-NBD196 con los 96 barcodes. --trim_barcodes Una vez realizado el demultiplexación elimina los barcodes de cada lectura. 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 12/21 En el siguiente bloque de comandos puede ejecutar el guppy-basecaller . Recuerde dejar su computador sin hibernar ni apagar, este proceso tardará varios minutos. Ejecutar Guppy También podemos crear una carpeta que en esta ocasión llamaremos virus_hac en nuestro Drive donde guardemos nuestros resultados del guppy , esto será necesario para poder encontrar nuestros archivos más adelante. %%shell inputPath="/content/drive/MyDrive/curso_virus_2023/fast5/" outputPath="/content/drive/MyDrive/virus_hac/" guppy_bc=/content/ont-guppy/bin/guppy_basecaller # Configurar ruta de los bin guppy_cfg_fast=/content/ont-guppy/data/dna_r9.4.1_450bps_fast.cfg # fast model calling guppy_cfg_hac=/content/ont-guppy/data/dna_r9.4.1_450bps_hac.cfg # high accuracy calling guppy_cfg_mod=/content/ont-guppy/data/dna_r9.4.1_450bps_modbases_5mc_hac.cfg # base modification calling $guppy_bc -i $inputPath -s $outputPath \ --recursive \ --config $guppy_cfg_hac \ --gpu_runners_per_device 16 \ --cpu_threads_per_caller 2 \ --device cuda:0 \ --barcode_kits EXP-NBD196 ONT Guppy basecalling software version 6.5.7+ca6d6af, minimap2 version 2.24-r1122 config file: /content/ont-guppy/data/dna_r9.4.1_450bps_hac.cfg model file: /content/ont-guppy/data/template_r9.4.1_450bps_hac.jsn input path: /content/drive/MyDrive/curso_virus_2023/fast5/ save path: /content/drive/MyDrive/virus_hac/ chunk size: 2000 chunks per runner: 256 minimum qscore: 9 records per file: 4000 num basecallers: 4 gpu device: cuda:0 kernel path: runners per device: 16 Use of this software is permitted solely under the terms of the end user license agreement (EULA). By running, copying or accessing this software, you are demonstrating your acceptance of the EULA. The EULA may be found in /content/ont-guppy/bin Found 107 input read files to process. Init time: 2205 ms 0% 10 20 30 40 50 60 70 80 90 100% |----|----|----|----|----|----|----|----|----|----| *************************************************** Caller time: 1250975 ms, Samples called: 1987820059, samples/s: 1.58902e+06 Finishing up any open output files. Basecalling completed successfully. Ahora puede explorar sus resultados. !ls /content/drive/MyDrive/virus_hac/ 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 13/21 !head /content/drive/MyDrive/virus_hac/sequencing_summary.txt !ls /content/drive/MyDrive/virus_hac/pass !ls /content/drive/MyDrive/virus_hac/pass/barcode93 !ls /content/drive/MyDrive/virus_hac/pass/barcode93 | wc !ls -R /content/drive/MyDrive/virus_hac/pass PycoQC calcula métricas y genera grá�cos de control de calidad interactivos para datos de secuenciación de tecnologías Oxford Nanopore. PycoQC se basa en el archivo sequencing_summary.txt generado por Albacore y/o Guppy, pero si es necesario, también puede generar un archivo de resumen desde la base llamado archivos fast5. El paquete admite ejecuciones 1D y 1D2 generadas con dispositivos Minion, Gridion y Promethion, llamadas base con Albacore 1.2.1+ o Guppy 2.1.3+. PycoQC está escrito en Python3 puro. Python2 no es compatible. Para una introducción rápida, consulte el tutorial de Tim Kahlke disponible en https://timkahlke.github.io/LongRead_tutorials/QC_P.html La documentación completa está disponible en https://a-slide.github.io/pycoQC PycoQC pip install pycoQC Collecting pycoQC Downloading pycoQC-2.5.2.tar.gz (35 kB) Preparing metadata (setup.py) ... done Requirement already satisfied: numpy>=1.19 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (1. Requirement already satisfied: scipy>=1.5 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (1.1 Requirement already satisfied: pandas>=1.1 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (1. Collecting plotly==4.1.0 (from pycoQC) Downloading plotly-4.1.0-py2.py3-none-any.whl (7.1 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.1/7.1 MB 66.5 MB/s eta 0:00:00 Requirement already satisfied: jinja2>=2.10 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (3 Requirement already satisfied: h5py>=3.1 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (3.9. Requirement already satisfied: tqdm>=4.54 in /usr/local/lib/python3.10/dist-packages (from pycoQC) (4.6 Collecting pysam>=0.16 (from pycoQC) Downloading pysam-0.22.0-cp310-cp310-manylinux_2_28_x86_64.whl (21.9 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 21.9/21.9 MB 26.9 MB/s eta 0:00:00 Collecting retrying>=1.3.3 (from plotly==4.1.0->pycoQC) Downloading retrying-1.3.4-py3-none-any.whl (11 kB) Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from plotly==4.1.0->pyco Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from jinja2> Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1. Building wheels for collected packages: pycoQC Building wheel for pycoQC (setup.py) ... done Created wheel for pycoQC: filename=pycoQC-2.5.2-py3-none-any.whl size=39528 sha256=62f840e2ae76ce1871 Stored in directory: /root/.cache/pip/wheels/0a/d4/66/b00822fb375efe771b0a547850990f2930d0cc6ac2f6c78 Successfully built pycoQC Installing collected packages: retrying, pysam, plotly, pycoQC Attempting uninstall: plotly Found existing installation: plotly 5.15.0 Uninstalling plotly-5.15.0: https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Ftimkahlke.github.io%2FLongRead_tutorials%2FQC_P.html https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fa-slide.github.io%2FpycoQC 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 14/21 Successfully uninstalled plotly-5.15.0 ERROR: pip's dependency resolver does not currently take into account all the packages that are install lida 0.0.10 requires fastapi, which is not installed. lida 0.0.10 requires kaleido, which is not installed. lida 0.0.10 requires python-multipart, which is not installed. lida 0.0.10 requires uvicorn, which is not installed. cufflinks 0.17.3 requires plotly>=4.1.1, but you have plotly 4.1.0 which is incompatible. Successfully installed plotly-4.1.0 pycoQC-2.5.2 pysam-0.22.0 retrying-1.3.4 cd /content/drive/MyDrive/virus_hac/ /content/drive/MyDrive/virus_hac ls fail/ guppy_basecaller_log-2023-11-27_09-05-41.log guppy_basecaller_log-2023-11-27_08-53-55.log guppy_basecaller_log-2023-11-27_09-06-32.log guppy_basecaller_log-2023-11-27_08-55-56.log guppy_basecaller_log-2023-11-27_09-07-21.log guppy_basecaller_log-2023-11-27_08-56-44.log guppy_basecaller_log-2023-11-27_09-08-10.log guppy_basecaller_log-2023-11-27_08-57-33.log guppy_basecaller_log-2023-11-27_09-09-01.log guppy_basecaller_log-2023-11-27_08-58-23.log guppy_basecaller_log-2023-11-27_09-09-53.log guppy_basecaller_log-2023-11-27_08-59-10.log guppy_basecaller_log-2023-11-27_09-10-45.log guppy_basecaller_log-2023-11-27_09-00-01.log guppy_basecaller_log-2023-11-27_09-11-33.log guppy_basecaller_log-2023-11-27_09-00-50.log guppy_basecaller_log-2023-11-27_09-12-24.log guppy_basecaller_log-2023-11-27_09-01-37.log guppy_basecaller_log-2023-11-27_09-13-15.log guppy_basecaller_log-2023-11-27_09-02-24.log guppy_basecaller_log-2023-11-27_09-14-07.log guppy_basecaller_log-2023-11-27_09-03-11.log pass/ guppy_basecaller_log-2023-11-27_09-04-00.log sequencing_summary.txt guppy_basecaller_log-2023-11-27_09-04-50.log sequencing_telemetry.js !pycoQC --summary_file sequencing_summary.txt --html_outfile run_1.html Checking arguments values Check input data files Parse data files Merge data Cleaning data Discarding lines containing NA values 0 reads discarded Filtering out zero length reads 0 reads discarded Sorting run IDs by decreasing throughput Run-id order ['3c88f8fc73dec1ee0c7ed511637a7c399fa0fdf0'] Reordering runids Processing reads with Run_ID 3c88f8fc73dec1ee0c7ed511637a7c399fa0fdf0 / time offset: Cleaning up low frequency barcodes 1,286 reads with low frequency barcode unset Cast value to appropriate type Reindexing dataframe by read_ids 428,000 Final valid reads Loading plotting interface Found 428,000 total reads Found 332,761 pass reads (qual >= 7.0 and length >= 0) Generating HTML report Parsing html config file Running method run_summary Computing plot Running method basecall_summary Computing plot Running method alignment_summary No Alignment information available 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 15/21 Running method read_len_1D Computing plot Running method align_len_1D No Alignment information available Running method read_qual_1D Computing plot Running method identity_freq_1D No identity frequency information available Running method read_len_read_qual_2D Computing plot Running method read_len_align_len_2D No Alignment information available Running method align_len_identity_freq_2D No identity frequency information available Running method read_qual_identity_freq_2D No identity frequency information available Running method output_over_time Computing plot Running method read_len_over_time Computing plot Running method read_qual_over_time Computing plot Running method align_len_over_time No Alignment information available Running method identity_freq_over_time No identity frequency information available Running method barcode_counts C ti l t Usted puede descargar el archivo run_1.html y ver los resultados generados directamente en cualquier navegador, o verlo directamente en este cuaderno con el siguiente comando. import IPython IPython.display.HTML(filename='/content/drive/MyDrive/virus_hac/run_1.html') 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 16/21 PycoQC report Generated on 27/11/23 with pycoQC 2.5.2 General run summary Basecall summary Basecalled reads length All Reads Pass Status 18.03 18.03 Run Duration (h) 494 492 Active Channels 1 1 Number of Runids 16 16 Number of Barcodes All Reads Status 4.280000e+5 Reads 1.583856e+8 Bases 358 N50 276 Median Read Length 8.968 Median PHRED score PycoQC Docs GitHub Menu https://a-slide.github.io/pycoQC/ https://a-slide.github.io/pycoQC/pycoQC_usage.html https://github.com/a-slide/pycoQC 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 17/21 Basecalled reads PHRED quality Basecalled reads length vs reads 10% 193.00 25% 224.00 Median 277.00 75% 367.00 90% 508.00 100 2 5 1000 2 5 10k 2 5 100k 2 0 1000 2000 3000 4000 5000 Density 10%Aa 25%Aa MedianAa 75%Aa 90%Aa Basecalled length Re ad d en si ty All Reads Pass Reads 10% 5.34 25% 7.27 Median 8.97 75% 10.84 90% 12.58 10 20 30 0 500 1000 1500 2000 2500 3000 Density 10%Aa 25%Aa MedianAa 75%Aa 90%Aa Read quality scores Re ad d en si ty All Reads Pass Reads 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 18/21 g PHRED quality Output over experiment time 100 2 5 1000 2 5 10k 2 5 100k 2 5 10 15 20 25 30 Density Median 0 33.9 67.8 101.7 135.6 169.6 203.5 237.4 Basecalled length PH R ED q ua lit y sc or es All Reads Pass Reads 50% 8.88h 214,312 reads 75% 13.4h 320,965 reads 90% 16.22h 385,075 reads 99% 17.88h 424,088 re 100% 18.06h 427,999 re 100k 200k 300k 400k Cumulative Interval 50%Aa 75%Aa 90%Aa 99%Aa 100%Aa C ou nt All Reads Pass Reads All Bases 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 19/21 Read length over experiment time Read quality over experiment time 5 10 15 0 Experiment time (h) Pass Bases 5 10 15 100 2 5 1000 2 5 10k 2 5 100kMax Min 75% 25% Median Experiment time (h) A lig nm en t le ng th All Reads Pass Reads 10 15 20 Max Min 75% 25% Median n re ad P H R ED q ua lit y 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 20/21 Number of reads per barcode Channel activity over time 5 10 15 5 Experiment time (h) M ea n All Reads Pass Reads 15.4% 1.53% 7.67% 1.29% 12.2% 3.16% 8.1% 1.08% 2.6 % 16% 0.588% 0.528% 0.539% 5.74% 1.52% 22% barcode02 barcode03 barcode04 barcode05 barcode06 barcode07 barcode08 barcode09 barcode10 barcode12 barcode13 barcode14 barcode15 barcode16 barcode17 unclassified All Reads Pass Reads 14 16 18 40 45 50 ) 1/19/24, 9:48 AM Curso_virus_parte1_basecaller.ipynb - Colaboratory https://colab.research.google.com/drive/1A6W8uGM23mrXFbiRLEccZowkjN4RR64W?authuser=2#scrollTo=NYwW8ojgpgQD&printMode=true 21/21 Source summary files c 1 c 27 c 53 c 79 c 105 c 131 c 157 c 183 c 209 c 235 c 261 c 287 c 313 c 339 c 365 c 391 c 417 c 443 c 469 c 495 2 4 6 8 10 12 5 10 15 20 25 30 35 Channel id Ex pe ri m en t tim e (h ) All Reads Pass Reads All Bases Pass Bases Como ustedes han observado los resultadops del basecaller son archivos en formato FASTQ , lo cuales se han almacenado en un número variable de archivos dentro de disintas carpetas (una por barcode). Para continuar con nuestro análisis de beremos concatenar (unir) estos archivos para obtener uno solo en formato FASTQ. Para esto cat se implementa de la siguiente manera: Por ejemplo: cat barcode01/*.fastq > cat01.fastq Archivos fastq Primero creemos una carpeta donde almacenaremos los archivos concatenados. !mkdir /content/drive/MyDrive/virus_hac/concatenated/ Concatene los archivos .fastq del primer barcode. %%shell 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 1/40 Parte 2 Instalación de programas y paquetes Instación conda La intalación de conda mas kraken2, braken y krakentools tarda aproximadamente 3 minutos. pip install -q condacolab WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with import condacolab condacolab.install() ✨🍰✨ Everything looks OK! Para que reconozca el conda puede tardar unos segundos adicionales, espere que de clink en el siguiente comando para con�rmar que funcione de forma apropiada. !conda --version conda 23.1.0 0.0.1. Instalar programas para Mapeo #!conda install -c bioconda fastqc==0.11.9 -y !conda install -c bioconda minimap2 -y Collecting package metadata (current_repodata.json): done Solving environment: done ==> WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 ## Package Plan ## 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 2/40 environment location: /usr/local added / updated specs: - minimap2 The following packages will be downloaded: package | build ---------------------------|----------------- ca-certificates-2023.11.17 | hbcca054_0 151 KB conda-forge certifi-2023.11.17 | pyhd8ed1ab_0 155 KB conda-forge k8-0.2.5 | hdcf5f25_4 1.7 MB bioconda minimap2-2.26 | he4a0461_2 1.1 MB bioconda openssl-3.2.0 | hd590300_0 2.7 MB conda-forge zlib-1.2.13 | h166bdaf_4 92 KB conda-forge ------------------------------------------------------------ Total: 5.8 MB The following NEW packages will be INSTALLED: k8 bioconda/linux-64::k8-0.2.5-hdcf5f25_4 minimap2 bioconda/linux-64::minimap2-2.26-he4a0461_2 zlib conda-forge/linux-64::zlib-1.2.13-h166bdaf_4 The following packages will be UPDATED: ca-certificates 2022.12.7-ha878542_0 --> 2023.11.17-hbcca054_0 certifi 2022.12.7-pyhd8ed1ab_0 --> 2023.11.17-pyhd8ed1ab_0 openssl 3.1.0-h0b41bf4_0 --> 3.2.0-hd590300_0 Downloading and Extracting Packages k8-0.2.5 | 1.7 MB | : 0% 0/1 [00:00 WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 ## Package Plan ## environment location: /usr/local added / updated specs: 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 3/40 - qualimap The following packages will be downloaded: package | build ---------------------------|----------------- openjdk-8.0.382 | hd590300_0 88.4 MB conda-forge python_abi-3.10 | 4_cp310 6 KB conda-forge qualimap-2.2.2a | 1 25.8 MB bioconda ------------------------------------------------------------ Total: 114.2 MB The following NEW packages will be INSTALLED: openjdk conda-forge/linux-64::openjdk-8.0.382-hd590300_0 qualimap bioconda/linux-64::qualimap-2.2.2a-1 The following packages will be UPDATED: python_abi 3.10-3_cp310 --> 3.10-4_cp310 Downloading and Extracting Packages python_abi-3.10 | 6 KB | : 0% 0/1 [00:00 WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 ## Package Plan ## environment location: /usr/local added / updated specs: - samtools The following packages will be downloaded: package | build ---------------------------|----------------- curl-7.88.1 | hdc1c0ab_1 86 KB conda-forge samtools-1.6 | hc3601fc_10 501 KB bioconda ------------------------------------------------------------ Total: 587 KB The following NEW packages will be INSTALLED: curl conda-forge/linux-64::curl-7.88.1-hdc1c0ab_1 samtools bioconda/linux-64::samtools-1.6-hc3601fc_10 Downloading and Extracting Packages samtools-1.6 | 501 KB | : 0% 0/1 [00:00 WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 ## Package Plan ## environment location: /usr/local added / updated specs: - blast The following packages will be downloaded: package | build ---------------------------|----------------- blast-2.12.0 | hf3cf87c_4 22.2 MB bioconda entrez-direct-16.2 | he881be0_1 7.5 MB bioconda gettext-0.21.1 | h27087fc_0 4.1 MB conda-forge libidn2-2.3.4 | h166bdaf_0 157 KB conda-forge libunistring-0.9.10 | h7f98852_0 1.4 MB conda-forge pcre-8.45 | h9c3ff4c_0 253 KB conda-forge perl-5.32.1 | 4_hd590300_perl5 12.7 MB conda-forge perl-archive-tar-2.40 | pl5321hdfd78af_0 33 KB bioconda perl-carp-1.38 | pl5321hdfd78af_4 17 KB bioconda perl-common-sense-3.75 | pl5321hdfd78af_0 13 KB bioconda perl-compress-raw-bzip2-2.201| pl5321h87f3376_1 48 KB bioconda perl-compress-raw-zlib-2.105| pl5321h87f3376_0 75 KB bioconda perl-encode-3.19 | pl5321hec16e2b_1 2.1 MB bioconda perl-exporter-5.72 | pl5321hdfd78af_2 16 KB bioconda perl-exporter-tiny-1.002002| pl5321hdfd78af_0 25 KB bioconda perl-extutils-makemaker-7.70| pl5321hd8ed1ab_0 154 KB conda-forge perl-io-compress-2.201 | pl5321hdbdd923_2 84 KB bioconda perl-io-zlib-1.14 | pl5321hdfd78af_0 12 KB bioconda perl-json-4.10 | pl5321hdfd78af_0 56 KB bioconda perl-json-xs-2.34 | pl5321h4ac6f70_6 66 KB bioconda perl-list-moreutils-0.430 | pl5321hdfd78af_0 32 KB bioconda perl-list-moreutils-xs-0.430| pl5321h031d066_2 50 KB bioconda perl-parent-0.236 | pl5321hdfd78af_2 8 KB bioconda perl-pathtools-3.75 | pl5321hec16e2b_3 42 KB bioconda perl-scalar-list-utils-1.62| pl5321hec16e2b_1 44 KB bioconda perl-types-serialiser-1.01 | pl5321hdfd78af_0 13 KB bioconda wget-1.20.3 | ha35d2d1_1 815 KB conda-forge ------------------------------------------------------------ !conda install -c bioconda bbmap Collecting package metadata (current_repodata.json): done Solving environment: done ==> WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 6/40 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 # All requested packages already installed. 0.0.3. Instalación paquetes kraken2 Instale mediante conda los programas de kraken2 , krakentools y braken los cuales quedarán en un ambiente llamado tutorial. !conda create -yqn tutorial -c conda-forge -c bioconda kraken2 krakentools bracken Collecting package metadata (current_repodata.json): ...working... done Solving environment: ...working... done ## Package Plan ## environment location: /usr/local/envs/tutorial added / updated specs: - bracken - kraken2 - krakentools The following packages will be downloaded: package | build ---------------------------|----------------- biopython-1.81 | py310h2372a71_1 2.7 MB conda-forge blast-2.15.0 | pl5321h6f7f691_1 146.0 MB bioconda bracken-2.9 | py310h0dbaff4_0 396 KB bioconda bzip2-1.0.8 | hd590300_5 248 KB conda-forge c-ares-1.22.1 | hd590300_0 146 KB conda-forge curl-8.4.0 | hca28451_0 94 KB conda-forge icu-70.1 | h27087fc_0 13.5 MB conda-forge kraken2-2.1.3 | pl5321hdcf5f25_0 251 KB bioconda krakentools-1.2 | pyh5e36f6f_0 24 KB bioconda krb5-1.21.2 | h659d440_0 1.3 MB conda-forge libblas-3.9.0 |20_linux64_openblas 14 KB conda-forge libcblas-3.9.0 |20_linux64_openblas 14 KB conda-forge libcurl-8.4.0 | hca28451_0 377 KB conda-forge libgcc-ng-13.2.0 | h807b86a_3 755 KB conda-forge libgfortran-ng-13.2.0 | h69a702a_3 23 KB conda-forge libgfortran5-13.2.0 | ha4646dd_3 1.4 MB conda-forge libgomp-13.2.0 | h807b86a_3 412 KB conda-forge liblapack-3.9.0 |20_linux64_openblas 14 KB conda-forge libnghttp2-1.58.0 | h47da74e_0 617 KB conda-forge libnsl-2.0.1 | hd590300_0 33 KB conda-forge libopenblas-0.3.25 |pthreads_h413a1c8_0 5.3 MB conda-forge libsqlite-3.44.2 | h2797004_0 826 KB conda-forge libssh2-1.11.0 | h0841786_0 265 KB conda-forge 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 7/40 libstdcxx-ng-13.2.0 | h7e041cc_3 3.7 MB conda-forge libxml2-2.9.14 | h22db469_4 771 KB conda-forge libzlib-1.2.13 | hd590300_5 60 KB conda-forge ncbi-vdb-3.0.8 | hdbdd923_0 10.5 MB bioconda ncurses-6.4 | h59595ed_2 864 KB conda-forge numpy-1.26.2 | py310hb13e2d6_0 6.6 MB conda-forge ossuuid-1.6.2 | hf484d3e_1000 56 KB conda-forge perl-alien-build-2.48 | pl5321hec16e2b_0 172 KB bioconda perl-alien-libxml2-0.17 | pl5321hec16e2b_0 13 KB bioconda perl-business-isbn-3.007 | pl5321hd8ed1ab_0 18 KB conda-forge perl-business-isbn-data-20210112.006| pl5321hd8ed1ab_0 21 KB conda-forge perl-capture-tiny-0.48 | pl5321ha770c72_1 16 KB conda-forge perl-carp-1.50 | pl5321hd8ed1ab_0 22 KB conda-forge perl-common-sense-3.75 | pl5321hd8ed1ab_0 20 KB conda-forge perl-compress-raw-bzip2-2.201| pl5321h166bdaf_0 54 KB conda-forge perl-compress-raw-zlib-2.202| pl5321h166bdaf_0 83 KB conda-forge perl-constant-1.33 | pl5321hd8ed1ab_0 15 KB conda-forge perl encode 3 19 | pl5321h166bdaf 0 2 2 MB conda forge Especi�que la ruta de instalación del ambiente en el PATH . # in a shell you would run the following, however we want the binaries to be available in this notebook # source activate tutorial conda_path = !source activate tutorial && echo $CONDA_PREFIX:$PATH %set_env PATH=$conda_path env: PATH=['/usr/local/envs/tutorial:/usr/local/envs/tutorial/bin:/usr/local/condabin:/opt/bin:/usr/loc Con�rme que la instalación este correcta llamando kraken2 . !kraken2 --version Kraken version 2.1.3 Copyright 2013-2023, Derrick Wood (dwood@cs.jhu.edu) Instale scipy #@title **Instale scipy** !pip install scipy Collecting scipy Downloading scipy-1.11.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (60 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 60.4/60.4 kB 2.7 MB/s eta 0:00:00 Requirement already satisfied: numpy<1.28.0,>=1.21.6 in /usr/local/envs/tutorial/lib/python3.10/site-pa Downloading scipy-1.11.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (36.4 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 36.4/36.4 MB 34.7 MB/s eta 0:00:00 Installing collected packages: scipy Successfully installed scipy-1.11.4 WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with ##@title **Instale kraken-biom** !pip install kraken-biom mailto:dwood@cs.jhu.edu 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 8/40 Collecting kraken-biom Downloading kraken_biom-1.0.1-py2.py3-none-any.whl (12 kB) Collecting biom-format>=2.1.5 (from kraken-biom) Downloading biom_format-2.1.15-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 12.0/12.0 MB 76.3 MB/s eta 0:00:00 Collecting click (from biom-format>=2.1.5->kraken-biom) Downloading click-8.1.7-py3-none-any.whl.metadata (3.0 kB) Requirement already satisfied: numpy>=1.9.2 in /usr/local/envs/tutorial/lib/python3.10/site-packages (f Requirement already satisfied: scipy>=1.3.1 in /usr/local/envs/tutorial/lib/python3.10/site-packages (f Collecting pandas>=0.20.0 (from biom-format>=2.1.5->kraken-biom) Downloading pandas-2.1.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (18 kB) Collecting h5py (from biom-format>=2.1.5->kraken-biom) Downloading h5py-3.10.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (2.5 kB) Collecting python-dateutil>=2.8.2 (from pandas>=0.20.0->biom-format>=2.1.5->kraken-biom) Downloading python_dateutil-2.8.2-py2.py3-none-any.whl (247 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 247.7/247.7 kB 19.0 MB/s eta 0:00:00 Collecting pytz>=2020.1 (from pandas>=0.20.0->biom-format>=2.1.5->kraken-biom) Downloading pytz-2023.3.post1-py2.py3-none-any.whl.metadata (22 kB) Collecting tzdata>=2022.1 (from pandas>=0.20.0->biom-format>=2.1.5->kraken-biom) Downloading tzdata-2023.3-py2.py3-none-any.whl (341 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 341.8/341.8 kB 23.9 MB/s eta 0:00:00 Collecting six>=1.5 (from python-dateutil>=2.8.2->pandas>=0.20.0->biom-format>=2.1.5->kraken-biom) Downloading six-1.16.0-py2.py3-none-any.whl (11 kB) Downloading pandas-2.1.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (12.3 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 12.3/12.3 MB 91.4 MB/s eta 0:00:00 Downloading click-8.1.7-py3-none-any.whl (97 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 97.9/97.9 kB 7.7 MB/s eta 0:00:00 Downloading h5py-3.10.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (4.8 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.8/4.8 MB 90.3 MB/s eta 0:00:00 Downloading pytz-2023.3.post1-py2.py3-none-any.whl (502 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 502.5/502.5 kB 28.1 MB/s eta 0:00:00 Installing collected packages: pytz, tzdata, six, h5py, click, python-dateutil, pandas, biom-format, kr Successfully installed biom-format-2.1.15 click-8.1.7 h5py-3.10.0 kraken-biom-1.0.1 pandas-2.1.3 python WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with Instale Krona #@title **Instale Krona** !conda install -c bioconda krona -y Collecting package metadata (current_repodata.json): done Solving environment: done ==> WARNING: A newer version of conda exists. <== current version: 23.1.0 latest version: 23.10.0 Please update conda by running $ conda update -n base -c conda-forge conda Or to minimize the number of packages updated during conda update use conda install conda=23.10.0 ## Package Plan ## environment location: /usr/local 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 9/40 added / updated specs: - krona The following packages will be downloaded: package | build ---------------------------|----------------- krona-2.8.1 | pl5321hdfd78af_1 155 KB bioconda ------------------------------------------------------------ Total: 155 KB The following NEW packages will be INSTALLED: krona bioconda/noarch::krona-2.8.1-pl5321hdfd78af_1 Downloading and Extracting Packages Preparing transaction: done Verifying transaction: done Executing transaction: / Krona installed. You still need to manually update the taxonomy databases before Krona can generate taxonomic reports. The update script is ktUpdateTaxonomy.sh. The default location for storing taxonomic databases is /usr/local/opt/krona/taxonomy If you would like the taxonomic data stored elsewhere, simply replace this directory with a symlink. For example: rm -rf /usr/local/opt/krona/taxonomy mkdir /path/on/big/disk/taxonomy ln -s /path/on/big/disk/taxonomy /usr/local/opt/krona/taxonomy ktUpdateTaxonomy.sh done 1. Monte sus datos 1.1. Vincule con google drive from google.colab import drive drive.mount('/content/drive') Mounted at /content/drive 1.2. Descargue las muestras 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 10/40 !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/008/SRR13697208/SRR13697208_1.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/008/SRR13697208/SRR13697208_2.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/054/SRR13697154/SRR13697154_1.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/054/SRR13697154/SRR13697154_2.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/036/SRR13697136/SRR13697136_1.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/036/SRR13697136/SRR13697136_2.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/069/SRR13697069/SRR13697069_1.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/069/SRR13697069/SRR13697069_2.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/067/SRR13697067/SRR13697067_1.fastq.gz !wget -qnc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR136/067/SRR13697067/SRR13697067_2.fastq.gz for sample in ["SRR13697208", "SRR13697154", "SRR13697136", "SRR13697069", "SRR13697067"]: !gunzip {sample}_1.fastq.gz !gunzip {sample}_2.fastq.gz 1. Mapeando con Minimap2 Realizaremos un mapeo de las lecturas para determinar en qué parte del genoma se originaron nuestras lecturas. Hay una serie de herramientas para elegir y, aunque no existe un estándar de oro, existen algunas herramientas que son más adecuadas para análisis NGS particulares. Usaremos el programa Minimap2 , el cual es un programa rápido y versátil de alineación de propósito general para mapear secuencias de ADN y ARNm largas, muy adecuado para mapear lecturas de Oxford Nanopore y PacBio contra una secuencia de referencia. Se puede utilizar para: Mapeo de lecturas cortas precisas (preferiblemente más largas que 100 bases). Mapeo de lecturas genómicas de 1 kb con una tasa de error del 15 % (por ejemplo, lecturas genómicas de PacBio o Oxford Nanopore). Mapeo de lecturas directas ruidosas de ARN o ADNc de longitud completa. Mapeo y comparación de contigs de ensamblaje o cromosomas completos estrechamente relacionados de cientos de megabases de longitud. Obtener secuencia de referencia 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 11/40 Debido a que realizaremos el mapeo contra un secuencia de referencia debemos proporcionarla para poder correr el programa. Podemos usar un repositorio de secuencias del NCBI para descargarla, o proporcionar referencias propias. En este caso vaya a la página del NCBI y realice la búsqueda y descarga de los siguiente accesiones: AH015299.2 NC_001658.1 KC505249.1 KC505250.1 KC505251.1 KC505252.1 Descargue los archivos en formato .fastq https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fwww.ncbi.nlm.nih.gov%2F 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 12/40 Una vez descargado arrastre el archivo "sequence.fasta", suba el archivo al directorio /content/drive/MyDrive/virus_tax/minimap2/ . Este archivo servirá como sus genomas de referencia. Con el comando head o less podremos echar un vistazo a este genoma de referencia. !head /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta >AH015299.2 Cassava frogskin virus segment 4 putative RNA dependent RNA polymerase genes, partial cds CGCCCAGGCTATCACGAAAAACTCAACCAAGAACAAAGAAACCTCAACAACTACGACCTAATTTTCAAGC TACTGAAAAGTGGAAAGGTATTATCCCAAGATGAGCATTTGGCCGGGAGAACAATAATAGGGACAGATTT AACGCATATTTGGAAGGAACCAGAAATTTTAGTTCGGCCACCAGCGTATGACCTCTCAAAATCATACTAT ACTTTCTCACAGAGCAGTTTAGGATTTCAGGATCATAGTCATCAATATTTGGTTTCAGTTGATGTGAAGG ATATTCTAAAGCATTTCTTCTTACTGAAGGAGTATGAAAAAGGTAAAATCGAGGCGACGTATTCAGAACA GCAAAGAATAGATATCTTAATTCTATTTTTCTTAATTACTGAGGCTAACTATCAGACAAGATTCAGCAAT ATTTGGAGACAGTCATTGGCATTTTTATTATCCTATCTCGCAAAATACTGTGATCAAAAAGGAATTGACC CGAAATGTTTGGAACAAATCAACGAGATATATCCGACATTAATGCCTTTCCGATTAATGGCGACGCTTTT GAAGAGTTATCAAGAGATTCCTTTTATGGAGAAAGCGAATATGATCCACTGGCGAGAGAATACGAGTGGT !grep ">" /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta >AH015299.2 Cassava frogskin virus segment 4 putative RNA dependent RNA polymerase genes, partial cds >NC_001658.1 Cassava common mosaic virus, complete genome >KC505249.1 Cassava Polero-like virus isolate CM5460-10 ORF1 gene, partial cds; ORF2 gene, complete cds >KC505250.1 Cassava Torrado-like virus isolate CM5460-10 segment RNA1 polyprotein gene, partial cds >KC505251.1 Cassava Torrado-like virus isolate CM5460-10 segment RNA2 polyprotein gene, partial cds >KC505252.1 Cassava Colombian symptomless virus isolate CM5460-10 replicase, TGB1, TGB2, and coat prote >MN385565.1 Sinsheimervirus phiX174, complete genome >OR753134.1 Human mastadenovirus D isolate HAdV-D/USA/84282/2010[P101H26F44], complete genome >NC_029119.1 Staphylococcus phage SPbeta-like, complete genome 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 13/40 Minimap2 es un programa de alineación de propósito general para mapear el ADN o secuencias largas de ARNm contra una gran base de datos de referencia. Es muy bueno con genomas de referencia. Minimap2 Primero inspecionemos las opciones que disponemos con minimap . !minimap2 --help Usage: minimap2 [options] | [query.fa] [...] Options: Indexing: -H use homopolymer-compressed k-mer (preferrable for PacBio) -k INT k-mer size (no larger than 28) [15] -w INT minimizer window size [10] -I NUM split index for every ~NUM input bases [8G] -d FILE dump index to FILE [] Mapping: -f FLOAT filter out top FLOAT fraction of repetitive minimizers [0.0002] -g NUM stop chain enlongation if there are no minimizers in INT-bp [5000] -G NUM max intron length (effective with -xsplice; changing -r) [200k] -F NUM max fragment length (effective with -xsr or in the fragment mode) [800] -r NUM[,NUM] chaining/alignment bandwidth and long-join bandwidth [500,20000] -n INT minimal number of minimizers on a chain [3] -m INT minimal chaining score (matching bases minus log gap penalty) [40] -X skip self and dual mappings (for the all-vs-all mode) -p FLOAT min secondary-to-primary score ratio [0.8] -N INT retain at most INT secondary alignments [5] Alignment: -A INT matching score [2] -B INT mismatch penalty (larger value for lower divergence) [4] -O INT[,INT] gap open penalty [4,24] -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1] -z INT[,INT] Z-drop score and inversion Z-drop score [400,200] -s INT minimal peak DP alignment score [80] -u CHAR how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n] -J INT splice mode. 0: original minimap2 model; 1: miniprot model [1] Input/Output: -a output in the SAM format (PAF by default) -o FILE output alignments to FILE [stdout] -L write CIGAR with >65535 ops at the CG tag -R STR SAM read group line in a format like '@RG\tID:foo\tSM:bar' [] -c output CIGAR in PAF --cs[=STR] output the cs tag; STR is 'short' (if absent) or 'long' [none] --MD output the MD tag --eqx write =/X CIGAR operators -Y use soft clipping for supplementary alignments -t INT number of threads [3] -K NUM minibatch size for mapping [500M] --version show version number Preset: -x STR preset (always applied before other options; see minimap2.1 for details) [] - map-pb/map-ont - PacBio CLR/Nanopore vs reference mapping - map-hifi - PacBio HiFi reads vs reference mapping - ava-pb/ava-ont - PacBio/Nanopore read overlap - asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5% sequence divergence - splice/splice:hq - long-read/Pacbio-CCS spliced alignment 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 14/40 - sr - genomic short-read mapping See `man ./minimap2.1' for detailed description of these and other advanced command-line options. %%shell for i in {01..09}; do refDB="/content/drive/MyDrive/virus_tax/minimap2/sequence.fasta" minimap2 -ax map-ont $refDB /content/drive/MyDrive/virus_hac/concatenated/cat${i}.fastq > /content/drive/ done %%shell for i in {10..16}; do refDB="/content/drive/MyDrive/virus_tax/minimap2/sequence.fasta" minimap2 -ax map-ont $refDB /content/drive/MyDrive/virus_hac/concatenated/cat${i}.fastq > /content/drive/ done !minimap2 -ax map-ont sequence.fasta SRR13697136_1.fastq > SRR13697069_ref.sam [M::mm_idx_gen::0.019*0.73] collected minimizers [M::mm_idx_gen::0.030*0.84] sorted minimizers [M::main::0.030*0.84] loaded/built the index for 9 target sequence(s) [M::mm_mapopt_update::0.031*0.85] mid_occ = 10 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 9 [M::mm_idx_stat::0.032*0.85] distinct minimizers: 34980 (98.69% are singletons); average occurrences: 1 [M::worker_pipeline::12.408*1.48] mapped 1602714 sequences [M::main] Version: 2.26-r1175 [M::main] CMD: minimap2 -ax map-ont sequence.fasta SRR13697136_1.fastq [M::main] Real time: 12.413 sec; CPU: 18.372 sec; Peak RSS: 0.631 GB El archivo SAM es un archivo de texto delimitado por tabulaciones que contiene información para cada lectura individual y su alineación con el genoma. La versión binaria comprimida de SAM se denomina archivo BAM. Usamos esta versión para reducir el tamaño y permitir la indexación, lo que permite un acceso aleatorio e�ciente de los datos contenidos en el archivo. El archivo comienza con un encabezado, que es opcional. El encabezado se usa para describir la fuente de datos, secuencia de referencia, método de alineación, etc., esto cambiará dependiendo del alineador que se use. Después del encabezado está la sección de alineación. Cada línea que sigue corresponde a la información de alineación para una sola lectura. Cada línea de alineación tiene 11 campos obligatorios para información de mapeo esencial y un número variable de otros campos para información especí�ca del alineador. A continuación se muestra una entrada de ejemplo de un archivo SAM con los diferentes campos resaltados. Revisión del mapeo con samtools 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 15/40 Convertiremos el archivo SAM a formato BAM usando el programa samtools usando el comando view y le diremos a este comando que la entrada está en formato SAM (-S) y que salga en formato BAM (-b). %%shell for i in {01..09}; do samtools view -S -b /content/drive/MyDrive/virus_tax/minimap2/cat${i}.sam > /content/drive/MyDrive/virus_ samtools sort /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam -o /content/drive/MyDrive/virus_tax/m samtools index /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam done %%shell for i in {10..16}; do samtools view -S -b /content/drive/MyDrive/virus_tax/minimap2/cat${i}.sam > /content/drive/MyDrive/virus_ samtools sort /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam -o /content/drive/MyDrive/virus_tax/m samtools index /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam done %%shell samtools view -S -b SRR13697069_ref.sam > SRR13697069_ref.bam samtools sort SRR13697069_ref.bam -o SRR13697069_ref.sorted.bam samtools index SRR13697069_ref.sorted.bam Los archivos SAM / BAM se pueden clasi�car de varias maneras, por ejemplo, por ubicación de alineación en el cromosoma, por nombre de lectura, etc. Es importante tener en cuenta que las diferentes herramientas de alineación generarán SAM / BAM ordenadas de manera diferente, y las diferentes herramientas posteriores requieren de manera diferente archivos de alineación ordenados como entrada. También puede utilizar samtools para obtener más información sobre este archivo bam. 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 16/40 %%shell for i in {02..09}; do samtools flagstat /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam done %%shell for i in {10..16}; do samtools flagstat /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam done !samtools flagstat SRR13697069_ref.sorted.bam 1604194 + 0 in total (QC-passed reads + QC-failed reads) 1359 + 0 secondary 121 + 0 supplementary 0 + 0 duplicates 59499 + 0 mapped (3.71% : N/A) 0 + 0 paired in sequencing 0 + 0 read1 0 + 0 read2 0 + 0 properly paired (N/A : N/A) 0 + 0 with itself and mate mapped 0 + 0 singletons (N/A : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5) La NGS ha proporcionado una variedad de nuevos desafíos signi�cativos para la visualización de ensamblajes de secuencias. Estos incluyen el gran volumen de datos que se generan, longitudes de lectura variables y diferentes tipos de datos y formatos de datos asociados con la diversidad de nuevas tecnologías de secuenciación. La visualización de ensamblajes y asignaciones de lectura, juega un papel importante en el análisis de estos datos. Entre las utilidades el aseguramiento de la calidad y el descubrimiento cientí�co, a través de funciones como descripciones generales de cobertura de referencias completas, resaltado de variantes, marcado de lectura emparejada, pistas de características basadas en GFF3 y proteínas. Podemos hacer uso de visualizadores como Tablet , IGV o pavian (usaremos este último). Para esto deberemos descargar desde el drive el archivo catNN.bam y el catNN.bam.bai a visualizar. Visualización del mapeo Tablet es un visor grá�co liviano y de alto rendimiento para ensamblajes y alineaciones de secuencias de próxima generación. Puede descargar le programa desde la página o�cial: https://ics.hutton.ac.uk/tablet/ Tablet https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Ffbreitwieser.shinyapps.io%2Fpavian%2F https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fics.hutton.ac.uk%2Ftablet%2F 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 17/40 Integrative Genomics Viewer (IGV) es una herramienta interactiva de alto rendimiento y fácil de usar para la exploración visual de datos genómicos. Admite la integración �exible de todos los tipos comunes de datos y metadatos genómicos, generados por investigadores o disponibles públicamente, cargados desde fuentes locales o en la nube. IGV 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 18/40 Pavian es una aplicación de navegador interactiva para analizar y visualizar resultados de clasi�cación metagenómica de clasi�cadores como Kraken, KrakenUniq, Kraken 2, Centrifuge y MetaPhlAn. Pavian también proporciona un visor de alineación para validar coincidencias con un genoma particular. Pueden probar Pavian en https://fbreitwieser.shinyapps.io/pavian/. Pavian Ingrese a la página de pavian y suba los archivos .bam y .bam.bai al mismo tiempo en la pestaña Alignment viewer . https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Ffbreitwieser.shinyapps.io%2Fpavian%2F 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 19/40 Ejemplo de visualización de genoma. QC del mapeo 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 20/40 Qualimap (profundidad) La herramienta Qualimap toma archivos BAM como entrada y explora las características de las lecturas asignadas y sus propiedades genómicas. Proporciona una vista general de la calidad de los datos en un formato de informe HTML. Los informes son completos con estadísticas de mapeo resumidas junto con cifras útiles y archivos delimitados por tabulaciones de datos de métricas. Ejecutaremos Qualimap desde nuestro archivo .bam y luego analizaremos detalladamente el informe. %%shell for i in {02..09}; do qualimap bamqc -bam /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam -c -nw 400 -hm 3 --java-mem-siz done %%shell for i in {10..16}; do qualimap bamqc -bam /content/drive/MyDrive/virus_tax/minimap2/cat${i}.bam -c -nw 400 -hm 3 --java-mem-siz done !qualimap bamqc -bam /content/SRR13697069_ref.sorted.bam -c -nw 400 -hm 3 --java-mem-size=4G -nt 3 Display: Java memory size is set to 4G Launching application... QualiMap v.2.2.2-dev Built on 2016-12-11 14:41 Selected tool: bamqc Available memory (Mb): 32 Max memory (Mb): 3817 Starting bam qc.... Loading sam header... Loading locator... Loading reference... Number of windows: 400, effective number of windows: 408 Chunk of reads size: 1000 Number of threads: 3 Processed 50 out of 408 windows... Processed 100 out of 408 windows... Processed 150 out of 408 windows... 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 21/40 Processed 200 out of 408 windows... Processed 250 out of 408 windows... Processed 300 out of 408 windows... Processed 350 out of 408 windows... Processed 400 out of 408 windows... Total processed windows:408 Number of reads: 1602835 Number of valid reads: 58140 Number of correct strand reads:0 Inside of regions... Num mapped reads: 58140 Num mapped first of pair: 0 Num mapped second of pair: 0 Num singletons: 0 Time taken to analyze reads: 5 Computing descriptors... numberOfMappedBases: 7070595 referenceSize: 189443 numberOfSequencedBases: 7070513 numberOfAs: 1700689 Computing per chromosome statistics... Computing histograms... Overall analysis time: 5 end of bam qc Computing report... Writing HTML report... HTML report created successfully Finished import IPython IPython.display.HTML(filename='/content/drive/MyDrive/virus_tax/minimap2/cat07_stats/qualimapReport.html') import IPython IPython.display.HTML(filename='/content/drive/MyDrive/virus_tax/SRR13697069_ref.sorted_stats/qualimapReport 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 22/40 Logo Qualimap Report: BAM QC Input data and parameters QualiMap command line qualimap bamqc -bam /content/SRR13697069_ref.sorted.bam -c -nw 400 -hm 3 Alignment Command line: minimap2 -ax map-ont sequence.fasta SRR13697136_1.fastq Draw chromosome limits: yes Analyze overlapping paired-end reads: no Program: minimap2 (2.26-r1175) Analysis date: Tue Nov 21 04:34:01 UTC 2023 Size of a homopolymer: 3 Skip duplicate alignments: no Number of windows: 400 BAM file: /content/SRR13697069_ref.sorted.bam Summary Globals Reference size 189,443 Number of reads 1,602,835 Mapped reads 58,140 / 3.63% Unmapped reads 1,544,695 / 96.37% Mapped paired reads 0 / 0% Secondary alignments 1,359 Supplementary alignments 121 / 0.01% Read min/max/mean length 0 / 126 / 123.99 Duplicated reads (estimated) 51,386 / 3.21% Duplication rate 86.41% Clipped reads 2,093 / 0.13% ACGT Content Number/percentage of A's 1,700,689 / 24.05% Number/percentage of C's 1,516,618 / 21.45% Number/percentage of T's 2,229,273 / 31.53% Number/percentage of G's 1,623,933 / 22.97% Number/percentage of N's 19 / 0% GC Percentage 44.42% Coverage Mean 37.3231 Standard Deviation 212.7464 Mapping Quality Mean Mapping Quality 2.38 Mismatches and indels General error rate 0.14% Insertions 152 Mapped reads with at least one insertion 0.26% Quast (Cobertura) Herramienta de evaluación del ensamblaje del genoma. QUAST evalúa los ensamblajes del genoma calculando varias métricas. Funciona con y sin genomas de referencia. La herramienta acepta múltiples ensamblajes, por lo http://qualimap.bioinfo.cipf.es/ http://qualimap.bioinfo.cipf.es/ https://localhost:8080/ 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 23/40 Deletions 21 Mapped reads with at least one deletion 0.04% Homopolymer indels 11.56% Chromosome stats Name Length Mapped bases Mean coverage Standard deviation AH015299.2 1464 0 0 0 NC_001658.1 6376 0 0 0 KC505249.1 2990 0 0 0 KC505250.1 1978 0 0 0 KC505251.1 2920 0 0 0 KC505252.1 5419 0 0 0 MN385565.1 5386 6791893 1,261.0273 216.7479 OR753134.1 35184 0 0 0 NC_029119.1 127726 278702 2.182 12.4946 Coverage across reference Coverage Histogram Coverage Histogram (0-50X) Genome Fraction Coverage Duplication Rate Histogram Mapped Reads Nucleotide Content Mapped Reads GC-content Distribution Mapped Reads Clipping Profile Homopolymer Indels Mapping Quality Across Reference Mapping Quality Histogram Contents Input data & parameters Summary que es adecuada para la comparación. !quast.py --threads 2 -r /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta /content/drive/MyDrive/vi !quast.py --threads 2 /content/SRR13697069_ref.sorted.bam -r sequence.fasta /usr/local/bin/quast.py --threads 2 /content/SRR13697069_ref.sorted.bam -r sequence.fasta Version: 5.2.0 System information: OS: Linux-5.15.120+-x86_64-with-glibc2.35 (linux_64) Python version: 3.10.10 CPUs number: 2 Started: 2023-11-21 04:42:12 Logging to /content/quast_results/results_2023_11_21_04_42_12/quast.log CWD: /content Main parameters: MODE: default, threads: 2, min contig length: 500, min alignment length: 65, min alignment IDY: 95.0, ambiguity: one, min local misassembly length: 200, min extensive misassembly length: 1000 WARNING: Can't draw plots: python-matplotlib is missing or corrupted. Reference: /content/sequence.fasta ==> sequence Contigs: Pre-processing... 'utf-8' codec can't decode byte 0x8b in position 1: invalid start byte Traceback (most recent call last): File "/usr/local/bin/quast.py", line 309, in return_code = main(sys.argv[1:]) File "/usr/local/bin/quast.py", line 107, in main contigs_fpaths, old_contigs_fpaths = qutils.correct_contigs(contigs_fpaths, corrected_dirpath, labe File "/usr/local/lib/python3.10/site-packages/quast_libs/qutils.py", line 181, in correct_contigs old_fpaths, corr_fpaths, broken_scaffold_fpaths, logs, is_fatal_error = run_parallel(parallel_corre File "/usr/local/lib/python3.10/site-packages/quast_libs/qutils.py", line 1090, in run_parallel results_tuples = Parallel(**parallel_args)(delayed(_fn)(*args) for args in fn_args) File "/usr/local/lib/python3.10/site-packages/joblib/parallel.py", line 1863, in __call__ return output if self.return_generator else list(output) File "/usr/local/lib/python3.10/site-packages/joblib/parallel.py", line 1792, in _get_sequential_outp res = func(*args, **kwargs) File "/usr/local/lib/python3.10/site-packages/quast_libs/qutils.py", line 240, in parallel_correct_co lengths = get_lengths_from_fasta(contigs_fpath, label) File "/usr/local/lib/python3.10/site-packages/quast_libs/qutils.py", line 140, in get_lengths_from_fa lengths = fastaparser.get_chr_lengths_from_fastafile(contigs_fpath).values() File "/usr/local/lib/python3.10/site-packages/quast_libs/fastaparser.py", line 103, in get_chr_length for raw_line in fasta_file: File "/usr/local/lib/python3.10/codecs.py", line 322, in decode (result, consumed) = self._buffer_decode(data, self.errors, final) UnicodeDecodeError: 'utf-8' codec can't decode byte 0x8b in position 1: invalid start byte ERROR! exception caught! In case you have troubles running QUAST, you can write to quast.support@cab.spbu.ru or report an issue on our GitHub repository https://github.com/ablab/quast/issues Please provide us with quast.log file from the output directory. mailto:quast.support@cab.spbu.ru https://github.com/ablab/quast/issues 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 24/40 Summary Coverage across reference Coverage Histogram Coverage Histogram (0-50X) Genome Fraction Coverage Duplication Rate Histogram Mapped Reads Nucleotide Content Mapped Reads GC-content Distribution Mapped Reads Clipping Profile Homopolymer Indels Mapping Quality Across Reference Mapping Quality Histogram 2023/11/21 04:34:02 Generated by QualiMap v.2.2.2-dev import IPython IPython.display.HTML(filename='quast_results/latest/report.html') 2. Búsqueda con Blast La herramienta básica de búsqueda de alineación local (BLAST) encuentra regiones de similitud local entre secuencias. El programa compara secuencias de nucleótidos o proteínas con bases de datos de secuencias y calcula la importancia estadística de las coincidencias. BLAST se puede utilizar para identi�car secuencias nuevas o desconocidas, inferir relaciones funcionales y evolutivas entre secuencias e identi�car miembros de familias de genes. El resultado de una búsqueda BLAST será un conjunto alineado de secuencias potencialmente relacionadas clasi�cadas según su similitud. Figura. Búsqueda con blast. Una coincidencia su�cientemente cercana entre subsecuencias (indicada por �echas en la �gura anterior, aunque las coincidencias suelen ser más largas de lo que se ilustra aquí) se denomina par de puntuación alta (HSP), mientras que se dice que una secuencia de consulta alcanza una secuencia objetivo si comparten una. o más PAS. A veces, sin embargo, el término "golpe" se utiliza de manera vaga, sin diferenciar entre los dos. Cada HSP está asociado con una "puntuación de bits" que se basa en la similitud de las subsecuencias determinada por un conjunto particular de reglas. Debido a que en conjuntos de temas más grandes es probable que se encuentren algunas buenas coincidencias por casualidad, cada PAS también está asociado con un “valor E”, que representa el número esperado de coincidencias que uno podría encontrar por casualidad en un conjunto de temas de ese tamaño con esa puntuación o mejor. Por ejemplo, un valor E de 0,05 signi�ca que podemos esperar una coincidencia por casualidad en 1 de cada 20 búsquedas similares, mientras que un valor E de 2,0 signi�ca que podemos esperar 2 coincidencias por casualidad para cada búsqueda similar. Usted puede realizar la búsqueda mediante blast directamente en un entorno grá�co en la página web del NCBI en el siguiente link: https://blast.ncbi.nlm.nih.gov/Blast.cgi https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fblast.ncbi.nlm.nih.gov%2FBlast.cgi 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 25/40 La versión de línea de comando (UNIX o Windows) de BLAST se denomina BLAST+, el cual dispone de muchas aplicaciones y opciones, algunas de las cuales usaremos en este guía. El uso de comandos nos facilita realizar esta labor de una manera sistemática, más fácil y de forma más reproducible. 2.1. Cree la base de datos La mayoría de las veces, cuando usamos BLAST+, deseamos crear una base de datos personalizada o una versión local de una base de datos prefabricada para usar BLAST. NCBI proporciona varias bases de datos prediseñadas que podrá descargar y usar en el siguiente enlace ftp://ftp.ncbi.nlm.nih.gov/blast/db/ . Suelen ser muy grandes y están divididos en varios archivos (observe la cantidad de archivos necesarios para nr). Otra opcines es la base de datos pdb. Este es un conjunto de datos que contiene todas las secuencias de proteínas asociadas con archivos PDB (archivos de estructura de proteínas). A menudo no queremos utilizar todo NR, ni ninguna de las bases de datos prefabricadas proporcionadas por NCBI. Es posible que queramos usar una base de datos personalizada de solo las especies o genes que nos interesan (Este puede ser nuestro caso). Si tenemos un archivo en formato fasta (no alineado) de estas secuencias, podemos crear una base de datos a partir de esto con el comando makeblastdb . Podemos hacer uso de los archivos de refencias que usamos en el apartado de minimap2 , el cual se llama "sequence.fasta". Ahora cree una base de datos de proteínas escribiendo: !mkdir /content/drive/MyDrive/virus_tax/blast/ ftp://ftp.ncbi.nlm.nih.gov/blast/db/ 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 26/40 mkdir: cannot create directory ‘/content/drive/MyDrive/virus_tax/blast/’: No such file or directory !makeblastdb -in /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta -out /content/drive/MyDrive/virus_t Building a new DB, current time: 11/27/2023 14:17:27 New DB name: /content/drive/MyDrive/virus_tax/blast/targetvirusDB New DB title: /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta Sequence type: Nucleotide Keep MBits: T Maximum file size: 3000000000B BLAST options error: File /content/drive/MyDrive/virus_tax/minimap2/sequence.fasta does not exist Usted puede también descargar la base de datos para genomas de referencia (incluida virus) directamente de la página del NCBI en el siguiente enlace: https://ftp.ncbi.nlm.nih.gov/refseq/release/ !cd /content/drive/MyDrive/virus_tax/blast/ !wget https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz --2023-11-21 23:46:03-- https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.10, 130.14.250.11, 2607:f220:41e:25 Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.10|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 168645964 (161M) [application/x-gzip] Saving to: ‘viral.1.1.genomic.fna.gz’ viral.1.1.genomic.f 100%[===================>] 160.83M 19.7MB/s in 6.7s 2023-11-21 23:46:10 (24.1 MB/s) - ‘viral.1.1.genomic.fna.gz’ saved [168645964/168645964] !gunzip /content/viral.1.1.genomic.fna.gz !makeblastdb -in /content/viral.1.1.genomic.fna -out targetvirusDB2 -dbtype nucl Building a new DB, current time: 11/21/2023 23:47:07 New DB name: /content/targetvirusDB2 New DB title: /content/viral.1.1.genomic.fna Sequence type: Nucleotide Keep MBits: T Maximum file size: 1000000000B Adding sequences from FASTA; added 18719 sequences in 4.40241 seconds. makeblastdb : Comando para crear una base de datos de voladuras a partir de un archivo FASTA. dbtype nucl: Especi�ca que el tipo de base de datos es de nucleótidos. parse_seqids : Comando de análisis de identi�cación de secuencia para asociar una secuencia con un nodo taxonómico. https://colab.research.google.com/corgiredirector?site=https%3A%2F%2Fftp.ncbi.nlm.nih.gov%2Frefseq%2Frelease%2F https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 27/40 Una vez creada la base de datos econtrará varios archivos, todos los cuales comienzan con 'pdbaa'. Cada uno de ellos permite realizar un tipo diferente de búsqueda. Notarás que cada uno tiene una 'p' después del punto. Esto indica que se ha creado una base de datos de proteínas. Veremos más adelante cómo crear nosotros mismos esta base de datos. 2.2. Consulte sus secuencias en la base de datos Las estrategias de búsqueda BLAST+ se ejecutan escribiendo el tipo que desea en la línea de comando seguido de las opciones de entrada. Esto incluye blastn, blastp, blastx, tblastn y tblastx. Asegúrese de saber qué estrategia de búsqueda es adecuada para sus datos y tipo de base de datos. Antes de realizar la búsqueda deberá tener sus lecturas en formato .fasta . Por consguiente, cambie el formato de sus archivos, puede hacer uso del script reformat.sh de bbmap . %%shell for i in {02..09}; do reformat.sh in=/content/drive/MyDrive/virus_hac/concatenated/cat${i}.fastq out=/content/drive/MyDrive/vir done %%shell for i in {10..16}; do reformat.sh in=/content/drive/MyDrive/virus_hac/concatenated/cat${i}.fastq out=/content/drive/MyDrive/vir done !reformat.sh in=SRR13697069_1.fastq out=SRR13697069_1.fasta 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 28/40 java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697069_1.fastq ou Executing jgi.ReformatReads [in=SRR13697069_1.fastq, out=SRR13697069_1.fasta] Input is being processed as unpaired Input: 839990 reads 83775245 bases Output: 839990 reads (100.00%) 83775245 bases (100.00%) Time: 4.066 seconds. Reads Processed: 839k 206.61k reads/sec Bases Processed: 83775k 20.61m bases/sec Hay muchas maneras de modi�car cómo se ejecuta el blast. Con la opción help imprimirá en la pantalla una gran cantidad de texto, detallando todas las opciones que podemos cambiar durante la búsqueda. %%shell for sample in *.fastq do reformat.sh in=$sample out=${sample%.fastq}.fasta done java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697067_1.fastq Executing jgi.ReformatReads [in=SRR13697067_1.fastq, out=SRR13697067_1.fasta] Input is being processed as unpaired Input: 3455219 reads 412775190 bases Output: 3455219 reads (100.00%) 412775190 bases (100.00%) Time: 24.104 seconds. Reads Processed: 3455k 143.35k reads/sec Bases Processed: 412m 17.12m bases/sec java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697067_2.fastq Executing jgi.ReformatReads [in=SRR13697067_2.fastq, out=SRR13697067_2.fasta] Input is being processed as unpaired Input: 3455219 reads 379845987 bases Output: 3455219 reads (100.00%) 379845987 bases (100.00%) Time: 12.428 seconds. Reads Processed: 3455k 278.02k reads/sec Bases Processed: 379m 30.56m bases/sec java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697069_1.fastq Executing jgi.ReformatReads [in=SRR13697069_1.fastq, out=SRR13697069_1.fasta] Exception in thread "main" java.lang.AssertionError: File SRR13697069_1.fasta exists and overwrite=f at shared.Tools.testOutputFiles(Tools.java:128) at jgi.ReformatReads.(ReformatReads.java:394) at jgi.ReformatReads.main(ReformatReads.java:52) java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697069_2.fastq Executing jgi.ReformatReads [in=SRR13697069_2.fastq, out=SRR13697069_2.fasta] Input is being processed as unpaired Input: 839990 reads 82736548 bases Output: 839990 reads (100.00%) 82736548 bases (100.00%) Time: 4.760 seconds. Reads Processed: 839k 176.47k reads/sec Bases Processed: 82736k 17.38m bases/sec java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697136_1.fastq Executing jgi.ReformatReads [in=SRR13697136_1.fastq, out=SRR13697136_1.fasta] 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N6hxWf&printMode=true 29/40 Input is being processed as unpaired Input: 1602714 reads 198721889 bases Output: 1602714 reads (100.00%) 198721889 bases (100.00%) Time: 5.477 seconds. Reads Processed: 1602k 292.63k reads/sec Bases Processed: 198m 36.28m bases/sec java -ea -Xms300m -cp /usr/local/opt/bbmap-39.01-1/current/ jgi.ReformatReads in=SRR13697136_2.fastq Executing jgi.ReformatReads [in=SRR13697136_2.fastq, out=SRR13697136_2.fasta] Input is being processed as unpaired Input: 1602714 reads 194103252 bases Output: 1602714 reads (100.00%) 194103252 bases (100.00%) Time: 5.169 seconds. Reads Processed: 1602k 310.04k reads/sec !blastn -help USAGE blastn [-h] [-help] [-import_search_strategy filename] [-export_search_strategy filename] [-task task_name] [-db database_name] [-dbsize num_letters] [-gilist filename] [-seqidlist filename] [-negative_gilist filename] [-negative_seqidlist filename] [-taxids taxids] [-negative_taxids taxids] [-taxidlist filename] [-negative_taxidlist filename] [-entrez_query entrez_query] [-db_soft_mask filtering_algorithm] [-db_hard_mask filtering_algorithm] [-subject subject_input_file] [-subject_loc range] [-query input_file] [-out output_file] [-evalue evalue] [-word_size int_value] [-gapopen open_penalty] [-gapextend extend_penalty] [-perc_identity float_value] [-qcov_hsp_perc float_value] [-max_hsps int_value] [-xdrop_ungap float_value] [-xdrop_gap float_value] [-xdrop_gap_final float_value] [-searchsp int_value] [-sum_stats bool_value] [-penalty penalty] [-reward reward] [-no_greedy] [-min_raw_gapped_score int_value] [-template_type type] [-template_length int_value] [-dust DUST_options] [-filtering_db filtering_database] [-window_masker_taxid window_masker_taxid] [-window_masker_db window_masker_db] [-soft_masking soft_masking] [-ungapped] [-culling_limit int_value] [-best_hit_overhang float_value] [-best_hit_score_edge float_value] [-subject_besthit] [-window_size int_value] [-off_diagonal_range int_value] [-use_index boolean] [-index_name string] [-lcase_masking] [-query_loc range] [-strand strand] [-parse_deflines] [-outfmt format] [-show_gis] [-num_descriptions int_value] [-num_alignments int_value] [-line_length line_length] [-html] [-sorthits sort_hits] [-sorthsps sort_hsps] [-max_target_seqs num_sequences] [-num_threads int_value] [-mt_mode int_value] [-remote] [-version] DESCRIPTION Nucleotide-Nucleotide BLAST 2.12.0+ Use '-help' to print detailed descriptions of command line arguments A menudo, si trabajamos con muchas secuencias, queremos que sea más fácil obtener los mejores resultados en un formato fácil de leer. Podemos hacer esto limitando la cantidad de resultados devueltos. A menudo, esto se realiza cambiando el número de alineaciones mostradas y/o el límite del valor e. %cd /content/drive/MyDrive/virus_tax/blast/ 1/19/24, 9:49 AM Curso_virus_parte2_taxonomía.ipynb - Colaboratory https://colab.research.google.com/drive/1xTCn82RkoV4Qn6OOX3XcKPLEbidf_3Dx?authuser=2#scrollTo=BMWw68N