clustal one-liner con parallel

Hola,
hoy comparto un comando que a veces utilizo cuando necesito calcular muchos alineamientos múltiples a partir de una colección de ficheros de secuencias en formato FASTA. Como mi máquina, igual que la de casi todos, tiene amplia RAM y muchos cores, es un trabajo ideal para parallel. Supongamos que los archivos de salida están en la carpeta 'entrada' y queremos guardar los ficheros de salida en la carpeta 'path/to/salida', y que tenemos 20 cores disponibles:

$ mkdir /path/to/salida/

$ cd entrada

$ ls -1 *fasta | parallel --gnu -j 20 ~/soft/clustal-omega-1.2.1/src/clustalo \

--threads=1 -i {} -o /path/to/salida/{} :::

Este comando pondrá a trabajar 20 cores del sistema hasta que todos los archivos FASTA de la carpeta entrada estén alineados, con ganancias de tiempo de ejecución importantes en un experimento con 100 ficheros:

| cores (-j) | time(real) | time(user) | time(sys) |
|   1        | 4m34.440s  | 4m5.180s   | 0m2.168s  |
|  10        | 0m29.358s  | 3m57.768s  | 0m2.400s  |
|  20        | 0m23.248s  | 5m6.204s   | 0m3.364s  |

Un saludo,
Bruno

mapeo fino de genes por NGS

Buenas,
esta semana copio aquí una reseña de un trabajo recientemente publicado de Carlos P Cantalapiedra, autor habitual de este blog y próximo doctor del grupo, donde se explica el proceso para localizar un loci responsable de una resistencia a infección por parte de hongos, combinando genética clásica y secuenciación de nueva generación: http://www.eead.csic.es/spreading/showspreading?Id=416

Pongo aquí una de las figuras del artículo:

Genotipo de varias líneas de cebada en torno al locus que confiere resistencia. En naranja, genotipos como los del parental resistente. En verde, genotipos como los del parental susceptible. La captura de exoma permite reducir la zona de búsqueda al punto donde se unen ambos genotipos (punto 211721 dentro del recuadro). Adaptada de https://dl.sciencesocieties.org/publications/tpg/first-look/pdf/plantgenome2015.10.0101.pdf.

La referencia del artículo completo, en inglés, es:

Cantalapiedra CP, Contreras-Moreira B, Silvar C, Perovic D, Ordon F, Gracia MP, Igartua E, Casas A. (2016) A cluster of NBS-LRR genes resides in a barley powdery mildew resistance QTL on 7HL. The Plant Genome. Early access. DOI: 10.3835/plantgenome2015.10.0101. URL.

Hasta luego,

Bruno

Calculando experimentos de secuenciación

Buenas,
hoy necesitábamos calcular cuántos individuos (de una especie monocotiledónea) podríamos secuenciar con cierta profundidad en un secuenciador Illumina, pensando en el HiSeq2500 en concreto. Al final decidimos apostar por una profundidad promedio de 80x, para is sobre seguro, usando parejas de lecturas de 2x125b. Buscando en Internet encontré rápidamente una calculadora del propio fabricante que igual algunos no conocéis y puede ayudar a hacer esto rápidamente.

Figura tomada de http://www.danielecook.com/calculate-depth-coverage-bam-file.

Vayamos con un ejemplo con la calculadora
[ http://support.illumina.com/downloads/sequencing_coverage_calculator.html ]:

0. tipo de secuenciación: DNA             [se puede elegir RNA también]
1. protocolo: whole-genome sequencing  [otras: Nextera, Truseq, custom]
2. tamaño del genoma: 320Mbp
3. profundidad deseada: 80x
4. % de duplicados: 2%                [valor por defecto]
5. instrumento: HiSeq 1500/2500

Volumen total de secuenciación requerido: 26,1Gb   [26.122.448.980b]

En mi ejemplo, usando el protocolo v4, esto equivale a 0.42 líneas o lanes, lo que significa que podría poner hasta 2 muestras por línea.

Hasta luego,
Bruno

PD Me dicen mis colegas Dave Des Marais y Pat Edger que la longitud de un genoma (de plantas en este caso)  puede estimarse aproximadamente a partir del contenido en DNA del núclo usando la fórmula long = 1C * 980.