# Nanopore bacterial genome assembly pipeline ## Transfer genome data to server with gdown `Gdown` downloads a public file/folder from Google Drive. `Gdown` provides what curl/wget doesn't for Google Drive: * Skip the security notice allowing you to download large files (`curl`/`wget` fails) * Recursive download of files in a folder (maximum 50 files per folder) * Specify download file format for Google Slides/Sheet/Docs like PDF/XML/CSV. ### Instala gdown si no lo tienes `pip install gdown` Si aparece error: ``` WARNING: The script gdown is installed in '/home/aguerra/.local/bin' which is not on PATH. ``` **¿Qué significa?** El programa gdown se instaló en la carpeta `/home/aguerra/.local/bin`, pero esa carpeta no está en la lista de lugares donde tu sistema busca comandos (la variable PATH). La consecuencia: Si escribes `gdown` ahora mismo en tu terminal, te dirá "command not found" (comando no encontrado). Para solucionar esto, necesitas añadir esa carpeta a tu **PATH**. Para hacerlo de forma permanente para que no tengas que repetirlo cada vez que abras una nueva terminal, abre tu archivo de configuración de la terminal. El más común es `.bashrc`. Ejecuta este comando para abrirlo con el editor de texto `nano`: `nano ~/.bashrc` Añade la línea al final del archivo. Usa las flechas para bajar hasta el final del archivo y pega la siguiente línea: `export PATH="$HOME/.local/bin:$PATH"` Guarda y cierra el archivo. Presiona `Ctrl + X` para salir. Finalmente, te preguntará si quieres guardar, presiona `Y` (de Yes). Presiona Enter para confirmar el nombre del archivo. "Refresca" tu terminal. Para que los cambios se apliquen inmediatamente, ejecuta uno de estos dos comandos (el primero es el más común): ``` source ~/.bashrc ``` **Recomendacion:** Cierra y abre la terminal, iniciando nuevamente sesion en el serve -si se esta corriendo alli- para que permita observar los cambios Para verificar que todo está correcto, escribe: ``` which gdown ``` Si todo salió bien, te debería devolver la ruta: `/home/aguerra/.local/bin/gdown`. Ahora ya puedes usar el comando `gdown` para descargar tus archivos desde Google Drive. ### Descarga el archivo usando su ID de Google Drive ``` gdown LINK_DRIVE_FOLDER -O OUTPUTDIR_NAME --folder ``` Se generara un folder llamado `OUTPUTDIR_NAME` con los archivos que estan dentro de ese folder del drive. LISTO PARA USAR. # Genome assembly using flye Despues de verificar que los reads estan correctos y se colocan en un unico archivo, utilizaremos el programa `flye` para emsamblar el genoma. ## Installation ``` conda create -n flye conda activate flye conda install flye ``` ## Analysis ``` flye --nano-hq all_reads.fastq --out-dir flye_assembly --genome-size 5.0m ``` # Genome assembly refirement using medaka Despues del emsamblaje, se debe hacer un tipo de control de calidad del emsamblaje mediante el programa `medaka`, el cual generar un assembly `consensus.fasta` correspondiente al fasta emsamblado pero con las respectivas correciones de errores. ## Installation ``` conda create -n medaka -c conda-forge -c nanoporetech -c bioconda medaka conda activate medaka ``` ## Analysis For native data with bacterial modifications, such as bacterial isolates, metagenomic samples, or plasmids expressed in bacteria, there is a research model that shows improved consensus accuracy. This model is compatible with several basecaller versions for the R10 chemistries. By adding the flag `--bacteria` the bacterial model will be selected if it is compatible with the input basecallers: `medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${OUTDIR} -t ${NPROC} --bacteria` The basecall input corresponde al archivo `all_read.fastq` que generamos cuando juntamos todos los fastq que se generaron del secuenciamiento. Se generaran multiples archivos en la carpeta `outdir` pero el que nos interesa es el consensus.fasta # Quality Assessment of the assembly genome using QUAST Una vez realizada el refinement del assembly del genoma, es necesario analizar la calidad mediante algunos parametros como: * Number of large contigs (i.e., longer than 500 bp) and total length of them. * Length of the largest contig. * N50 (length of a contig, such that all the contigs of at least the same length together cover at least 50% of the assembly). * Number of predicted genes, discovered either by GeneMark.hmm (for prokaryotes), GeneMark-ES or GlimmerHMM (for eukaryotes), or MetaGeneMark (for metagenomes). ## Installation ``` conda create -n quast conda activate quast conda install -c bioconda quast ``` ## Analysis ``` quast.py consensus.fasta -o quast_report ``` Se generara la carpeta con los quality reports # Contigs extraction in case genome assembly genero multiple contigs, incluyendo posibles plasmidos, etc. Utilizar el siguiente comando: `awk '/^><nombre_del_contig>$/{p=1} /^>/ && !/^><nombre_del_contig>$/{p=0} p' consensus.fasta > contig_extraido.fasta` ¿Cómo funciona? Le dice a `awk` que empiece a imprimir cuando encuentre el encabezado que buscas y que se detenga cuando encuentre el siguiente encabezado de otro contig. Por ejemplo, para extraer `contig_1`: ``` awk '/^>contig_1$/{p=1} /^>/ && !/^>contig_1$/{p=0} p' consensus.fasta > contig_1.fasta ``` Si la idea es extraer todos los contigs dentro del multifasta: ``` awk '/^>/ {out = substr($1, 2) ".fasta"} {print >> out}' consensus.fasta ``` Esto generara en la carpeta donde esta el consensus los multiples contigs con el header del fasta. # Evaluate contigs size ALGO ASI ``` for f in *.fasta; do awk '/^>/ {next} {l += length($0)} END {print FILENAME, l}' $f; done ``` Generara algo asi: ``` consensus_chromo.fasta 4459216 contig_1.fasta 4194028 contig_2.fasta 146552 contig_3.fasta 107201 contig_5.fasta 4878 contig_7.fasta 4673 contig_9.fasta 1884 ``` # Reorients complete circular genome using DNApler https://github.com/gbouras13/dnaapler # creates empty conda environment conda create -n dnaapler_env # activates conda environment conda activate dnaapler_env # installs dnaapler conda install -c bioconda dnaapler # runs dnaapler all dnaapler all -i input_mixed_contigs.fasta -o output_directory_path -p my_bacteria_name -t 8 # runs dnaapler all with a gfa file from e.g. Flye, Unicycler or Autocycler dnaapler all -i assembly.gfa -o output_directory_path -p my_bacteria_name -t 8 Se generarn multiples output incluyendo los archivos: * dnaapler_1767404886.451651.log * logs * UFO-VIO_all_reorientation_summary.tsv * UFO-VIO_MMseqs2_output.txt * UFO-VIO_reoriented.fasta El reoriented.fasta seria el fasta que debemos usar y que esta ya reorientado. # Genome annotation using bakta ## Installation ``` conda create -n bakta conda activate bakta conda install -c conda-forge -c bioconda bakta ``` ### Database download Bakta requires a mandatory database which is publicly hosted at Zenodo: DOI We provide 2 types: full and light. To get best annotation results and to use all features, we recommend using the full (default). If you seek for maximum runtime performance or if download time/storage requirements are an issue, please try the light version. Further information is provided in the database section below. List available DB versions (available as either full or light): `bakta_db list` To download the most recent compatible database version we recommend to use the internal database download & setup tool: `bakta_db download --output <output-path> --type [light|full]` Of course, the database can also be downloaded and installed manually: `wget https://zenodo.org/record/14916843/files/db-light.tar.xz` `bakta_db install -i db-light.tar.xz` If required, or desired, the AMRFinderPlus DB can also be updated manually: `amrfinder_update --force_update --database db-light/amrfinderplus-db/` ## Analysis ### For a single genome: ``` bakta --db path/to/dbdir --prefix PREFIXname --output outputdir/name --verbose "$file" ``` ### For multiple genomes: ``` #!/bin/bash outputdir="bakta" fnadir="/path/to/fna/files/fna_coded" extension=".fna" cpus=4 dbdir="/path/to/bakta/db" # Create output directory if it doesn't exist if [ ! -d "$outputdir" ]; then mkdir "$outputdir" fi allfna=("${fnadir}"/*"${extension}") for file in "${allfna[@]}" do bname=$(basename "$file" "${extension}") bakta --db "$dbdir" --prefix "$bname" --output "$outputdir/$bname" --verbose "$file" done ``` # Extraction of 16S rRNA ribosomal gene https://github.com/tseemann/barrnap ## Installation ``` conda create -n barrnap conda activate barrnap conda install -c bioconda -c conda-forge barrnap ``` ## Analysis Extraction of ribosomal RNA sequences: ``` barrnap --kingdom bac --threads 8 genome.fasta --outseq rrna_sequences.fasta > rrna_locations.gff ``` Dependiendo del modelo se extraeran los diferentes genes, incluyendo para bacteria por ejemplo 16s y 18s usando `bac` El archivo rrna_sequences.fasta: ![image](https://hackmd.io/_uploads/rJZ3rjbRlx.png) Como se observa multiples 16s y 18s fueron identificados ya que generalmente existen diferentes copias en el genoma de estos genes. El archivo rrna_locations.gff: ![image](https://hackmd.io/_uploads/r1qeLs-0xx.png) Este archivo contiene las informaciones de la localizacion de estos genes en el genoma. En este caso solo usamos el contig 1. Para renombrar las secuencias utilizando solo el nombre de la especie, segun lo que ya ncbi da que es que el header del fasta viene en el orden: ACCESSION GENERO EPITETO OTRAS INFORMACIONES ADICIONALES. El siguiente comando tomara unicamente el GENERO y EPITETO para renombrar. Siempre verificar manualmente que todo este bien. Solo colocar el siguiente comando en el CL cambiando los nombres del input y output: ``` sed -E '/^>/ s/^>[^ ]+ ([^ ]+) ([^ ]+).*/ >\1_\2 /' 16s_sequences_raw.fasta > 16s_sequences_final.fasta ``` NOTA: esto genera el header pero con un espacio al comiendo, **lo que podria dar problemas**, es necesario modificar esto o con otro comando que permita quitar el espacio, o que genere los nombres de los header sin esos espacios. # Sequence aligment using MAFFT ## Installation ``` sudo apt install mafft ``` - **--op # :** Gap opening penalty, default: 1.53 - **--ep # :** Offset (works like gap extension penalty), default: 0.0 - **--maxiterate # :** Maximum number of iterative refinement, default: 0 - **--clustalout :** Output: clustal format, default: fasta - **--reorder :** Outorder: aligned, default: input order - **--quiet :** Do not report progress - **--thread # :** Number of threads (if unsure, --thread -1) - **--dash :** Add structural information (Rozewicki et al, submitted) **If unsure which option to use:** ``` mafft --auto in > out ``` **To run standard aligment:** ``` mafft ---reorder ---adjustdirection ---auto input.fasta > output.fas ``` - **--reorder**: Outorder: aligned, default: input order - **--adjustdirection:** organize the direction of the sequences based on the first sequence - **--auto:** selection of auto model for the aligment - **in > out:** input.fas and output.fas # Generate the phylogenetic tree using IQ-Tree ## Installation ### Create a enviroment for IQ-Tree using conda: ``` conda create -n iqtree ``` ### Activate the enviroment: ``` conda activate iqtree ``` ### Install the packages using bioconda: ``` conda install -c bioconda iqtree ``` ## General ML tree analysis Combine ModelFinder, tree search, SH-aLRT test and ultrafast bootstrap with 1000 replicates: ``` iqtree -s example.phy -B 1000 -alrt 1000 ``` # Roary https://github.com/sanger-pathogens/Roary ## Installation ## Analysis To run the analysis using gff files from PROKKA annotation: `roary -f ./roary_output -e -r -n -v -p 7 *.gff` ## Plots ### Script para generar curva de rarefaccion ``` library(ggplot2) # Archivos de datos archivo_conserved = "number_of_conserved_genes.Rtab" archivo_total = "number_of_genes_in_pan_genome.Rtab" # Leer datos y calcular medias conserved = colMeans(read.table(archivo_conserved)) total = colMeans(read.table(archivo_total)) # Crear un dataframe combinado datos = data.frame( Genes = c(conserved, total), Genomas = c(1:length(conserved), 1:length(total)), Tipo = c(rep("Genes Conservados", length(conserved)), rep("Genes Totales", length(total))) ) # Crear la gráfica p <- ggplot(data = datos, aes(x = Genomas, y = Genes, color = Tipo)) + geom_line(size = 1.2) + theme_classic() + xlab("Number of genomes") + ylab("Number of genomes") + theme_bw(base_size = 16) + theme(legend.position = "right") # Guardar la gráfica en formato PNG ggsave(filename = "pangenoma_unico.png", plot = p, width = 10, height = 10) # Guardar la gráfica en formato SVG ggsave(filename = "pangenoma_unico.svg", plot = p, width = 10, height = 10) ``` ### Script para generar grafico UpSet enfocado en un genoma, en este caso UFO-VIO ``` #install.packages(c("UpSetR", "tidyverse")) #install.packages("svglite") # --- 1. Cargar las librerías --- library(tidyverse) library(UpSetR) library(svglite) # Usamos svglite que es más robusto # (Asegúrate de tenerlo instalado: install.packages("svglite")) # --- 2. Cargar y Preparar los Datos --- # (Tu código de preparación de datos no cambia) genome_cols <- c( "A56AF", "BASUSDA_45", "CV1", "CV8", "Cv1", "DSM26508", "SAM215", "chromo_utfo" ) genome_of_interest <- "UFO_VIO" tab20_colors <- c( "#1F77B4", "#FF7F0E", "#2CA02C", "#D62728", "#8C564B", "#E377C2", "#7F7F7F", "#9467BD" ) tryCatch({ df <- read_csv("gene_presence_absence.csv") df_binary <- df %>% select(all_of(genome_cols)) %>% mutate(across(everything(), ~ as.integer(!is.na(.) & . != ""))) %>% rename(!!genome_of_interest := chromo_utfo) df_filtered <- df_binary %>% filter(!!sym(genome_of_interest) == 1) if(nrow(df_filtered) == 0) { stop(paste("No se encontró ningún gen para '", genome_of_interest, "'."), call. = FALSE) } # --- 3. Generar el Gráfico y guardarlo en una Variable --- print("Generando objeto de gráfico UpSet...") # ¡CAMBIO CLAVE 1: Asignar el gráfico a una variable 'p'! # Esto NO lo dibujará todavía, solo crea el objeto. p <- upset( as.data.frame(df_filtered), nsets = ncol(df_filtered), order.by = "freq", decreasing = TRUE, main.bar.color = tab20_colors[8], sets.bar.color = tab20_colors, matrix.color = tab20_colors[8], mainbar.y.label = "Genes en Intersección (con UFO-VIO)", sets.x.label = "Genes Compartidos con UFO-VIO", text.scale = 1.2 ) # --- 4. Mostrar el Gráfico en la Ventana "Plots" --- # ¡CAMBIO CLAVE 2: Imprimir 'p' en la pantalla! # Esto fuerza a R a dibujarlo en la ventana de Plots de RStudio. print("Mostrando gráfico en la ventana de 'Plots'...") print(p) # --- 5. Guardar el Gráfico en SVG --- print("Guardando gráfico en 'UFO_VIO_UpSet.svg'...") # Abrir el dispositivo SVG svglite::svglite("UFO_VIO_UpSet.svg", width = 12, height = 7) # ¡CAMBIO CLAVE 3: Imprimir 'p' de nuevo! # Esta vez, se dibuja en el archivo SVG que está "abierto". # Esta es la corrección que evita el SVG en blanco. print(p) # Cerrar el dispositivo (¡muy importante!) dev.off() print("¡Gráfico guardado exitosamente!") }, error = function(e) { print(paste("Error:", e$message)) if (grepl("Could not find function", e$message)) { print("Asegúrate de que los paquetes 'tidyverse', 'UpSetR' y 'svglite' están instalados y cargados.") } else { print("Asegúrate de que 'gene_presence_absence.csv' esté en tu directorio.") } })# --- INICIO DE LA MODIFICACIÓN --- # 1. Abrir el dispositivo SVG usando svglite (más robusto) svglite::svglite("UFO_VIO_UpSet.svg", width = 12, height = 7) # 2. Ejecutar el comando del gráfico (sin cambios) upset( as.data.frame(df_filtered), nsets = ncol(df_filtered), order.by = "freq", decreasing = TRUE, main.bar.color = tab20_colors[8], sets.bar.color = tab20_colors, matrix.color = tab20_colors[8], mainbar.y.label = "Genes en Intersección (con UFO-VIO)", sets.x.label = "Genes Compartidos con UFO-VIO", text.scale = 1.2 ) # 3. Cerrar el dispositivo SVG (sin cambios) dev.off() # --- FIN DE LA MODIFICACIÓN --- print("¡Gráfico guardado exitosamente con svglite!") ``` ### Script para generar graficos: modificacion del py Crear el archivo y correr el script usando: `python roary_plots_custom.py [tu_arbol.newick] [tu_matriz.csv]` El arbol newick puede ser generado usando IQTREE y el archivo `core_gene_alignment.aln` A continuacion el script ``` #!/usr/bin/env python # Copyright (C) <2015> EMBL-European Bioinformatics Institute # Modified for enhanced aesthetics and SVG output by Gemini. # Este script es una versión modificada de roary_plots.py para generar # gráficos de mayor calidad visual, listos para una publicación. # Para usarlo, ejecuta en la terminal: # python roary_plots_custom.py [tu_arbol.newick] [tu_matriz.csv] __author__ = "Marco Galardini & Gemini" __version__ = '0.9.0' def get_options(): """Configura y parsea los argumentos de la línea de comandos.""" import argparse description = "Create enhanced plots from Roary outputs." parser = argparse.ArgumentParser(description=description, prog='roary_plots_custom.py') parser.add_argument('tree', action='store', help='Newick Tree file (e.g., accessory_binary_genes.fa.newick)') parser.add_argument('spreadsheet', action='store', help='Roary gene presence/absence spreadsheet (e.g., gene_presence_absence.csv)') parser.add_argument('--labels', action='store_true', default=False, help='Add node labels to the tree (up to 10 chars)') # El formato por defecto ahora es SVG parser.add_argument('--format', choices=('png', 'tiff', 'pdf', 'svg'), default='svg', help='Output format [Default: svg]') parser.add_argument('-N', '--skipped-columns', action='store', type=int, default=14, help='First N columns of Roary\'s output to exclude [Default: 14]') parser.add_argument('--version', action='version', version='%(prog)s ' + __version__) return parser.parse_args() if __name__ == "__main__": options = get_options() import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import seaborn as sns import pandas as pd import numpy as np from Bio import Phylo # Importamos una función para crear mapas de color personalizados from matplotlib.colors import ListedColormap # --- Configuración de Estilo Global --- # --- MODIFICACIÓN: Cambiamos el estilo a 'ticks' para un look más limpio --- # Esto elimina la cuadrícula de fondo pero mantiene las marcas de los ejes. sns.set_theme(style="ticks") # --- MODIFICACIÓN: Aseguramos que todo el texto sea negro por defecto --- plt.rc('text', color='black') plt.rc('axes', labelcolor='black', edgecolor='black') plt.rc('xtick', color='black') plt.rc('ytick', color='black') plt.rc('figure', dpi=300) # --- Carga y Preparación de Datos --- print("-> Loading tree and Roary spreadsheet...") t = Phylo.read(options.tree, 'newick') roary = pd.read_csv(options.spreadsheet, low_memory=False) roary.set_index('Gene', inplace=True) roary.drop(list(roary.columns[:options.skipped_columns-1]), axis=1, inplace=True) roary.replace('.{2,100}', 1, regex=True, inplace=True) roary.replace(np.nan, 0, regex=True, inplace=True) print("-> Data loaded successfully.") # --- Cálculo de Categorías del Pangenoma (se hace una vez para todos los gráficos) --- num_genomes = roary.shape[1] core_threshold = num_genomes * 0.99 soft_core_threshold = num_genomes * 0.95 shell_threshold = num_genomes * 0.15 gene_sums = roary.sum(axis=1) # --- Gráfico 1: Histograma de Frecuencia del Pangenoma (Mejorado) --- print("-> Generating pangenome frequency plot...") freq_df = pd.DataFrame({'frequency': gene_sums}) def assign_category_hist(freq): if freq >= core_threshold: return 'Core genes' elif freq >= soft_core_threshold: return 'Soft-core genes' elif freq >= shell_threshold: return 'Shell genes' else: return 'Cloud genes' freq_df['Category'] = freq_df['frequency'].apply(assign_category_hist) comp_palette = sns.color_palette("Pastel1", 4) color_map = { 'Cloud genes': comp_palette[0], 'Shell genes': comp_palette[1], 'Soft-core genes': comp_palette[2], 'Core genes': comp_palette[3] } plt.figure(figsize=(10, 6)) sns.histplot(data=freq_df, x='frequency', hue='Category', multiple="stack", bins=num_genomes, palette=color_map, edgecolor="white", linewidth=0.5, hue_order=['Cloud genes', 'Shell genes', 'Soft-core genes', 'Core genes']) plt.xlim(-0.5, num_genomes + 0.5) # Ajuste para centrar las barras correctamente plt.title('Pangenome Gene Frequency', fontsize=16, fontweight='bold') plt.xlabel('Number of Genomes', fontsize=12, fontweight='bold') plt.ylabel('Number of Genes', fontsize=12, fontweight='bold') # --- MODIFICACIÓN: Eliminamos las líneas superiores y derechas del recuadro --- sns.despine() plt.tight_layout() plt.savefig(f'pangenome_frequency.{options.format}') plt.clf() # --- Gráfico 2: Árbol Filogenético y Matriz (Mejorado) --- print("-> Generating tree and presence/absence matrix plot...") roary_sorted = roary.loc[roary.sum(axis=1).sort_values(ascending=False).index] roary_sorted = roary_sorted[[x.name for x in t.get_terminals()]] fig = plt.figure(figsize=(17, 10)) gs = fig.add_gridspec(1, 40, wspace=0) ax_tree = fig.add_subplot(gs[0, 0:10]) with plt.rc_context({'lines.linewidth': 0.5}): Phylo.draw(t, axes=ax_tree, show_confidence=False, label_func=lambda x: None, xticks=([],), yticks=([],), ylabel=('',), xlabel=('',), axis=('off',), title=(f'Phylogenetic Tree\n({roary.shape[1]} strains)',), do_show=False) ax_matrix = fig.add_subplot(gs[0, 10:40]) # --- MODIFICACIÓN: Nueva paleta de colores blanco y verde pastel --- # Usamos blanco para la ausencia (0) y un verde pastel para la presencia (1). # El segundo color de la paleta "Pastel1" (#b3e2cd) es un verde claro. pastel_cmap = ListedColormap(["#e9f7f1ff", "#51be8eff"]) ax_matrix.matshow(roary_sorted.T, cmap=pastel_cmap, aspect='auto', interpolation='none') ax_matrix.set_yticks([]) ax_matrix.set_xticks([]) ax_matrix.axis('off') ax_matrix.set_title(f'Gene Presence/Absence\n({roary.shape[0]} gene clusters)') plt.savefig(f'pangenome_matrix.{options.format}') plt.clf() # --- Cálculo de Categorías para Gráfico de Dona --- core_genes_count = (gene_sums >= core_threshold).sum() soft_core_genes_count = ((gene_sums >= soft_core_threshold) & (gene_sums < core_threshold)).sum() shell_genes_count = ((gene_sums >= shell_threshold) & (gene_sums < soft_core_threshold)).sum() cloud_genes_count = (gene_sums < shell_threshold).sum() pangenome_counts = { 'Core genes': core_genes_count, 'Soft-core genes': soft_core_genes_count, 'Shell genes': shell_genes_count, 'Cloud genes': cloud_genes_count } pangenome_df = pd.DataFrame(list(pangenome_counts.items()), columns=['Category', 'Count']) total = pangenome_df['Count'].sum() # --- Gráfico 3: Gráfico de Dona (Pie Chart Mejorado) --- print("-> Generating pangenome donut chart...") plt.figure(figsize=(10, 8)) colors = sns.color_palette("Pastel1", 4) # --- MODIFICACIÓN: Aseguramos que el texto de los porcentajes sea negro --- pie_wedges, pie_labels, pie_texts = plt.pie(pangenome_df['Count'], labels=pangenome_df['Category'], autopct=lambda p: '{:.0f}'.format(p * total / 100), colors=colors, wedgeprops={'edgecolor': 'white', 'linewidth': 1.5}, pctdistance=0.85, startangle=90, textprops={'color':"black"}) # <-- Texto en negro my_circle = plt.Circle((0, 0), 0.7, color='white') p = plt.gcf() p.gca().add_artist(my_circle) plt.title('Pangenome Composition', fontsize=16, fontweight='bold') plt.savefig(f'pangenome_pie.{options.format}') plt.clf() print("\nAll plots generated successfully!") ``` ## Extract singletons of a interest genome ``` ``` # EggNOGMapper ## eggnogCOGextractor https://github.com/vsmicrogenomics/eggnogCOGextractor/tree/main ## Instalation Solo es descargar el archivo `eggnogCOGextractor.py` para correr el programa ## Analysis Crear una carpeta llama `input` y otra `output`. En la input debe estar el`MM_rqbj1453.emapper.annotations.tsv`, sin embargo, es importante que este no tenga la extension tsv si no que termine en **`.emapper.annotations`** Para correr el analisis correr el siguiente comando: `python eggnogCOGextractor.py ` Se generara un archivo tsv en el directorio output/ de la siguiente forma: ![image](https://hackmd.io/_uploads/SkkuYPwHeg.png) **Este archivo es importante para hacer algunos analisis como barplot**, entre otros. ## Vertical barplot de la frecuencia relativa de COGs based on eggnogCOGextractor output tsv (.emapper.cog.tsv) #### Script: ``` import pandas as pd import matplotlib.pyplot as plt import seaborn as sns # --- CONFIGURACIÓN --- INPUT_FILE = '/home/aguerra/chromobacterium/eggnogmapper/output/MM_wd3sxapu.emapper.cog.tsv' # Nombres para los archivos de salida OUTPUT_PLOT_PNG = 'cog_bar_chart_final.png' OUTPUT_PLOT_SVG = 'cog_bar_chart_final.svg' # --- SCRIPT PRINCIPAL --- print(f"Paso 1: Cargando datos desde '{INPUT_FILE}'...") try: df = pd.read_csv(INPUT_FILE, sep='\t') except FileNotFoundError: print(f"\nERROR: No se encontró el archivo '{INPUT_FILE}'.") exit() except pd.errors.EmptyDataError: print(f"\nERROR: El archivo '{INPUT_FILE}' está vacío.") exit() df = df[df['COGCount'] > 0] df = df.sort_values('COGCount', ascending=False) print("Paso 2: Generando la gráfica de barras...") # Crear la figura (fig) y los ejes (ax) fig, ax = plt.subplots(figsize=(12, 8)) # Crear la gráfica de barras horizontales sns.barplot(data=df, x='COGCount', y='COGDescription', hue='COGDescription', palette='tab20', orient='h', legend=False, ax=ax) # Añadir el número exacto al final de cada barra. for index, value in enumerate(df['COGCount']): ax.text(value, index, f' {value}', va='center', ha='left') # Añadir la letra de la categoría (COGID) en el lado derecho. ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax2 = ax.twinx() ax2.spines['right'].set_visible(False) ax2.set_ylim(ax.get_ylim()) ax2.set_yticks(ax.get_yticks()) ax2.set_yticklabels(df['COGID'], fontsize=10, weight='bold') ax.set_title('Número de Genes por Categoría Funcional COG', fontsize=16, weight='bold') ax.set_xlabel('Número de Genes', fontsize=12) ax.set_ylabel('Categoría COG', fontsize=12) plt.tight_layout() # --- MODIFICACIÓN FINAL: Guardar en ambos formatos --- plt.savefig(OUTPUT_PLOT_PNG, dpi=300, bbox_inches='tight') plt.savefig(OUTPUT_PLOT_SVG, bbox_inches='tight') # SVG es un vector, no necesita dpi print(f"¡Éxito! Gráficas guardadas como '{OUTPUT_PLOT_PNG}' y '{OUTPUT_PLOT_SVG}'") plt.show() ``` Se debe de generar un output svg y png como el siguiente: ![cog_bar_chart_final](https://hackmd.io/_uploads/HyhmhFwrll.png) Realizar las modificaciones especificas para las necesidades de los datos. Se puede modificar la paleta, entre otros detalles visuales de la figura. ## Vertical barplot de los genes especificos de algun COGs category (Ej: V and Q) based on eggnogCOGextractor output tsv (.emapper.cog.tsv) #### Script ``` import pandas as pd import matplotlib.pyplot as plt import seaborn as sns # --- CONFIGURATION --- EGGNOG_ANNOTATIONS_FILE = 'MM_wd3sxapu.emapper.annotations' # Color map for categories of interest COLOR_MAP = { 'Q': '#2eabb8ff', # Color for 'Secondary metabolites' 'V': '#d1d197ff' # Color for 'Defense mechanisms' } CATEGORIES_TO_PLOT = list(COLOR_MAP.keys()) # --- PLOTTING FUNCTION --- def generate_detailed_plot(df, cog_letter, plot_color): """ Filters the dataframe by a COG category and generates a bar plot using a specific color. """ print(f"\n--- Processing Category: {cog_letter} with color {plot_color} ---") category_df = df[df['Cleaned_COG'] == cog_letter].copy() if category_df.empty: print(f"No genes found for category '{cog_letter}'.") return plot_data = category_df.loc[category_df.groupby('Description')['score'].idxmax()] plot_data = plot_data.sort_values('score', ascending=False).head(25) print(f"Found {len(plot_data)} unique functions to plot.") plt.figure(figsize=(12, 10)) sns.barplot(data=plot_data, x='score', y='Description', color=plot_color, orient='h') # --- Plot labels in English --- plt.title(f'Detailed Functions for COG Category: "{cog_letter}"', fontsize=16, weight='bold') plt.xlabel('Annotation Score (Bit-score)', fontsize=12) plt.ylabel('Function Description', fontsize=12) plt.tight_layout() output_filename_svg = f'detailed_plot_category_{cog_letter}.svg' plt.savefig(output_filename_svg, bbox_inches='tight') print(f"Plot saved as '{output_filename_svg}'") plt.close() # --- MAIN SCRIPT --- print(f"Step 1: Loading data from '{EGGNOG_ANNOTATIONS_FILE}'...") try: header_row = 0 with open(EGGNOG_ANNOTATIONS_FILE, 'r') as f: for i, line in enumerate(f): if line.startswith('#query'): header_row = i break main_df = pd.read_csv(EGGNOG_ANNOTATIONS_FILE, sep='\t', header=header_row) except FileNotFoundError: print(f"\nERROR: File not found '{EGGNOG_ANNOTATIONS_FILE}'.") exit() print("Step 2: Cleaning COG categories...") main_df.dropna(subset=['COG_category'], inplace=True) main_df['Cleaned_COG'] = main_df['COG_category'].apply(lambda x: str(x)[0] if str(x) != '-' else None) # Generate a plot for each category of interest for category, color in COLOR_MAP.items(): generate_detailed_plot(main_df, category, color) print("\nAnalysis complete.") ``` Se debe de generar un output svg y png como el siguiente: ![final_cohesive_plot](https://hackmd.io/_uploads/rkAFgRnBee.png) # KEGGaNOG: visualize KEGG Pathways using output results from EggNOG-mapper https://github.com/iliapopov17/KEGGaNOG ## Installation ``` # Linux / WSL / macOS conda create -n kegganog pip -y conda activate kegganog pip install kegganog ``` ## Analysis ### Single mode plot for KEGGs completeness: `KEGGaNOG -i MM_rqbj1453.emapper.annotations.tsv -o output_dir` By the key -i the input (eggNOG-mapper output file) must be provided By the key -o the path to output directory (it will be created if it does not exist) must be provided ### Barplot for KEGGs completeness: **PENDIENTE MODIFICAR EL SCRIPT PARA QUE GENERE EL SVG. ASI COMO LA PALETA** Cambiar la ubicacion del archivo **`SAMPLE_pathways.tsv`** generado por el programa. ``` import kegganog as kgn import pandas as pd # Cargar los datos df = pd.read_csv("/Users/abraham/Downloads/eggNOG_mapper_analysis_CLAS/grouped_clas/SAMPLE_pathways.tsv", sep="\t") # Preparar el objeto de la gráfica kgnbar = kgn.barplot(df, figsize = (8, 10), sort_order="descending", yticks_fontsize = 8, cmap = "Blues", ) # ---- CORRECCIÓN FINAL ---- # 1. Generar la figura y guardarla en una variable 'fig' fig = kgnbar.plotfig() # 2. Usar el método savefig() de la figura para guardarla en un archivo fig.savefig("KEGG_barplot.png", # Nombre del archivo de salida dpi=300, # Para alta resolución bbox_inches='tight') # Asegura que las etiquetas no se corten print("¡Éxito! Gráfica guardada como KEGG_barplot.png") ``` # Orthovenn # PPanGGolin Los spots generados vendrian siendo las islas genomicas que pueden estar presentes en multiples genomas. # CLinker La idea es tomar los gbk que genero bakta y utilizarlos como input para CLinker pero se tomarian solos los ranges de cada isla genomica. Para eso se usan los archivos spot_X_identical_rgps.tsv que continen las informaciones de las coordenadas de los spots correspondientes a los RGP. Una vez que clinker genera el archivo HTML, lo abres en tu navegador. Verás los genomas completos con sus coordenadas. Tú simplemente te desplazas y haces zoom a la región que te interesa (ej. te vas a la marca de 1,179,962 en el genoma gxpsy) para ver tu isla. # PGPg_finder https://github.com/tpellegrinetti/PGPg_finder ## Installation To use the tool, you first need to clone the repository and ensure that all dependencies are installed. ``` git clone https://github.com/tpellegrinetti/PGPg_finder/ ``` After this procedure, you need to performe the installation of dependancies with conda (recomended). The conda will create a separated environment called PGPg_finder. ``` bash install.sh ``` ## Analysis for a single genome ``` python PGPg_finder/PGPg_finder.py -w genome_wf -i GENOMEFILE.fna -o pgp_chromo -t 20 ``` Obtener el py del heatmap que esta en el github o copiar y pegar para crear un nuevo python. Para generar el heatmap y el resumen ya por categoria: `python heatmap.py gene_counts.txt heatmap_output/ pathways_plabase.txt summary.txt ` Es necesario obtener los archivos `pathways_plabase.txt ` `summary.txt` mediante wget con el link que esta en el github o segun la ubicacion en la carpeta donde fue instalado