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