Uso de snakemake para realizar operaciones simples con comodines en muchos, muchos, muchos archivos

Recientemente co-enseñé otra lección de hacer serpientes (con la Dra. Abhijna Parigi), y me acordé de uno de mis usos favoritos fuera de etiqueta de la creación de serpientes: reemplazar complicado bash for bucles con flujos de trabajo de creación de serpientes simples y robustos.

Un ejemplo

Como investigador de bioinformática, con frecuencia necesito realizar operaciones simples en muchos archivos. Como parte de esto, normalmente quiero cambiar el nombre del archivo para representar el cambio en el contenido del archivo.

Por ejemplo, supongamos que tengo un montón de archivos FASTQ (digamos, los que aquí), y quiero subdividirlos en las primeras 400 líneas. Todos los nombres de archivo tienen la forma NAME.fastqy quiero agregar .subset.fastq al final de los nombres de archivo del subconjunto para distinguirlos. (Ver esta lección de scripts de shell para obtener más antecedentes y motivación para esta operación en particular).

Usando bashla ronda 1

Durante muchos años hice esto con bash for bucles El siguiente código funciona, asumiendo que los archivos fastq originales están en un data/ subdirectorio:

:::bash
mkdir subset
for i in data/*.fastq
do
    head -400 $i > subset/$(basename $i).subset.fastq
done

A partir de un montón de archivos,

>data/F3D0_S188_L001_R1_001.fastq
>data/F3D0_S188_L001_R2_001.fastq
>...

este bucle producirá

>subset/F3D0_S188_L001_R1_001.fastq.subset.fastq
>subset/F3D0_S188_L001_R2_001.fastq.subset.fastq

Mejorando la solución bash

Los nombres de los archivos de salida son un poco feos, porque fastq se repite. Eso es solo porque bash hace que sea tan fácil agregar nombres de archivos; podemos solucionar esto agregando .fastq en el $(basename ...) llamar:

:::bash
mkdir subset2
for i in data/*.fastq
do
    head -400 $i > subset2/$(basename $i .fastq).subset.fastq
done

Así que… no es difícil de leer y es bastante sencillo. ¿Por qué usaría algo más?

Usando snakemake en su lugar

tl; dr El código bash anterior es frágil cuando lo modifico; no es lo suficientemente robusto para un trabajo importante.

En mi (extensa ) experiencia, el enfoque anterior falla en un porcentaje razonable de las veces. Por lo general, lo hago bien la primera vez que lo escribo, y luego lo modifico y retoco, y se produce el caos porque omito algo en el ciclo for.

Entonces, hace un año o dos, decidí probar la fabricación de serpientes para una de estas operaciones.

Aquí está el contenido de un archivo llamado Snakefile.subsetque hace lo mismo que el bucle for anterior:

:::python
# pull in all files with .fastq on the end in the #data
FILES = glob_wildcards('data/.fastq')

# extract the  values into a list
NAMES = FILES.name

rule all:
    input:
        # use the extracted name values to build new filenames
        expand("subset3/.subset.fastq", name=NAMES)

rule subset:
    input:
        "data/.fastq"
    output:
        "subset3/.subset.fastq"
    shell: """
        head -400  > 
    """

y puedes ejecutarlo con snakemake -s Snakefile.subset -j 1.

Con este Snakefile, snakemake extrae todos los archivos que coinciden con el patrón global y extrae sus nombres, y luego construye un conjunto de “objetivos” en la regla. all que debe crear. El subset regla especifica cómo construir objetivos de ese nombre.

¿Por qué me gusta más SnakeMake que Bash para esto?

Entonces, ¿por qué me gusta más hacer serpientes? Algunas razones que creo que son intrínsecas a snakemake vs bash:

• snakemake falla si algo no está bien con los nombres de los archivos, antes ¡haciendo algo!
• Si alguna de las operaciones falla, ¡Snakemake se detiene y me alerta de manera predeterminada!
• Puedo hacer las operaciones en paralelo especificando por ejemplo snakemake -j 4 para usar 4 núcleos.
• el lenguaje de plantillas para usar es agradable, simple y estándar de Python (ver esta publicación de blog sobre f-strings y también el la referencia del minilenguaje de plantillas).
• a medida que las operaciones se vuelven más complicadas, la creación de serpientes no necesita complicarse más, mientras que la solución bash tiende a complicarse hasta convertirse en ilegibilidad…
• ¡Creo que la solución de fabricación de serpientes es más fácil de entender y modificar!

Por encima de todo, la estructura general de snakemake es declarativo en vez de procesal. Declaramos cómo queremos que se vea el resultado, y Snakemake usa las reglas disponibles para crear el conjunto general de pasos que deben ejecutarse y lo hace posible. Esto es lo que hace posible la comprobación de errores y la paralelización.

Otra “característica” de esta solución es que hay más comentarios porque comento Snakefiles más que bash scripts. Este es probablemente un problema mío causado por la creación de serpientes. forzando yo para editar un archivo :).

No he reutilizado mucho Snakefiles, pero creo que puede reutilizar Snakefiles con bastante facilidad; consulte la siguiente sección.

¿Hay alguna desventaja? La principal es que la solución de creación de serpientes me parece más pesada: implica crear un archivo, obtener el espacio/la sangría correctos, etc., etc. Así que todavía no lo uso tanto como probablemente debería.

Pensamientos bienvenidos!

–tito

Apéndice: un Snakefile más reutilizable

A continuación hay un Snakefile que es un poco más reutilizable para situaciones en las que sus directorios de entrada y salida no coinciden con los nombres que usé anteriormente; puede anular PREFIX y OUTPUT ejecutando snakemake -C prefix=PREFIX output=OUTPUT.

(Realmente no me gusta la sintaxis de usar f-strings aquí, pero es más limpio que cualquier otra cosa que haya encontrado. Se aceptan sugerencias).

:::python
# pull in all files with .fastq on the end in the 'data' directory.             
PREFIX = config.get('prefix', 'data')
print(f"looking for FASTQ files under ''/")

OUTPUT = config.get('output', 'subset5')
print(f"subset results will go under ''/")

FILES = glob_wildcards(f'/.fastq')

# extract the  values into a list                                         
NAMES = FILES.name

# request the output files                                                      
rule all:
    input:
        # use the extracted 'name' values to build new filenames                
        expand("/.subset.fastq", output=OUTPUT, name=NAMES)

# actually do the subsetting                                                    
rule subset_wc:
    input:
        f"/.fastq"
    output:
        "/.subset.fastq"
    shell: """                                                                  
        head -400  >                                             
    """