instalación de Grid Engine en servidor multi-core Debian

Hola,

esta es mi primera entrada tras volver a la EEAD-CSIC.

Si hace unos años explicaba en este blog cómo instalar el gestor de colas de cálculo gridengine en un servidor  Ubuntu, esta semana lo actualizo para un servidor multicore con Debian 11:

Instalación de los diferentes componentes gridengine en la misma máquina:

   sudo apt install gridengine-master gridengine-qmon gridengine-exec

Si tuvieras que empezar de cero puedes hacerlo con:

    sudo apt purge gridengine-common
    sudo rm -rf /var/spool/gridengine/spooldb/sge
    sudo apt install gridengine-master gridengine-qmon gridengine-exec

Tras leer esto edité el fichero /etc/hosts así:

   127.0.0.1 localhost.localdomain localhost
   127.0.1.1 master master
   121.xxx.yyy.zzz myhost
 

Ojo que la instalación crea el usuario sgeadmin. Para que tu propio usuario tenga privilegios de administración, crear y configurar una cola llamada all.q deberás hacer los siguiente:  

    sudo -u sgeadmin qconf -am myuser

    # and to a userlist:
    qconf -au myuser users

   # Add a submission host:
   qconf -as myhost

   # Add an execution host, you will be prompted for information about the execution host
   qconf -ae
 
   # Add a new host group:
   qconf -ahgrp @allhosts

   # Add the exec host to the @allhosts list:
   qconf -aattr hostgroup hostlist myhost @allhosts

   # Add and configure queue, set the slots to CPU/cores, check parallel env (pe_list)
   qconf -aq all.q

   # Add the host group to the queue:
   qconf -aattr queue hostlist @allhosts  all.q

   # Allocate slots in this queue:
   qconf -aattr queue slots "[myhost=12]" all.q

Con esto estaría, y puedes comprobarlo con qstat -f o lanzando un trabajo con qsub

Con los siguientes comandos puedes reiniciar los respectivos servicios en Debian:

   sudo systemctl restart gridengine-master.service
   sudo systemctl restart gridengine-exec.service

Hasta pronto,

Bruno

Cortar secuencias desde un fichero GFF y un FASTA del genoma

Hola,

hoy cuento cómo cortar subsecuencias de un genoma correspondientes a elementos anotados en un fichero GFF asociado, que tienen este aspecto, con columnas separadas por tabuladores. A diferencia del formato BED, las coordinadas aquí empiezan a contar en 1:

1  proveedor  gene  16399  20144  .  +  .  ID=500;...
1  proveedor  mRNA  16399  20144  .  +  .  ID=500-01;Parent=500;...
1  proveedor  exon  16399  16976  .  +  .  Parent=500-01;Name=500-01-E1;...
1  proveedor  CDS   16599  16976  .  +  0  ID=CDS:500-01;Parent=500-01;...
1  proveedor  exon  17383  17474  .  +  .  Parent=500-01;Name=O500-01-E2;...
1  proveedor  CDS   17383  17474  .  +  0  ID=CDS:500-01;Parent=500-01;...
...

Como se puede ver, cada gen es un intervalo dentro de un cromosoma o contig. En este ejemplo, adaptado de Ensembl Plants, el gen con identificador 500 está en el segmento 16399-20144 del cromosoma 1 en la hebra directa (+).

Para cortar la secuencia de los genes en un fichero GFF podemos hacer lo siguiente, con ayuda de bedtools getfasta, que deberás instalar previamente. Los identificadores de cromosomas/contigs deben coincidir en ambos ficheros:

# por acortar las URLs
ENSEMBL="http://ftp.ebi.ac.uk/ensemblgenomes/pub/release-51/plants"

# descargar fichero FASTA del genoma, o del chr en este caso
wget ${ENSEMBL}/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa.gz .
gunzip Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa.gz

# descargar y descomprimir fichero GFF
wget ${ENSEMBL}/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3.gz .
gunzip Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3.gz

# extraer solamente los elementos de tipo "gene"
perl -lane 'print if($F[2] eq "gene")' Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3 > \ 
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.gff3

# cortar los genes y guardarlos en un fichero FASTA
bedtools getfasta -fi Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    -bed Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.gff3 \
    -fo Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.fa

Si necesitamos cortar las secuencias codificantes (CDS) necesitamos instalar gffread desde https://github.com/gpertea/gffread (mirar artículo aquí). Con este programa puedes extraer los CDS como nucleótidos o aminoácidos, pero tiene muchas más opciones:

# secuencia de exones codificantes concatenados, una por tránscrito/mRNA
/path/to/gffread-0.12.7.Linux_x86_64/gffread -x cds.fna \
    -g Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3

# traducción de exones codificantes concatenados, una por tránscrito/mRNA
/path/to/gffread-0.12.7.Linux_x86_64/gffread -y cds.faa \
    -g Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3

Hasta pronto,

Bruno

Análisis de secuencias eficiente con SeqAn

Hola,

hoy he participado en un curso de introducción a SeqAn, una biblioteca Open Source escrita en lenguaje C++ para el análisis eficiente de secuencias biológicas, una tarea de la que ya habíamos hablado en este blog. La versión original se publicó en BMC Bioinformatics y ha evolucionado hasta la actual SeqAn3, que puede obtenerse en https://github.com/seqan/seqan3 y require un compilador g++ con versión >= 7 y el estándar C++17  (g++-7 -std=c++17).

Para programar tus propias aplicaciones deberás tener algo de conocimiento de C++ moderno, en concreto de rangos, funciones lambda y la STL, pero afortunadamente el material del curso es muy ameno y fácil de seguir. Lo puedes encontrar en:

https://seqan.github.io/learning-resources/doc/biocpp/presentation

 

Os dejo un ejemplo que permite filtrar lecturas de secuencia de un fichero FASTQ que superen un cierto umbral de calidad y require leer la documentación de la API: 

#include <string_view>
#include <seqan3/std/algorithm>
#include <seqan3/io/sequence_file/all.hpp>
#include <seqan3/alphabet/quality/all.hpp>

//ask user a PHRED quality cutoff value
auto ask_quality() -> seqan3::phred42
{
   std::cout<<" Please type read quality cutoff [0,41]\n";
   uint32_t cutoff{};
   std::cin>>cutoff;
        
   return seqan3::assign_rank_to( cutoff, seqan3::phred42{});
}

int main(int argc, const char * argv[])
{

   //note there's no user interface help nor proper arg checking
   std::string_view fastq_path = argv[1];
   std::string_view fasta_path = argv[2];

   // note file formats are deduced from file extensions
   seqan3::sequence_file_input fastq_file_in{fastq_path};
   seqan3::sequence_file_output fasta_file_out{fasta_path};

   // filter reads by quality
   seqan3::quality_alphabet qual_cutoff = ask_quality();
   
   // this is where you need ti understand ranges & lambda functions,
   // plus a little bit about the SeqAnAPI
   for (const auto & [sequence,id,quality] : fastq_file_in)
   {
      if (std::ranges::all_of(quality, [qual_cutoff] (auto const & quality) 
          { return quality >= qual_cutoff; }))
          fasta_file_out.emplace_back(sequence, id);
   }

   return 0;
}

Si te interesa, el próximo curso gratuito será el 14 de septiembre de 16:00 a 19:00 CEST (más información en https://ogy.de/ts51). Hasta pronto,

Bruno