Algunas aplicaciones en Bioinformática

Los ejemplos previos pueden ser de utilidad en diferentes aplicaciones; los siguientes se orientan ya a tareas habituales en la Bioinformática, y algunas se inspiran en cádigo del repositorio Scriptome:

1. cuenta las secuencias de un archivo FASTA (equivalente a grep -c "^>" archivo.fasta | wc -l :

   perl -lne 'if(/^>/){ $total++ } END{ print "secuencias=$total"}'

2. lee un archivo FASTA y formatea las secuencias en una sola línea (ademas las ordena alfabéticamente por sus descripciones):

   perl -lne 'if(/^(>.*)/){ $head=$1 } else { $fa{$head} .= $_ } END{ foreach $s (sort(keys(%fa))){ print "$s\n$fa{$s}\n" }}'

3. elimina secuencias idénticas en un archivo FASTA:

   perl -lne 'if(/^(>.*)/){ $head=$1 } else { $fa{$head} .= $_ } END{ foreach $s (keys(%fa)){ print "$s\n$fa{$s}\n" if(!$uniq{$fa{$s}}); $uniq{$fa{$s}}=1 }}'

4. elimina secuencias de cierta longitud de un archivo FASTA:

   perl -lne 'if(/^(>.*)/){ $head=$1 } else { $fa{$head} .= $_ } END{ foreach $s (keys(%fa)){ print "$s\n$fa{$s}\n" if(length($fa{$s})>100) }}'

5. descarga un archivo GenBank, con ayuda de BioPerl (ver sección 4.2 y ejemplo en el taller de (bio)perl):

   $ perl -MBio::Perl -e '$gb=get_sequence("genbank","CP000524"); write_sequence(">CP000524.gb","genbank",$gb)'

6. descarga una secuencia de UniProt en formato FASTA por servicios web (recurre a módulo SOAP::Lite):

   $ perl -MSOAP::Lite -e '$soap=SOAP::Lite->service("http://www.ebi.ac.uk/Tools/webservices/wsdl/WSDbfetch.wsdl");print $soap->fetchData("uniprot:Q06546","fasta","raw")'

7. convierte un archivo de formato GenBank a FASTA (de nuevo usa BioPerl):

   $ perl -MBio::SeqIO -e '$in=Bio::SeqIO->newFh(-file=>$ARGV[0],-format=>"genbank"); $out=Bio::SeqIO->newFh(-format=>"fasta"); while(<$in>){ print $out $_ }'

8. calcula el N50 de un ensamblaje de secuencias:

   $ perl -lne 'if(/^>(\S+)/){$h=$1} else {$l=length($_);$TL+=$l;$L{$h}+=$l} \
      END{ foreach $s (sort {$L{$b}<=>$L{$a}} keys(%L)){ $t+=$L{$s}; \
      if($t>$TL/2){ print "N50=$L{$s}";exit }}}' ensamblaje.fasta
   

9. crea un histograma (archivo llamado 'hist.pdf') con ayuda de R:

   $ perl -lane '$F[3]=~s/,/./; print $F[3]' all.blast | Rscript -e 'data=scan(file="stdin"); hist(data)'

10. convierte las calidades de un fichero FASTQ a codificación Sanger/Illumina 1.8+:

    $ perl -lne 'tr/\x40-\xff\x00-\x3f/\x21-\xe0\x21/ if($.%4==0); print' reads.fq > reads.sanger.fq