Mostrando entradas con la etiqueta secuencias. Mostrar todas las entradas
Mostrando entradas con la etiqueta secuencias. Mostrar todas las entradas

15 de septiembre de 2021

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

11 de febrero de 2014

Jugando a leer y extraer secuencias de un fichero FASTQ comprimido con GZIP (.fq.gz)

Un fichero con extensión '.fq.gz' es un fichero en formato FASTQ comprimido en formato GZIP. Los ficheros FASTQ con datos de los nuevos métodos de secuenciación (NGS) suelen ocupar decenas de Gigabytes de espacio de disco (una cantidad indecente para la mayoría de los mortales) y comprimidos con GZIP se reducen.

Para empezar intentaremos saber cuanto espacio ocupa realmente un archivo 'fq.gz' comprimido, para ello usaremos el mismo comando 'gzip':

> gzip --list reads.fq.gz
         compressed        uncompressed  ratio uncompressed_name
        18827926034          1431825024 -1215.0% reads.fq


Parece que GZIP en vez de comprimir expande, pero no es verdad, simplemente que la opción '--list' de 'gzip' no funciona correctamente para archivos mayores de 2GB. Así que hay que recurrir a un método más lento y esperar unos minutos:


> zcat reads.fq.gz | wc --bytes
61561367168


Si queremos echar un vistazo al contenido del archivo podemos usar el comando 'less' o 'zless':

> less reads.fq.gz
> zless reads.fq.gz

Y para saber el número de secuencias que contiene simplemente hay que contar el número de líneas y dividir por 4 (le costará unos minutos):

> zcat reads.fq.gz | echo $((`wc -l`/4))
256678360

Podemos buscar una determinada secuencia en el archivo y contar cuántas veces aparece, por ej. ATGATGATG:

> zcat reads.fq.gz | awk '/ATGATGATG/ {nlines = nlines + 1} END {print nlines}'
398065

A veces nos interesará tomar un trozo del fichero para hacer pruebas, por ejemplo las primeras 1000 secuencias (o 4000 líneas):
> zcat reads.fq.gz | head -4000 > test_reads.fq
> zcat reads.fq.gz | head -4000 | gzip > test_reads.fq.gz

O extraer un rango de líneas (1000001-1000004):

> zcat reads.fq.gz | sed -n '1000001,1000004p;1000005q' > lines.fq 

También nos puede interesar dividir el fichero en varios más pequeños de por ejemplo 1 miĺlón de secuencias (o 4 millones de líneas):

> zcat reads.fq.gz | split -d -l 4000000 - reads.split
> gzip reads.split*
> rename 's/\.split(\d+)/.$1.fq/' reads.split*

Y posteriormente reunificarlos:

> zcat reads.*.fq.gz | gzip > reads.fq.gz


Finalmente os recomiendo leer un post similar en inglés:
http://darrenjw.wordpress.com/2010/11/28/introduction-to-the-processing-of-short-read-next-generation-sequencing-data/

5 de junio de 2013

FASTQ sort + parallel

Buenas,
recientemente en el laboratorio hemos estado manipulando archivos de pares de secuencias (pair-end reads) en formato FASTQ. Una de las tareas habituales ha sido limpiar los archivos de secuencias de baja calidad, ya sea recortando o eliminando directamente, para luego volver a definir parejas entre las secuencias que superaron el corte. Una estrategia posible para esta tarea es simplemente linearizar las secuencias, de manera que cada una ocupe ahora una sola línea, separando con tabuladores la cabecera, la secuencia, el separador y las calidades. Por ejemplo, la siguiente secuencia:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTA...
+
!''*((((***+))%%%++)(...
 
quedaría así:

@SEQ_ID  GATTTGGGGTTCAAAGCAGTA...   +   !''*((((***+))%%%++)(...

Tras esta transformación ya sí es posible ordenar un archivo FASTQ con GNU sort, que viene instalado en cualquier sistema linux (y se puede instalar en Windows). GNU sort es ideal para ordenar conjuntos de datos que no caben en memoria (M), como a menudo ocurre con los archivos FASTQ, porque de manera implícita divide el problema inicial, de tamaño N, en N/M trozos que luego mezcla (merge) de manera externa.

En nuestra experiencia GNU sort es significativemente más eficiente que nuestros scripts para este tipo de problemas, puesto que ya trae de fábrica toda la lógica para partir el problema en trozos y luego mezclar las soluciones parciales. Sólo hay que tener cuidado de asignar la variable de ambiente LC_ALL, por ejemplo con:

$ export LC_ALL=POSIX

y echar a andar. Muy bien. Pero te quedas con la duda de si estás sacando el máximo partido a tu CPU multicore, podremos optimizar sort en paralelo? Y si invocamos a GNU parallel (mira el vídeo)?

Nos ponemos manos a la obra y hacemos pruebas con un archivo FASTQ linearizado real, de 576Mb:
 
$ ls -lh /tmp/unsortedXXpJ6CaB
-rw-------. 1 576M Jun  5 10:15 /tmp/unsortedXXpJ6CaB
 
Ahora lo ordenamos con GNU sort, dándole un máximo de 500Mb de área en RAM para trabajar (la M de antes):
 
$ time sort -k 1,1 -u -S 500M /tmp/unsortedXXpJ6CaB > /tmp/unsortedXXpJ6CaB.S
real 0m7.628s
user 0m5.143s
sys 0m2.373s

Finalmente probamos ahora con parallel, agrupando las secuencias en grupos de 100.000 elementos (ojo con esto, puedes llegar a obtener resultados parcialmente desordenados porque la segunda llamada a parallel puede recibir más argumentos de los que el shell soporta). Cambiando este valor a 10000 o a 10E6 los resultados son similares:

$ time cat /tmp/unsortedXXpJ6CaB | parallel -N 100000 --pipe --files sort -k 1,1 | 
   parallel -Xj1 sort -k 1,1 -u -m {} ';' rm {} > /tmp/unsortedXXpJ6CaB.P
real 0m15.451s
user 0m9.919s
sys 0m8.371s

Comprobamos que los resultados son idénticos:
 
$ diff /tmp/unsortedXXpJ6CaB.S /tmp/unsortedXXpJ6CaB.P

Conclusión de estas pruebas: no vale la pena complicarse con parallel para ordenar grandes archivos FASTQ, ya que probablemente el cuello de botella sea el merge final, y eso parece resolverlo mejor directamente GNU sort. De todos modos es posible que haya otras maneras de invocar a parallel más ventajosas
Si en vuestras pruebas obtenéis resultados distintos por favor escribid,
un saludo,
Bruno



30 de noviembre de 2012

comparando secuencias de igual longitud

Hola,
hoy quería compartir recetas en Perl para la comparación eficiente de secuencias de igual longitud. Las dos funciones que vamos a ver devuelven como resultado listas con las posiciones (contando desde 0) en las que difieren dos secuencias.

Para qué sirve esto? En general, para calcular distancias de edición (edit distances); en bioinformática es una manera, por ejemplo, de comparar oligos de igual longitud o secuencias extraídas de un alineamiento múltilple, incluyendo gaps. El siguiente código fuente incluye una subrutina en Perl (con el operador binario ^ y la función pos) y otra en C embebido para hacer la misma operación por dos caminos:

 use strict;  
 use warnings;  
 my $seq1 = 'ACTGGA';  
 my $seq2 = 'AGTG-A';  
   
 my @Pdiffs = find_diffs($seq1,$seq2);  
 my @Cdiffs = find_diffs_C($seq1,$seq2);  
   
 printf("# version Perl\ndiferencias = %d\n%s\n",scalar(@Pdiffs),join(',',@Pdiffs));  
 printf("# version C\ndiferencias = %d\n%s\n",scalar(@Cdiffs),join(',',@Cdiffs));  
   
 sub find_diffs  
 {  
    my ($seqA,$seqB) = @_;  
    my $XOR = $seqA ^ $seqB; 
    my @diff;   
    while($XOR =~ /([^\0])/g)  
    {   
       push(@diff,pos($XOR)-1);  
    }  
    return @diff;  
 }  
   
 use Inline C => <<'END_C';  
 void find_diffs_C(char* x, char* y)   
 {                        
   int i;   
    Inline_Stack_Vars;                               
   Inline_Stack_Reset;                                                                
   for(i=0; x[i] && y[i]; ++i) {                          
    if(x[i] != y[i]) {                              
     Inline_Stack_Push(sv_2mortal(newSViv(i)));                 
    }                                       
   }   
    Inline_Stack_Done;                                                                  
 }     
 END_C  
Podéis ver en stackoverflow otras versiones y su comparación en base a su tiempo de cálculo,
un saludo,
Bruno




13 de junio de 2012

Alineamiento con transformadas de Fourier

En 1992 el grupo de Ilya Vakser publicó en PNAS el método fundamental para poder hacer simulaciones de docking de biomoléculas de manera eficiente, disminuyendo considerablemente el tiempo de cálculo aplicando transformadas rápidas de Fourier. Diez años más tarde se publicó la primera versión del exitoso programa de alineamiento múltiple de secuencias MAFFT, del que ya hemos hablado en otras entradas, que aplica ideas similares para el problema de encajar secuencias similares, como se hace al alinear secuencias.

Figura tomada de http://www.ncbi.nlm.nih.gov/pubmed/12136088. donde se muestran dos picos del análisis de Fourier que se corresponden con dos subalineamientos locales.

En esta entrada mi intención es explicar el fundamento de esta técnica recurriendo al lenguaje Octave, la versión open source de MatLab. El código es el siguiente:

 % alineamiento (sin gaps) de secuencias proteicas usando la FFT  
 % escrito en Octave, que es compatible con Matlab   
 % requiere la libreria 'bioinfo'   
   
 clear all  
   
 %% Conversor binario de secuencias   
 %% Parametros:   
 %% 1) longitud de la secuencia S  
 %% 2) longitud del fragmento F (de S)  
 function [S , F] = aa2bits( sequence , fragment )  
   
 % cada residuo se representa por un numero del 1 al 20  
 % http://www.mathworks.com/help/toolbox/bioinfo/ref/aa2int.html  
 [ seq ]= aa2int( sequence )  
 [ frag ] = aa2int( fragment )   
   
 % codificamos residuos como cadenas binarias de ancho fijo (20 columnas)  
 S=[];  
 for i=1:length(seq)  
   % cada residuo ocp  
   S=[S ones(1,seq(i)) zeros(1,20-seq(i)) ];  
 end  
   
 F=[];  
 for i=1:length(frag)  
   F=[F ones(1,frag(i)) zeros(1,20-frag(i)) ];  
 end  
   
 % completa F hasta igualar la longitud de S  
 F=[F zeros(1,length(S)-length(F))];  
   
 endfunction % no se usa en Matlab  
   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   
 % define la secuencia S y el fragmento F que queremos alinear  
 sequence = 'SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY';  
 fragment = 'KMLNSVCRSWWWWC';  
   
 % convertimos ambas secuencias en dos segnales binarias,  
 % de manera que cada residuo se representa por un natural del 1 al 20  
 [S , F] = aa2bits(sequence,fragment);  
   
 % calculamos el alineamiento de F con S optimizando la funcion objetivo Fobj,  
 % que en este ejemplo es simplemente el producto (como un AND binario)  
 % obtenido al desplazar F de izq a der sobre S:   
 % producto(n) = sum { F(i) * S(i-n) }  
 %        i   
 % aprovechamos que Fobj se puede expresar por medio de transformadas  
 % de Fourier para reducir las operaciones necesarias (de N^2 a 4NlogN)  
 % https://sites.google.com/site/cartografiaygeodesia/prac7.pdf  
 % http://www.pnas.org/content/89/6/2195.abstract  
   
 FTsec = fft(S);  
 FTfrag = fft(F);  
   
 FTproducto = conj(FTfrag) .* FTsec; % producto termino a termino  
   
 % Deshacemos transformacion y localizamos valor maximo,   
 % la posicion que optimiza el alineamiento   
 producto = ifft(FTproducto);   
 [valor maxpos] = max(producto);  
   
 % imprime alineamiento sobre secuencicas binarias   
 tit = sprintf("Optimal alignment, optimal position = %d",maxpos);  
 plot(S*0.75)  
 title(tit)  
 xlabel('sequence position')  
 hold on  
 plot([zeros(1,maxpos) F*0.9] , 'r');  
 axis([0 length(S) 0 1]);  
 legend('sequence','fragment');  
 print -dpng figure_align.png;  
   
 % alineamiento sobre secuencias originales, cambiando maxpos de escala  
 align = printf("Alignment:\nS %s\nF %s%s\n",sequence,blanks((maxpos-1)/20),fragment);  

Como resultado obtendremos un alineamiento como éste, que :
S SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY
F                        KMLNSVCRSWWWWC
 

Que se corresponde con éste producto de Fourier:



En Perl podríamos usar el módulo Math::FFT para este fin, o la GNU Scientific Library que ya revisamos en otra entrada. Un saludo,
Bruno

15 de septiembre de 2011

Guía de campo de tecnologías de secuenciación

Hola,
ayer me encontré en la Red una revisión, publicada en Mayo de 2011 por TC Glenn, que contiene la siguiente tabla, muy útil  para comparar de un vistazo las plataformas de secuenciación de segunda generación disponibles actualmente:

Tabla original publicada en Molecular Ecology Resources
Esta tabla se complementa con otras disponibles en la 'NGS Field Guide', actualizadas regularmente, incluyendo por ejemplo los costes y los tipos de errores más frecuentes en cada una de ellas. De hecho habrá que esperar para tener datos empíricos de los errores típicos de la plataforma IonTorrent, que por ahora se basan en datos proporcionados por la compañía (previas a su publicación en Nature el pasado mes de Julio).
Hasta otra, Bruno


13 de abril de 2011

Bajar de ensembl todas las secuencias de proteína de los ortólogos de un gen (o de varios)

Dada una secuencia de proteína, puede interesarnos bajar de una sola vez todas las secuencias de proteínas ortólogas conocidas. Si aceptamos el criterio de ensembl para determinar qué es y qué no es ortólogo (una buena opción en una decisión a veces complicada), las siguientes instrucciones pueden servir de guía para esta tarea.

i) pasos previos

Tendréis que tener instalado el ensemble-api (ver instrucciones de instalación) y el módulo DBD::Mysql. En Mac OS X he conseguido instalar este módulo usando macports. En 5 líneas, sería:

~% sudo port selfupdate
~% sudo port install mysql5-server
~% sudo port install p5-dbd-mysql
~% sudo port install p5-dbi
~% sudo port -f activate p5-dbi



ii) bajar los ortólogos

En un fichero de texto, ponéis una lista con los identificadores de genes para los que queréis hacer la búsqueda.

(por ejemplo)
~% cat IDs.txt
ENSMUSG00000031628
ENSMUSG00000027997
ENSG00000003400

Ejecutáis get_orthologues.pl dándole este fichero como $ARGV[0]

~% perl get_orthologues.pl IDs.txt

Si veis algo así en pantalla, todo OK.

.
.
--
Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000031628
main IDs for each orthologue written to ENSMUSG00000031628.mainIDs.txt
all IDs for each orthologue written to ENSMUSG00000031628.allIDs.txt
sequeences for the main IDs dumped to ENSMUSG00000031628.ort.fa
--
Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000027997
main IDs for each orthologue written to ENSMUSG00000027997.mainIDs.txt
all IDs for each orthologue written to ENSMUSG00000027997.allIDs.txt
sequeences for the main IDs dumped to ENSMUSG00000027997.ort.fa
--
Obtaining 1-to-1 and 1-to-many orthologues for ENSG00000003400
main IDs for each orthologue written to ENSG00000003400.mainIDs.txt
all IDs for each orthologue written to ENSG00000003400.allIDs.txt
sequeences for the main IDs dumped to ENSG00000003400.ort.fa
.
.
(en mi ordenador, para estos 3 IDs tarda unos 5 minutos)


El resultado, para cada identificador ENSG.., son como véis 3 ficheros:
ENSG*.mainIDs.txt --> listado con los protein-IDs del transcript principal de cada ortólogo
ENSG*.allIDs.txt --> listado con los protein-IDs de todas las secuencias que ensembl guarda para cada ortólogo
ENSG*.ort.fa --> secuencias de proteina (formato fasta) correspondientes al transcript principal de cada ortólogo

Aclaración por si alguien se ha liado con esto de transcript_principal / todas_las_secuencias. Como sabéis, para cada gen (y para cada uno de sus ortólogos) hay descritos varios transcripts normalmente. ¿Cuál utiliza ensembl como transcript principal? El más largo. La secuencia de proteína de ese se graba en OUT3, su ID en OUT1 y las IDs de todos en OUT2.


El código de get_orthologues.pl es el siguiente:
#!/usr/bin/perl
# Script to retrieve protein sequences for 1-to-1(and 1-to-many)-orthologues for a given gene.
# (adaptado del que me pasó Bert Overduin de helpdesk@enseml.org)

# Definicion de librerias y modulos

use strict;
use warnings;

use lib '/Users/jramon/Soft/ensembl_API/ensembl/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-compara/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-variation/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-functgenomics/modules/';
use lib '/Users/jramon/Soft/BioPerl-1.6.1/';
use Bio::EnsEMBL::Registry;

my $reg = "Bio::EnsEMBL::Registry";
$reg->load_registry_from_url('mysql://anonymous@ensembldb.ensembl.org');

my $member_adaptor = $reg->get_adaptor("Multi", "compara", "Member");
my $homology_adaptor = $reg->get_adaptor("Multi", "compara", "Homology");


# Leer identificadores de genes y buscar ortologos para cada uno
# guardando IDs en OUT1,2 y secuencias en OUT3

open (INPUT_IDs,"$ARGV[0]");
while () {
my $gene_stable_id = $_;
chop $gene_stable_id;

my $member = $member_adaptor->fetch_by_source_stable_id("ENSEMBLGENE", $gene_stable_id);

open (OUT1,">$gene_stable_id.mainIDs.txt");
open (OUT2,">$gene_stable_id.allIDs.txt");
open (OUT3,">$gene_stable_id.ort.fa");
print "--\nObtaining 1-to-1 and 1-to-many orthologues for $gene_stable_id\n";
print "main IDs for each orthologue written to $gene_stable_id.mainIDs.txt\n";
print "all IDs for each orthologue written to $gene_stable_id.allIDs.txt\n";
print "sequences for the main IDs dumped to $gene_stable_id.ort.fa\n";

my $all_homologies = $homology_adaptor->fetch_all_by_Member($member);

foreach my $homology (@$all_homologies) {
if($homology->description eq (('ortholog_one2one')||('ortholog_one2many'))){
my $members = $homology->get_all_Members();
foreach my $member (@$members) {
if($member->stable_id ne $gene_stable_id){
my $canonical_peptide = $member->get_canonical_peptide_Member;
print OUT1
$canonical_peptide-> stable_id, "\n";

my $gene_adaptor = $reg->get_adaptor($member->genome_db->name, "core", "Gene");
my $gene = $gene_adaptor->fetch_by_stable_id($member->stable_id);
my @transcripts = @{$gene->get_all_Transcripts};
foreach my $transcript(@transcripts){
if($transcript->translation){
print OUT2 $transcript->translation->stable_id, "\n";
}
}

print OUT3
">", $canonical_peptide->stable_id, "\n",
$canonical_peptide->sequence, "\n";
}
}
}
}
close OUT1;
close OUT2;
close OUT3;
}
close INPUT_IDs;


De todas formas, para el volumen de datos que ha de bajar, se me hace un poco lento el código. Se agradecen mejoras en este sentido.

31 de marzo de 2011

Curso 'Next Generation Sequence Analysis: Practice and Departure to New Frontiers'

Hola,
en el 11 GABI Status Meeting, que tuvo lugar recientemente en Alemania,
nos anunciaron el siguiente taller de verano, que tendrá lugar en el mes de Agosto:

'Next Generation Sequence Analysis: Practice and Departure to New Frontiers'
organized by: R. Fries (Animal Breeding, TUM), T. Meitinger (Human Genetics, TUM), K.F.X. Mayer (MIPS, Helmholtz Zentrum München), T. Strom (Human Genetics, Helmholtz Zentrum München)

-> OPEN FOR APPLICATION! (see link below)

Dates:       08.08. – 15.08.2011
Location:  Herrsching am Ammersee, Germany
Description:
The Synbreed Summer School provides an introduction to next generation sequence analysis for PhD students and postdoctoral researchers in animal and plant breeding. The course consists of a practical part with hands-on exercises guided by experienced Synbreed scientists and guest lecturers covering the frontiers of next gen sequencing in animal and plant improvement.
Program Overview:
  • Practice the pipeline from next gen sequence reads to annotated variants (alignment, visualization, variant calling, imputing and annotation)
  • Visiting the sequencing facility of the Helmholtz Zentrum Munich
  • Guest lectures about advances in sequencing technology, genome assembly, genomic selection, next gen population genomics, efficient computation and agricultural genome projects
  • Download preliminary programme
Target Group and Requirements:
  • The course is directed towards PhD students and postdoctoral researchers in animal and plant breeding
  • Knowledge of Linux and a scripting language (e.g. Phyton)
  • Dual-core processor equipped laptop with a recent Linux distribution (e.g. Ubuntu 10.10) and at least 150 GB free storage. Linux setup assistance will be provided during the introductory session.
Course language: English
Costs:
Course fee                                    500 Euro    
Participant in single bedroom*         826 Euro
Participant in double bedroom*        756 Euro
(*includes: full accommodation and meals)
Services at Course Location: please visit http://www.hdbl-herrsching.de.
Application deadline:  01.06.2011
Number of participants: 20

Contact:

Project Coordination Synbreed
Natalie Ohl / Wolf-Christian Saul
Chair of Plant Breeding
Center of Life and Food Sciences Weihenstephan
Technische Universität München
Emil-Ramann-Str. 4
D-85354 Freising
Ph.:    +49 (0)8161/71-5226
Fax:    +49 (0)8161/71-4511
Email: synbreed[at]wzw.tum.de

9 de diciembre de 2010

Compresión de secuencias de ADN

Un problema que me ha surgido recientemente es el de manejar en memoria RAM grandes volúmenes de secuencias de ADN, en mi caso genes, pero que podrían ser también lecturas de secuenciación masiva. Cuando el tamaño total de las secuencias supera la RAM disponible el sistema operativo se verá forzado a hacer swapping/intercambio con el disco y el proceso se ralentiza considerablemente. Una de las posibles soluciones, la que nos ocupa hoy, es la compresión de las secuencias, de manera que ocupen menos espacio físico en RAM.
El siguiente programa muestra cómo podemos realizar esto en Perl, y nos sugiere que esta estrategia vale realmente la pena a medida que el tamaño de las secuencias se aproxima a las 200 bases:




 use strict;  
 use Compress::Zlib;  
   
 # http://perldoc.perl.org/Compress/Zlib.html  
   
 my (@S,$source,$compressed,$uncompressed);  
   
 print "size(source)\tsize(compressed)\tratio\n";  
 for(my $size=25;$size<=500;$size += 25)  
 {  
    # secuencia aleatoria de longitud $size  
    @S = map { ('A','C','G','T')[rand(4)] } 1 .. $size;  
    $source = join(/ /,@S); #print "$source\n";  
   
    $compressed = compress($source);  
    $uncompressed = uncompress($compressed) ;  
      
    printf("%d\t%d\t%g\n",  
       $size,length($compressed),length($compressed)/$size);  
 }  
   

Un saludo,
Bruno

13 de julio de 2010

Eliminar secuencias redundantes y transcritos repetidos (II, con CD-HIT)

Como complemento a la entrada anterior de Álvaro os paso un enlace al programa
CD-HIT, que es una opción muy cómoda y eficiente para eliminar redundancia dentro de un conjunto de secuencias de proteínas o de DNA. El artículo donde se describe la primera versión de CD-HIT se publicó en Bioinformatics en 2006.


Aunque hay una serie de servidores web que pueden ayudarnos en esta tarea, si tienes un volumen realmente grande de secuencias lo más probable es que necesites descargar el código fuente (en C++) y compilarlo con make. Una vez compilado, el programa tiene múltiples opciones, pero su uso es muy sencillo.  El algoritmo se basa en calcular la identidad de todas las secuencias (previamente ordenadas por longitud) por parejas para ir generando clusters de secuencias con una identidad superior al umbral deseado. Por defecto el umbral de identidad es del 90%, medida a lo largo de toda la secuencia (global), no sólo la alineada. Por supuesto se puede modificar este comportamiento y elegir por ejemplo umbrales más bajos o identidades locales.


Terminaré con un ejemplo de uso. Por ejemplo, dado el mismo archivo FASTA de la entrada anterior de Álvaro, para obtener el subconjunto no redundante al 80% sería necesario el siguiente comando en el terminal Linux:

$ cd-hit/cd-hit -i redundante.faa -o noredundante.faa -c 0.8 

Como resultado, además del fichero noredundante.faa, se obtienen otros dos archivos (noredundante.faa.clstr y noredundante.faa.bak.clstr) que contienen
los clusters encontrados y las secuencias redundantes observadas.

5 de julio de 2010

Eliminar secuencias redundantes y transcritos repetidos en genomas y proteomas

Una situación cada vez más habitual en genómica y proteómica es la redundancia de datos. Con las modernas técnicas de secuenciación de alto rendimiento se generan gigas y gigas de datos que es necesario filtrar y corregir para obtener resultados inteligibles.

El problema más habitual y sencillo es la redundancia de secuencias, cuando en un mismo fichero aparecen secuencias idénticas con identificadores diferentes. En otras ocasiones el problema es el contrario, identificadores idénticos apuntan a secuencias parciales de una misma secuencia (especialmente en datos de proteómica o transcriptómica).

El siguiente script filtra un fichero de sequencias en formato FASTA, dos tipos de filtro son posibles: 1) filtrar las secuencias redundantes (repetidas), 2) filtrar las secuencias de transcritos parciales (recuperando únicamente los transcritos de mayor secuencia en ficheros de proteínas)

 #!/usr/bin/perl  
   
 my ($infile) = @ARGV;  
 my $outfile = filter_fasta_file($infile,['longtrans','nonredundant']);  
   
 # Create a FASTA file with custom filtered sequences  
 sub filter_fasta_file {  
   
      my ($infile, $filter, $outfile) = @_;  
   
      if (!defined($outfile)){$outfile=$infile.'.filtered'}  
   
      open(INFILE,"$infile");  
      open(OUTFILE,">$outfile");  
   
      my ($sequences,$names,$headers) = read_FASTA_file($infile);  
   
      my (@previous_sequences,@previous_names,@previous_headers);  
   
      # my (@filtered_sequences,@filtered_names,@fitered_headers);  
   
      for (my $i=0; $i<=$#{$headers}; $i++){  
   
           my $write=1;  
   
           if (in_array($filter,'nonredundant')){  
                my ($found,$posic)=in_array($sequences, $sequences->[$i], 1);  
                foreach my $pos (@{$posic}){  
                     if ($i != $pos && $i > $pos){  
                          $write = 0; last;  
                     }  
                }  
           }  
   
           if (in_array($filter,'longtrans')){  
                $names->[$i] =~ /([^\.]+)\.?(\d?)/;  
                my $name = $1;  
                my $variant = $2;  
                my ($found,$posic)=in_array($names, "/$1/", 1);  
                foreach my $pos (@{$posic}){  
                     if ($i != $pos){  
                          if (length($sequences->[$i]) < length($sequences->[$pos])) {  
                               $write = 0; last;  
                          }  
                          if (length($sequences->[$i]) == length($sequences->[$pos]) && $sequences->[$i] eq $sequences->[$pos]){  
                               $names->[$pos] =~ /([^\.]+)\.?(\d?)/;  
                               my $name2 = $1;  
                               my $variant2 = $2;  
                               if ($variant2 < $variant){  
                                    $write = 0; last;  
                               }  
                          }  
                     }  
                }  
           }  
   
   
           if ($write) {  
                print OUTFILE $headers->[$i]."\n".$sequences->[$i]."\n";  
           }  
   
      }  
   
      close INFILE;  
      close OUTFILE;  
   
      return $outfile;  
 }  
   
   
 #################################################################################  
   
 # Create 3 array references with sequences, names and headers from a FASTA file  
 sub read_FASTA_file  
 {  
      my($FASTAfile,$full_header,$allow_empty_sequences,$allow_multi_rows) = @_;  
   
      if(!$full_header){ $full_header = $allow_empty_sequences = 0 }  
      if(!$allow_empty_sequences){ $allow_empty_sequences = 0 }  
      if(!$allow_multi_rows){ $allow_multi_rows = 0 }  
   
      my (@sequences,@names,@headers,$seq,$seqname,$header);  
      my $headerOK = 0;  
   
      open(FASTA,$FASTAfile) || die "# read_FASTA_file : cannot find $FASTAfile\n";  
      $seq = '';  
      while(<FASTA>){  
           next if(/^#/); # allow comments  
           s/\012\015?|\015\012?|^\s+|\s+$//;  
           if(/>/){  
                $headerOK++;  
                if($seq || ($allow_empty_sequences && $seqname)){  
                     if(!$allow_multi_rows){  
                          $seq =~ s/\d|\s//g;  
                     }  
                     push(@sequences,$seq); # store this FASTA sequence  
                     push(@names,$seqname);  
                     push(@headers,$header);  
                     $seq='';  
                }  
   
                if($full_header){  
                     $seqname = substr($_,1);  
                }else { # get next sequence name  
                     if (/^>\s*([^\t|^\||^;]+)/){  
                          $header = $_;  
                          $seqname = $1;  
                          if (length($seqname)>32 && $seqname =~ /^([^\s]+)+\s/){  
                                    $seqname = $1;  
                          }  
                          $seqname =~ s/^\s+|\s+$// ;  
                     }  
   
                }  
           }else{  
                if(!$allow_multi_rows){   
                     s/\s+//g;  
                }
                $seq .= $_;  
           }  
      }  
      close(FASTA);  
   
      # take care of last FASTA sequence  
      push(@sequences,$seq);  
      push(@names,$seqname);  
      push(@headers,$header);  
   
      if(!$headerOK){ @sequences = @names = @headers = (); } #Aug2007  
   
      return (\@sequences,\@names,\@headers);  
 }  
   
 #################################################################################  
   
 # Searchs a string inside an array  
 sub in_array {  
      my ($array,$search_for,$search_posic) = @_;  
 #      my %items = map {$_ => 1} @$array; # create a hash out of the array values  
 #      my $exist = (exists($items{$search_for}))?1:0;  
 #      return (exists($items{$search_for}))?1:0;  
   
      my @posic;  
   
      # If search a pattern  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @posic = grep { $array->[$_] =~ /$search_for/ } 0 .. $#{$array};  
      } else {  
           @posic = grep { $array->[$_] eq $search_for } 0 .. $#{$array};  
      }  
   
      if (defined($search_posic) && $search_posic){  
           (scalar @posic)?return (1,\@posic):return (0);  
      } else {  
           (scalar @posic)?return 1:return 0;  
      }  
 }  
   

An example FASTA file to test the script:

 >AT1G51370.2 | Symbols: | F-box family protein | chr1:19045615-19046748 FORWARD  
 MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG  
 LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH  
 LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK  
 LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK  
 NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR  
 PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILVMSSFNPS*  
 >AT1G70790.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G50920.1 | Symbols: | GTP-binding protein-related | chr1:18870555-18872570 FORWARD  
 MVQYNFKRITVVPNGKEFVDIILSRTQRQTPTVVHKGYKINRLRQFYMRKVKYTQTNFHA  
 KLSAIIDEFPRLEQIHPFYGDLLHVLYNKDHYKLALGQVNTARNLISKISKDYVKLLKYG  
 DSLYRCKCLKVAALGRMCTVLKRITPSLAYLEQIRQHMARLPSIDPNTRTVLICGYPNVG  
 KSSFMNKVTRADVDVQPYAFTTKSLFVGHTDYKYLRYQVIDTPGILDRPFEDRNIIEMCS  
 ITALAHLRAAVLFFLDISGSCGYTIAQQAALFHSIKSLFMNKPLVIVCNKTDLMPMENIS  
 EEDRKLIEEMKSEAMKTAMGASEEQVLLKMSTLTDEGVMSVKNAACERLLDQRVEAKMKS  
 KKINDHLNRFHVAIPKPRDSIERLPCIPQVVLEAKAKEAAAMEKRKTEKDLEEENGGAGV  
 YSASLKKNYILQHDEWKEDIMPEILDGHNVADFIDPDILQRLAELEREEGIREAGVEEAD  
 MEMDIEKLSDEQLKQLSEIRKKKAILIKNHRLKKTVAQNRSTVPRKFDKDKKYTTKRMGR  
 ELSAMGLDPSSAMDRARSKSRGRKRDRSEDAGNDAMDVDDEQQSNKKQRVRSKSRAMSIS  
 RSQSRPPAHEVVPGEGFKDSTQKLSAIKISNKSHKKRDKNARRGEADRVIPTLRPKHLFS  
 GKRGKGKTDRR*  
 >AT1G70795.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G70790.2 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G36960.1 | Symbols: | unknown protein | chr1:14014796-14015508 FORWARD  
 MTRLLPYKGGDFLGPDFLTFIDLCVQVRGIPLPYLSELTVSFIAGTLGPILEMEFNQDTS  
 TYVAFIRVKIRLVFIDRLRFFRREEAAASNTITDQTHMTSSNSSDISPASPISQPPLPAS  
 LPSHDSYFDAGIQASRLVNPRAISQHHFSSSYSDFKGKEKAKIKIGECSKRKKDKQVDSG  
 T*  
 >AT1G51370.1 | Symbols: | F-box family protein | chr1:19045615-19047141 FORWARD  
 MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG  
 LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH  
 LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK  
 LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK  
 NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR  
 PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILEMVKNQSTRRHGEKREPNVMVS  
 TVPWCLVSSLKFVELKRSIPRYEGEMELVRYVLTNSTVLKKLRLNVYYTKKAKCAFLTEL  
 VAIPRCSSTCVVLVL*