snakemake para hacer bioinformática: una guía para principiantes (parte 2)

lun 23 enero 2023

Por C. Titus Brown

En la ciencia.

Etiquetas: serpientehacedeslizándose

(La siguiente publicación contiene extractos de Deslizándose hacia la bioinformática con snakemakeHackmd Press, 2023.)

En la Sección 1, presentamos a snakemake como un sistema para (eficiente y efectivamente) ejecutar una serie de comandos de shell.

En la Sección 2, exploraremos una serie de características importantes de la creación de serpientes. Junto con la Sección 1, esta sección cubre el conjunto básico de la funcionalidad de snakemake que necesita saber para aprovechar de manera efectiva el uso de snakemake.

Después de esta sección, estará bien posicionado para escribir algunos flujos de trabajo propios y luego podrá regresar y explorar funciones más avanzadas a medida que las necesite.

Capítulo 4: ejecutar reglas en paralelo

Echemos un vistazo a la sketch_genomes regla desde el ultimo
Snakefile entrada:

rule sketch_genomes:
    input:
        "genomes/GCF_000017325.1.fna.gz",
        "genomes/GCF_000020225.1.fna.gz",
        "genomes/GCF_000021665.1.fna.gz",
    output:
        "GCF_000017325.1.fna.gz.sig",
        "GCF_000020225.1.fna.gz.sig",
        "GCF_000021665.1.fna.gz.sig"
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

Este comando funciona bien tal como está, pero es levemente incómodo, porque, dado que la bioinformática es bioinformática, es probable que queramos agregar más genomas a la comparación en algún momento, y en este momento cada genoma adicional tendrá que agregarse tanto a la entrada como a la salida. No es mucho trabajo, pero es innecesario.

Además, si añadimos un lote de los genomas, este paso podría convertirse rápidamente en un cuello de botella. sourmash sketch puede ejecutarse rápidamente en 10 o 20 genomas, ¡pero se ralentizará si le da 100 o 1000! (De hecho, sourmash sketch escala con la cantidad de genomas, por lo que tomará 100 veces más en 100 genomas que en 1). ¿Hay alguna manera de acelerar eso?

Sí, podemos escribir una regla que se pueda ejecutar para cada genoma y luego dejar que Snakemake la ejecute en paralelo para nosotros.

Comencemos dividiendo esta regla en tres separado normas:

rule sketch_genomes_1:
    input:
        "genomes/GCF_000017325.1.fna.gz",
    output:
        "GCF_000017325.1.fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

rule sketch_genomes_2:
    input:
        "genomes/GCF_000020225.1.fna.gz",
    output:
        "GCF_000020225.1.fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

rule sketch_genomes_3:
    input:
        "genomes/GCF_000021665.1.fna.gz",
    output:
        "GCF_000021665.1.fna.gz.sig"
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

# rest of Snakefile here!

Es prolijo, pero funcionará. Ejecute:

snakemake -j 1 --delete-all plot_comparison
snakemake -j 1 plot_comparison

Antes de seguir modificando el archivo, disfrutemos de los frutos de nuestro trabajo: ¡ahora podemos decirle a snakemake que ejecute más de una regla a la vez!

Intenta escribir esto:

snakemake -j 1 --delete-all plot_comparison
snakemake -j 3 plot_comparison

Si observa detenidamente, debería ver que Snakemake está ejecutando los tres
sourmash sketch dna comandos al mismo tiempo.

Esto es genial y es una de las características prácticas más poderosas de Snakemake: una vez que le dices a SnakeMake que quieres que haga (especificando la(s) salida(s) deseada(s)) y dale a snakemake el conjunto de recetas diciéndolo como hacer cada pasosnakemake descubrirá la forma más rápida de ejecutar todos los pasos necesarios con los recursos que le ha proporcionado.

En este caso, le dijimos a snakemake que podría ejecutar hasta tres trabajos a la vez, con -j 3. También podríamos haberle dicho que ejecutara más trabajos a la vez, pero por el momento solo hay tres reglas que se pueden ejecutar al mismo tiempo: sketch_genomes_1, sketch_genomes_2y
sketch_genomes_3. Esto se debe a que la regla compare_genomes necesita la salida de estas tres reglas para ejecutarse, y del mismo modo plot_genomes necesita la salida de compare_genomes correr. ¡Así que no se pueden ejecutar al mismo tiempo que otras reglas!

Capítulo 5: visualización de flujos de trabajo

¡Visualicemos lo que estamos haciendo! Aquí está la salida de snakemake
--dag plot_comparison
visualizado con el paquete graphviz:

gráfico interm2 de trabajos

Este diagrama muestra la relación entre las reglas que hemos puesto en el Snakefile: compare_genomes toma la salida del sketch_genome
reglas como su propia entrada, y luego plot_comparison utiliza la salida de
compare_genomes para construir su propia parcela.

Un aspecto clave de este gráfico es que le muestra dónde se pueden ejecutar varias reglas al mismo tiempo porque no requieren ni se requieren para las demás: aquí, las tres
sketch_genome normas. ¡Eso es lo que nos permitió ejecutar los tres simultáneamente en el capítulo anterior!

Nota: a veces debe tener una sola regla que se ocupe de todos los genomas; por ejemplo, compare_genomes tiene que comparar todas los genomas, y no hay una forma sencilla de evitarlo. Pero con sketch_genomes¡Tenemos la opción de romper la regla!

Capítulo 6: uso de comodines para hacer que las reglas sean más genéricas

Echemos otro vistazo a uno de esos sketch_genomes_ normas:

rule sketch_genomes_1:
    input:
        "genomes/GCF_000017325.1.fna.gz",
    output:
        "GCF_000017325.1.fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

Hay algo de redundancia allí: la adhesión GCF_000017325.1 aparece dos veces. ¿Podemos hacer algo al respecto?

¡Si podemos! Podemos usar una función de creación de serpientes llamada «comodines», que nos permitirá dar a la creación de serpientes un espacio en blanco para que se complete automáticamente.

Con comodines, le indica a SnakeMake que una parte particular de un nombre de archivo de entrada o salida es un juego justo para sustituciones usando
que rodea el nombre comodín. Vamos a crear un comodín llamado accession
y colóquelo en los bloques de entrada y salida de la regla:

rule sketch_genomes_1:
    input:
        "genomes/.fna.gz",
    output:
        ".fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  
            --name-from-first
    """

Lo que esto hace es decirle a snakemake que cada vez que desee un archivo de salida que termine con .fna.gz.sigdebe buscar un archivo con ese prefijo (el texto anterior .fna.gz.sig) en el genomes/ directorio, terminando en
.fna.gzy si existierautilice ese archivo como entrada.

(¡Sí, puede haber varios comodines en una regla! ¡Te lo mostraremos más adelante!)

Si revisa y usa los comodines en sketch_genomes_2 y
sketch_genomes_3notarás que las reglas terminan pareciéndose idéntico. Y resulta que solo necesita (y de hecho solo puede tener) una regla: ahora puede unir las tres reglas en una sketch_genome gobernar de nuevo.

Aquí está el completo Snakefile:

rule sketch_genome:
    input:
        "genomes/.fna.gz",
    output:
        ".fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  --name-from-first
    """

rule compare_genomes:
    input:
        "GCF_000017325.1.fna.gz.sig",
        "GCF_000020225.1.fna.gz.sig",
        "GCF_000021665.1.fna.gz.sig"
    output:
        "compare.mat"
    shell: """
        sourmash compare  -o 
    """

rule plot_comparison:
    message: "compare all input genomes using sourmash"
    input:
        "compare.mat"
    output:
        "compare.mat.matrix.png"
    shell: """
        sourmash plot 
    """

Se parece mucho al Snakefile con el que comenzamos, con la diferencia crucial de que ahora usamos comodines.

Aquí, a diferencia de la situación en la que nos encontrábamos al final de la última sección donde teníamos una regla que dibujaba tres genomas, ahora tenemos una regla que dibuja un genoma a la vez, ¡pero también se puede ejecutar en paralelo! Asi que snakemake -j 3 seguirá funcionando! Y continuará funcionando a medida que agregue más genomas y aumente la cantidad de trabajos que desea ejecutar al mismo tiempo.

Antes de hacer eso, echemos otro vistazo al flujo de trabajo ahora: notará que tiene la misma forma, ¡pero se ve ligeramente diferente! Ahora, en lugar de que las tres reglas para esbozar genomas tengan nombres diferentes, todas tienen el mismo nombre pero valores diferentes para el accession ¡comodín!

gráfico interm3 de trabajos

Capítulo 7: dar nombres de archivos de creación de serpientes en lugar de nombres de reglas

Agreguemos un nuevo genoma a la mezcla y comencemos generando un archivo de boceto (que termina en .sig) para ello.

Descargue el archivo de ensamblaje de RefSeq (el _genomic.fna.gz archivo) para GCF_008423265.1 de este enlace NCBIy colóquelo en el genomes/ subdirectorio como GCF_008423265.1.fna.gz. (También puede descargar una copia guardada con el nombre correcto de este enlace osf.io.)

Ahora, nos gustaría construir un boceto ejecutando sourmash sketch dna
(a través de la marca de serpiente).

¿Necesitamos añadir algo a la Snakefile ¿para hacer esto? ¡No, no, no lo hacemos!

Para crear un boceto para este nuevo genoma, puede pedirle a Snakemake que haga el nombre de archivo correcto así:

snakemake -j 1 GCF_008423265.1.fna.gz.sig

¿Por qué funciona esto? Funciona porque tenemos una regla de comodín genérica para construir .sig archivos de archivos en genomes/!

Cuando le pide a Snakemake que cree ese nombre de archivo, busca sus reglas en todos los bloques de salida y elige la regla con la salida coincidente; lo que es más importante, esta regla puede tiene comodines, y si los tiene, ¡extrae el comodín del nombre del archivo!

Advertencia: el sketch_genome ¡la regla ahora ha cambiado!

Como nota al margen, ya no puede pedirle a Snakemake que ejecute la regla por su nombre, sketch_genome – esto se debe a que la regla debe completar el comodín y no puede saber qué debería estar sin que le demos el nombre de archivo.

Si intentas correr snakemake -j 1 sketch_genomeobtendrá el siguiente error:

WorkflowError: las reglas de destino pueden no contener comodines. Especifique archivos concretos o una regla sin comodines en la línea de comando, o tenga una regla sin comodines en la parte superior de su flujo de trabajo (por ejemplo, la típica «regla para todos» que solo recopila todos los resultados que desea generar al final).

Esto le dice que snakemake no sabe cómo completar el comodín (y le brinda algunas sugerencias sobre cómo podría hacerlo, que exploraremos a continuación).

¡En este capítulo no necesitamos modificar el Snakefile en absoluto para hacer uso de la nueva funcionalidad!

Capítulo 8 – agregando nuevos genomas

Así que tenemos un nuevo genoma y podemos construir un boceto para él. Vamos a agregarlo a nuestra comparación, por lo que estamos creando una matriz de comparación para cuatro genomas, y no solo tres!

Para agregar este nuevo genoma a la comparación, todo lo que necesita hacer es agregar el boceto en el compare_genomes entrada, y SnakeMake localizará automáticamente el genoma asociado archivar y ejecutar
sketch_genome en él (como en el capítulo anterior), y luego ejecute
compare_genomes en eso. ¡Snakemake se encargará del resto!

rule sketch_genome:
    input:
        "genomes/.fna.gz",
    output:
        ".fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  --name-from-first
    """

rule compare_genomes:
    input:
        "GCF_000017325.1.fna.gz.sig",
        "GCF_000020225.1.fna.gz.sig",
        "GCF_000021665.1.fna.gz.sig",
        "GCF_008423265.1.fna.gz.sig",
    output:
        "compare.mat"
    shell: """
        sourmash compare  -o 
    """

rule plot_comparison:
    message: "compare all input genomes using sourmash"
    input:
        "compare.mat"
    output:
        "compare.mat.matrix.png"
    shell: """
        sourmash plot 
    """

Ahora cuando corres snakemake -j 3 plot_comparison obtendrás un
compare.mat.matrix.png archivo que contiene una matriz de 4×4! (Ver figura.)

Comparación de matriz 4x4 de genomas

Tenga en cuenta que el diagrama de flujo de trabajo ahora se ha ampliado para incluir también nuestro cuarto genoma.

gráfico interm3 de trabajos

Capítulo 9 – Uso expand hacer nombres de archivos

Puede notar que la lista de archivos en el compare_genomes todas las reglas comparten el mismo sufijo, y todas están construidas usando la misma regla. ¿Podemos usar eso de alguna manera?

¡Sí! Podemos usar una función llamada expand(...) y asígnele un nombre de archivo de plantilla para construir, y una lista de valores para insertar en ese nombre de archivo.

A continuación, construimos una lista de accesiones denominadas ACCESSIONSy luego usar
expand para construir la lista de archivos de entrada del formato .fna.gz.sig
de esa lista, creando un nombre de archivo para cada valor en ACCESSSIONS.

ACCESSIONS = ["GCF_000017325.1",
              "GCF_000020225.1",
              "GCF_000021665.1",
              "GCF_008423265.1"]

rule sketch_genome:
    input:
        "genomes/.fna.gz",
    output:
        ".fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  --name-from-first
    """

rule compare_genomes:
    input:
        expand(".fna.gz.sig", acc=ACCESSIONS),
    output:
        "compare.mat"
    shell: """
        sourmash compare  -o 
    """

rule plot_comparison:
    message: "compare all input genomes using sourmash"
    input:
        "compare.mat"
    output:
        "compare.mat.matrix.png"
    shell: """
        sourmash plot 
    """

Mientras que los comodines y expand usan la misma sintaxis, hacen cosas bastante diferentes.

expand genera una lista de nombres de archivo, basada en una plantilla y una lista de valores para insertar en la plantilla. Por lo general, se usa para hacer una lista de archivos que desea que Snakemake cree para usted.

Los comodines en las reglas proporcionan las reglas mediante las cuales se crearán realmente uno o más archivos. Son recetas que dicen: «cuando desee crear un archivo con un nombre que se parezca a ESTO, puede hacerlo a partir de archivos que se parezcan A ESTO, y esto es lo que debe ejecutar para que eso suceda.

expand le dice a Snakemake QUÉ quieres hacer, las reglas comodín le dicen a Snakemake CÓMO hacer esas cosas.

Capítulo 10: uso de reglas predeterminadas

El último cambio que haremos en el Snakefile para esta sección es agregar lo que se conoce como una regla predeterminada. ¿Qué es esto y por qué?

El ‘por qué’ es más fácil. Anteriormente, hemos tenido cuidado de proporcionar nombres de reglas o nombres de archivos específicos para snakemake, porque de lo contrario se ejecutará de manera predeterminada la primera regla en el archivo Snake. (No hay otra manera en la que el orden de las reglas en el archivo importe, pero Snakemake intentará ejecutar la primera regla en el archivo si no le da un nombre de regla o un nombre de archivo en la línea de comando).

Esto es menos que genial, porque es una cosa más para recordar y escribir. En general, es mejor tener lo que se llama una «regla predeterminada» que le permite simplemente ejecutar snakemake -j 1 para generar el archivo o archivos que desee.

Esto es sencillo de hacer, pero implica una sintaxis ligeramente diferente: una regla con solo un input, y sin shell ni bloques de salida. Aquí hay una regla predeterminada para nuestro archivo Snakefile que debe colocarse en el archivo como la primera regla:

rule all:
    input:
        "compare.mat.matrix.png"

Lo que dice esta regla es, «Quiero el archivo compare.mat.matrix.png.» No da ninguna instrucción sobre cómo hacer eso – ¡eso es lo que son el resto de las reglas en el archivo! – y no correr nada, porque no tiene bloque de shell, y tampoco crear nada, porque no tiene bloque de salida.

La lógica aquí es simple, si no directa: esta regla tiene éxito cuando existe esa entrada.

Si coloca eso en la parte superior de Snakefile, entonces ejecute
snakemake -j 1 Producirá compare.mat.matrix.png. Ya no necesita proporcionar un nombre de regla o un nombre de archivo en la línea de comando a menos que quiera hacer algo otro que generar ese archivo, en cuyo caso lo que pongas en la línea de comando ignorará el rule all:.

Capítulo 11 – nuestro Snakefile final – revisión y discusión

Aquí está el Snakefile final:

ACCESSIONS = ["GCF_000017325.1",
              "GCF_000020225.1",
              "GCF_000021665.1",
              "GCF_008423265.1"]

rule all:
    input:
        "compare.mat.matrix.png"

rule sketch_genome:
    input:
        "genomes/.fna.gz",
    output:
        ".fna.gz.sig",
    shell: """
        sourmash sketch dna -p k=31  --name-from-first
    """

rule compare_genomes:
    input:
        expand(".fna.gz.sig", acc=ACCESSIONS),
    output:
        "compare.mat"
    shell: """
        sourmash compare  -o 
    """

rule plot_comparison:
    message: "compare all input genomes using sourmash"
    input:
        "compare.mat"
    output:
        "compare.mat.matrix.png"
    shell: """
        sourmash plot 
    """

Esto Snakefile proporciona algunas características agradables.

Primero, es fácil agregar nuevos genomas a la comparación: descargamos el genoma, le ponemos el nombre de su acceso y lo agregamos a ACCESSIONS en la cima. ¡Voila!

En segundo lugar, no tenemos que recordar los nombres de ninguna regla para ejecutar todo el flujo de trabajo, porque el rule all: en la parte superior proporciona un valor predeterminado razonable.

En tercer lugar, es fácil cambiar los parámetros de boceto o comparación y luego volver a ejecutar todo el flujo de trabajo desde cero, lo que nos permite explorar rápidamente parámetros alternativos para bocetos y comparaciones si así lo deseamos.

En secciones futuras, revisaremos este Snakefile básico desde arriba y exploraremos algunos de los detalles de las reglas, los comodines y otras funciones.

Fuente del artículo

Deja un comentario