¡Snakemake para hacer bioinformática: entradas y salidas y más!

Como vimos antes, snakemake automáticamente “encadenará” reglas conectando entradas a salidas. Es decir, Snakemake se dará cuenta que correr para producir el resultado deseado, incluso si toma muchos pasos.

También vimos que snakemake completará y en el comando de shell basado en el contenido del input: y output: bloques Esto se vuelve aún más útil cuando se usan comodines para generalizar reglas, donde los valores de comodines se sustituyen correctamente en el y valores.

Los bloques de entrada y salida son componentes clave de los flujos de trabajo de snakemake. A continuación, discutiremos el uso de bloques de entrada y salida de manera un poco más completa.

Proporcionar entradas y salidas

Como vimos anteriormente, snakemake felizmente tomará múltiples valores de entrada y salida a través de listas separadas por comas y los sustituirá en cadenas en bloques de shell.

rule example:
   input:
       "file1.txt",
       "file2.txt",
   output:
       "output file1.txt",
       "output file2.txt",
   shell: """
       echo 
       echo 
       touch 
   """

Cuando estos se sustituyen en comandos de shell con y
se convertirán en listas ordenadas separadas por espacios: por ejemplo, el comando de shell anterior se imprimirá primero file1.txt
file2.txt y luego output file1.txt output file2.txt antes de usar touch para crear los archivos de salida vacíos.

En este ejemplo, también le pedimos a snakemake que cite los nombres de archivo para el comando de shell usando :q – esto significa que si hay espacios, caracteres como comillas simples o dobles, u otros caracteres con un significado especial, se escaparán correctamente usando
Función shlex.quote de Python. Por ejemplo, aquí ambos archivos de salida contienen un espacio, y así touch
crearía tres archivos — output, file1.txty
file2.txt — en lugar de los dos archivos correctos, output file1.txt
y output file2.txt.

Citando nombres de archivo con siempre debe usarse para cualquier cosa ejecutada en un bloque de shell – ¡No hace daño y puede prevenir errores graves!

Digresión: ¿Dónde podemos (y debemos) poner comas?

En el ejemplo de código anterior, notará que "file2.txt" y
"output file2.txt" tener comas después de ellos:

rule example:
   input:
       "file1.txt",
       "file2.txt",
   output:
       "output file1.txt",
       "output file2.txt",
   shell: """
       echo 
       echo 
       touch 
   """

¿Son necesarios? No. El código anterior es equivalente a:

rule example:
   input:
       "file1.txt",
       "file2.txt"
   output:
       "output file1.txt",
       "output file2.txt"
   shell: """
       echo 
       echo 
       touch 
   """

donde no hay comas después de la última línea en entrada y salida.

La regla general es esta: necesita comas internas para separar los elementos de la lista, porque de lo contrario las cadenas se concatenarán entre sí, es decir "file1.txt" "file2.txt" se convertirá "file1.txtfile2.txt", ¡incluso si hay una nueva línea entre ellos! Pero una coma después del último nombre de archivo es opcional (y se ignora).

¿¡Por qué!? Estos son tuplas de Python y puede agregar una coma final si lo desea: a, b, c, es equivalente a a, b, c.

Entonces, ¿por qué añadimos una coma final? Sugiero usar comas finales porque facilita agregar una nueva entrada o salida sin olvidar agregar una coma, ¡y este es un error que cometo con frecuencia! Este es un (pequeño y simple pero aún útil) ejemplo de programación defensivadonde podemos usar reglas de sintaxis opcionales para evitar errores comunes.

Las entradas y salidas son listas ordenadas

También podemos referirnos a entradas de entrada y salida individuales usando corchetes para indexarlas como listas, comenzando con la posición 0:

rule example:
   ...
   shell: """
       echo first input is 
       echo second input is 
       echo first output is 
       echo second output is 
       touch 
   """

Sin embargo, no recomendamos esto porque es frágil. Si cambia el orden de las entradas y salidas, o agrega nuevas entradas, debe revisar y ajustar los índices para que coincidan. ¡Confiar en el número y la posición de los índices en una lista es propenso a errores y hará que su Snakefile sea más difícil de cambiar más adelante!

Uso de palabras clave para archivos de entrada y salida

También puede nombrar entradas y salidas específicas usando el palabra clave
sintaxis, y luego refiérase a aquellos que usan input. y output. prefijos La siguiente regla de Snakefile hace esto:

rule example:
   input:
       a="file1.txt",
       b="file2.txt",
   output:
       a="output file1.txt",
       c="output file2.txt"
   shell: """
       echo first input is 
       echo second input is 
       echo first output is 
       echo second output is 
       touch 
   """

Aquí, a y b en el bloque de entrada, y a y c en el bloque de salida, hay nombres de palabras clave para los archivos de entrada y salida; en el comando de shell, se puede hacer referencia a ellos con , , y
respectivamente. Se puede usar cualquier nombre de variable válido, y el mismo nombre se puede usar en los bloques de entrada y salida sin colisión, como con input.a y output.aarriba, que son valores distintos.

Esta es nuestra forma recomendada de referirse a archivos de entrada y salida específicos. Es más claro de leer, resistente a las reorganizaciones o adiciones y (quizás lo más importante) puede ayudar a guiar al lector (incluido el “futuro usted”) al objetivo de cada entrada y salida.

Si usa nombres de palabras clave incorrectos en su código de shell, obtendrá un mensaje de error. Por ejemplo, este código:

rule example:
   input:
       a="file1.txt",
   output:
       a="output file1.txt",
   shell: """
       echo first input is 
   """

da este mensaje de error:

AttributeError: 'InputFiles' object has no attribute 'z', when formatting the following:

       echo first input is 

Ejemplo: escribir una línea de comando flexible

Un ejemplo en el que es particularmente útil poder hacer referencia a entradas específicas es cuando se ejecutan programas en archivos donde los nombres de los archivos de entrada deben especificarse como argumentos opcionales. Uno de esos programas es el megahit ensamblador cuando se ejecuta en lecturas de entrada emparejadas. Considere el siguiente Snakefile:

rule all:
    input:
        "assembly_out"

rule assemble:
    input:
        R1="sample_R1.fastq.gz",
        R2="sample_R2.fastq.gz",
    output:
        directory("assembly_out")
    shell: """
        megahit -1  -2  -o 
    """

En el comando de shell aquí, necesitamos proporcionar las lecturas de entrada como dos archivos separados, con -1 ante uno y -2 antes del segundo. ¡Como beneficio adicional, el comando de shell resultante es muy legible!

Funciones de entrada y características más avanzadas

Hay una serie de usos más avanzados de entrada y salida que se basan en la programación de Python; por ejemplo, uno puede definir una función de Python que se llama para generar un valor dinámicamente, como se muestra a continuación:

def multiply_by_5(w):
    return f"file.txt"


rule make_file:
    input:
        # look for input file.txt if asked to create output.txt
        filename=multiply_by_5,
    output:
        "output.txt"
    shell: """
        cp  
    """

Cuando se le pide que cree output5.txtesta regla buscará
file25.txt como entrada.

Dado que esta funcionalidad se basa en el conocimiento de los comodines, así como en algún conocimiento de Python, es demasiado avanzada para hablar de ella aquí. ¡Más sobre eso más tarde!

Referencias y enlaces

Como vimos anteriormente, los bloques de entrada y salida son clave para la forma en que funciona SnakeMake: permiten que SnakeMake conecte automáticamente las reglas en función de las entradas necesarias para crear la salida deseada. Sin embargo, los bloques de entrada y salida están limitados de ciertas maneras: más específicamente, cada entrada en los bloques de entrada y salida debe ser un nombre de archivo. Y, debido a la forma en que funciona SnakeMake, los nombres de archivo especificados en los bloques de entrada y salida deben existir para que el flujo de trabajo supere esa regla.

Con frecuencia, los comandos de shell necesitan tomar parámetros distintos a los nombres de archivo, y estos parámetros pueden ser valores que puede o debe calcular Snakemake. Por lo tanto, snakemake también admite una
params: bloque que se puede utilizar para proporcionar cadenas de parámetros que son no
nombres de archivo en el bloque de shell. Como verá a continuación, estos se pueden usar para una variedad de propósitos, incluidos los parámetros configurables por el usuario, así como los parámetros que el código de Python puede calcular automáticamente.

Un ejemplo simple de un bloque de parámetros

Considerar:

rule use_params:
    params:
        val = 5
    output: "output.txt"
    shell: """
        echo  > 
    """

Aquí, el valor 5 se le asigna el nombre val en el params: bloque, y luego está disponible bajo el nombre en el shell: bloquear. Esto es análogo al uso de palabras clave en bloques de entrada y salida, pero a diferencia de los bloques de entrada y salida, las palabras clave debe ser utilizado en bloques de parámetros.

En este ejemplo, no se gana en funcionalidad, pero se gana en legibilidad: la sintaxis deja claro que val es un parámetro ajustable que se puede modificar sin comprender los detalles del bloque de shell.

Los bloques de parámetros tienen acceso a comodines

Al igual que el input: y output: bloques, los valores comodín están directamente disponibles en params: bloques sin usar el wildcards
prefijo; por ejemplo, esto significa que puede usarlos en cadenas con el estándar operaciones de formateo de cadenas.

Esto es útil cuando un comando de shell necesita usar algo que no sea el nombre del archivo, por ejemplo, el bowtie El software de alineación de lectura toma la prefijo del archivo SAM de salida a través de -Slo que significa que no puede nombrar el archivo correctamente con bowtie ... -S . En su lugar, podrías usar al igual que:

rule all:
    input:
        "reads.sam"

rule use_params:
    input: ".fq",
    output: ".sam",
    params:
        prefix = ""
    shell: """
        bowtie index -U  -S 
    """

Si fueras a usar -S aquí, terminarías produciendo un archivo
reads.sam.sam!

Enlaces y referencias:

Los comodines de Snakemake facilitan la aplicación de reglas a muchos archivos, pero también crean un nuevo desafío: ¿cómo generar todos los nombres de archivo que desea?

Como ejemplo de este desafío, considere la lista de genomas necesarios para la regla compare_genomes desde antes –

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",

Esta lista es crítica. porque especifica los bocetos que se crearán mediante la regla de comodines. Sin embargo, escribir esta lista es molesto y propenso a errores, porque partes de cada nombre de archivo son idénticas y se repiten.

Peor aún, si necesita usar esta lista en varios lugares, o producir nombres de archivo ligeramente diferentes con las mismas accesiones, será propenso a errores: es probable que desee agregar, eliminar o editar elementos de la lista, y lo hará. necesita cambiarlo en varios lugares.

Anteriormente, mostramos cómo cambiar esto a una lista de las accesiones en la parte superior de Snakefile y luego usamos una función llamada
expand para generar la lista:

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

#...

rule compare_genomes:
    input:
        expand(".fna.gz.sig", acc=ACCESSIONS),

Usando expand generar listas de nombres de archivo es un patrón común en Snakefiles, ¡y lo exploraremos más a continuación!

Usando expand con un solo patrón y una lista de valores

En el ejemplo anterior, proporcionamos un solo patrón, .fna.gz.sigy pregunta expand para resolverlo en muchos nombres de archivo completando valores para el nombre del campo acc de cada elemento en ACCESSIONS. (Es posible que reconozca la sintaxis de palabras clave para especificar valores, acc=ACCESSIONSde los bloques de entrada y salida, arriba!

El resultado de expand('.fna.gz.sig', acc=...) aquí está
idéntico para escribir los cuatro nombres de archivo en forma larga:

"GCF_000017325.1.fna.gz.sig",
"GCF_000020225.1.fna.gz.sig",
"GCF_000021665.1.fna.gz.sig",
"GCF_008423265.1.fna.gz.sig"

Eso es, expand no hace ninguna coincidencia especial de comodines o inferencia de patrones, solo completa los valores y devuelve la lista resultante.

Aquí, ACCESSIONS puede ser cualquier Python iterable – por ejemplo, una lista, una tupla o un diccionario.

Usando expand con múltiples listas de valores

También puedes usar expand con varios nombres de campo. Considerar:

expand('.fna.`, acc=ACCESSIONS, extension=['.gz.sig', .gz'])

Esto producirá los siguientes ocho nombres de archivo:

"GCF_000017325.1.fna.gz.sig",
"GCF_000017325.1.fna.gz",
"GCF_000020225.1.fna.gz.sig",
"GCF_000020225.1.fna.gz",
"GCF_000021665.1.fna.gz.sig",
"GCF_000021665.1.fna.gz",
"GCF_008423265.1.fna.gz.sig",
"GCF_008423265.1.fna.gz"

sustituyendo todo posible combinaciones de acc y extension en el patrón proporcionado.

generando todo combinaciones contra por parejas combinaciones

Como vimos arriba, con múltiples patrones, expand generará todas las combinaciones posibles: es decir,

X = [1, 2, 3]
Y = ['a', 'b', 'c']

rule all:
   input:
      expand('.by.', x=X, y=Y)

generará 9 nombres de archivo: 1.by.a, 1.by.b, 1.by.c, 2.by.aetc. Y si agregaste un tercer patrón al expand cadena, expand ¡También agregaría eso a las combinaciones!

Entonces, ¿qué está pasando aquí?

Por defecto, expand hace una expansión de todo por todo que contiene todas las combinaciones posibles. (Esto a veces se denomina producto cartesiano, producto cruzado o combinación externa).

Pero no siempre quieres eso. ¿Cómo podemos cambiar este comportamiento?

El expand función toma un segundo argumento opcional, el combinador, que le dice expand cómo combinar las listas de valores que vienen después. Por defecto expand utiliza una función de Python llamada
itertools.productque crea todas las combinaciones posibles, pero puedes darle otras funciones.

En particular, puedes decir expand para crear combinaciones por pares usando zip en cambio, algo que hicimos en uno de los ejemplos de comodines.

Aquí hay un ejemplo:

X = [1, 2, 3]
Y = ['a', 'b', 'c']

rule all:
   input:
      expand('.by.', zip, x=X, y=Y)

que ahora generará solo tres nombres de archivo: 1.by.a, 2.by.by 3.by.c.

La gran advertencia aquí es que zip creará una lista de salida de la longitud de la lista de entrada más corta, por lo que si le da una lista de tres elementos y una lista de dos elementos, solo usará dos elementos de la primera lista.

por ejemplo, en el expand en esto Snakefile,

X = [1, 2, 3]
Y = ['a', 'b']

rule all_zip_short:
   input:
      expand('.by.', zip, x=X, y=Y)

solo 1.by.a y 2.by.b se generará, ya que no hay ningún socio para 3 en la segunda lista.

Para obtener más información, consulte el documentación de snakemake sobre el uso de zip en lugar de producto.

Obtener una lista de identificadores para usar en expand

El expand proporciona una solución efectiva cuando tiene listas de identificadores que usa varias veces en un flujo de trabajo, ¡un patrón común en bioinformática! Sin embargo, escribir estas listas en un Snakefile (como lo hacemos en los ejemplos anteriores) no siempre es práctico; ¡puede tener de docenas a cientos de identificadores!

Las listas de identificadores se pueden cargar desde otro archivos en una variedad de formas, y también se pueden generar a partir del conjunto de archivos reales en un directorio usando glob_wildcards.

Ejemplos de carga de listas de accesiones desde archivos o directorios

Cargar una lista de accesiones desde un archivo de texto

Si tiene una lista simple de accesiones en un archivo de texto
accessions.txtal igual que:

Archivo accessions.txt:

GCF_000017325.1
GCF_000020225.1
GCF_000021665.1
GCF_008423265.1

luego, el siguiente código cargará cada línea en el archivo de texto como una identificación separada.

Snakefile para cargar accessions.txt:

with open('accessions.txt', 'rt') as fp:
    ACCESSIONS = fp.readlines()
    ACCESSIONS = [ line.strip() for line in ACCESSIONS ]

print(f'ACCESSIONS is a Python list of length ')
print(ACCESSIONS)

rule all:
    input:
        expand(".sig", acc=ACCESSIONS)

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

y construya firmas sourmash para cada accesión.

Cargar una columna específica desde un archivo CSV

Si en lugar de un archivo de texto tiene un archivo CSV con varias columnas y los ID para cargar están todos en una columna, puede usar Python
biblioteca de pandas para leer en el CSV. En el código de abajo, pandas.read_csv carga el CSV en un objeto pandas DataFrame, y luego seleccionamos el accession columna y usar eso como un iterable.

Archivo accessions.csv:

accession,information
GCF_000017325.1,genome 1
GCF_000020225.1,genome 2
GCF_000021665.1,genome 3
GCF_008423265.1,genome 4

Snakefile para cargar accessions.csv:

import pandas

CSV_DATAFRAME = pandas.read_csv('accessions.csv')
ACCESSIONS = CSV_DATAFRAME['accession']

print(f'ACCESSIONS is a pandas Series of length ')
print(ACCESSIONS)

rule all:
    input:
        expand(".sig", acc=ACCESSIONS)

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

Cargando desde el archivo de configuración

Snakemake también admite el uso de archivos de configuración, donde el archivo de serpiente proporciona el nombre de un archivo de configuración predeterminado (que a su vez puede anularse en la línea de comando).

Un archivo de configuración también puede ser un buen lugar para colocar las accesiones. Considerar:

accessions:
- GCF_000017325.1
- GCF_000020225.1
- GCF_000021665.1
- GCF_008423265.1

que es utilizado por el siguiente Snakefile.

Snakefile para cargar accesiones desde config.yml:

configfile: "config.yml"

ACCESSIONS = config['accessions']

print(f'ACCESSIONS is a Python list of length ')
print(ACCESSIONS)

rule all:
    input:
        expand(".sig", acc=ACCESSIONS)

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

Aquí, config.yml es un archivo YAML, que es un formato legible por humanos que también puede ser leído por computadoras. ¡Hablaremos de los archivos de configuración más tarde!

Usando glob_wildcards para cargar identificaciones o accesiones desde un conjunto de archivos

presentamos el glob_wildcards Comando brevemente en la publicación sobre comodines: glob_wildcards hace coincidencia de patrones en los archivos realmente presente en el directorio.

Aquí hay un Snakefile que usa glob_wildcards para obtener las cuatro accesiones de los nombres de archivo reales:

GLOB_RESULTS = glob_wildcards("genomes/.fna.gz")
ACCESSIONS = GLOB_RESULTS.acc

print(f'ACCESSIONS is a Python list of length ')
print(ACCESSIONS)

rule all:
    input:
        expand(".sig", acc=ACCESSIONS)

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

Esta es una forma particularmente conveniente de obtener una lista de accesiones, aunque puede ser peligroso utilizarla. En particular, es fácil eliminar accidentalmente un archivo y no notar que falta una muestra. Por ese motivo, sugerimos proporcionar una lista independiente de archivos para cargar en muchas situaciones.

Comodines y expand – algunos pensamientos finales

Combinado con comodines, expand es extremadamente poderoso y útil. Sin embargo, al igual que los comodines, este poder viene con cierta complejidad. Aquí hay un breve resumen de cómo se combinan estas características.

El expand función hace un lista de archivos para crear de un patrón y una lista de valores para completar.

Los comodines en las reglas proporcionan recetas para crear archivos cuyos nombres coincidan con un patrón.

Típicamente en Snakefiles usamos expand para generar una lista de archivos que coincidan con un determinado patrón y luego escriba una regla que use comodines para generar esos archivos reales.

La lista de valores a usar con expand puede provenir de muchos lugares, incluidos archivos de texto, archivos CSV y archivos de configuración. Puede también viene de
glob_wildcardsque utiliza un patrón para extracto la lista de valores de los archivos que están realmente presentes.

Enlaces y referencias