230 likes | 335 Views
Introducción a Bioperl. Verónica Jiménez Jacinto vjimenez @ ibt.unam.mx Enero 2010. ¿Qué es Bioperl?. Bioperl es un esfuerzo comunitario para producir código en Perl, el cual sea útil, produciendo aplicaciones para la bioinformática, genómica y ciencias biológicas en general.
E N D
Introducción a Bioperl Verónica Jiménez Jacinto vjimenez @ ibt.unam.mx Enero 2010
¿Qué es Bioperl? • Bioperl es un esfuerzo comunitario para producir código en Perl, el cual sea útil, produciendo aplicaciones para la bioinformática, genómica y ciencias biológicas en general. • Es un proyecto de software libre y fue fundamental en el proyecto de secuenciación del genoma humano. • Consiste en un conjunto de módulos que facilita el desarrollo en Perl de herramientas bioinformáticas.
Clase Atributos Metodos Objeto Objeto Objeto ¿Por qué es atractivo BioPerl? • Paradigma de programación orientada a objetos.
Instalación • Por medio de la aplicación 'Sistema -> Administración -> gestor Synaptic', en realidad la manera sencilla de instalar cualquier software, Lo mismo se puede lograr desde la Terminal, por medio de: $ sudo apt-get install 'nombre:software' . • Usando BIOPERL BUNDLE Instala Bioperl usando CPAN.: >perl -MCPAN -e "install Bundle::BioPerl" Otra manera; >perl -MCPAN -e shell cpan >install Bundle::BioPerl • Instalando BIOPERL usando el shell de CPAN : >perl -MCPAN -e shell Then find the name of the Bioperl version you want: cpan>d /bioperl/ Ahora instala: cpan>install B/BI/BIRNEY/bioperl-1.4.tar.gz • Instalando BIOPERL usando 'make' Descarga, descomprime y desempaqueta el archivo: >gunzip bioperl-1.2.tar.gz >tar xvf bioperl-1.2.tar >cd bioperl-1.2 Luego usa el comando make: >perl Makefile.PL >make >make test • Para windows, se puede descargar una version en : http://www.activestate.com/activeperl/
¿Qué se puede hacer con Bioperl? • Acceder a secuencias locales o remotas • Transformar formatos de diferentes Bases de datos • Manipular secuencias individuales • Buscar secuencias similares • Crear y manipular alineamientos de secuencias • Buscar genes y estructuras geonómicas sobre DNA • Desarrollar código para leer anotaciones
Que se puede hacer con bioperl? Sequences • Bio::Seq es el principal objeto de la clase secuencia de Bioperl. • Bio::PrimarySeq es un objeto secuencia sin características • Bio::SeqIO Proporciona funciones para leer y escribir en secuencias desde archivos • Bio::Tools::SeqStats proporciona estadísticas sobre secuencias. • Bio::LiveSeq::* maneja cambio de secuencias. • Bio::Seq::LargeSeq proporciona soporte para manejar secuencias muy grandes.
Seq es el objeto cetral para manipular secuencias. Es una secuencia con características. use Bio::Perl; # this script will only work with an internet connection # on the computer it is run on $seq_object = get_sequence('swissprot',"ROA1_HUMAN"); my $seq_object2 = get_sequence('embl',"AI129902"); my $seq_object3 = get_sequence('genbank',"AI129902"); write_sequence(">roa1.fasta",'fasta',$seq_object2); cat roa1.fasta >unknown id qc41b07.x1 Soares_pregnant_uterus_NbHPU Homo sapiens cDNA clone IMAGE:1712149 3' similar to SW:ROA1_ SCHAM P21522 HETEROGENEOUS NUCLEAR RIBONUCLEOPROTEIN A1, A2/B1 HOMOLOG. ;contains MSR1.b2 MSR1 repetitive elemen t ;, mRNA sequence. CTCCGCGCCAACTCCCCCCACCCCCCCCCCACACCCC get_secuencias.pl
use Bio::Perl; # this script will only work with an # Internet connection # on the computer it is run on $seq_object = get_sequence('swissprot',"ROA1_HUMAN"); #uses the default database -nr in this case $blast_result = blast_sequence($seq_object); write_blast(">roa1.blast",$blast_result);
# gets a sequence from a file $seqio = Bio::SeqIO->new( '-format' => 'embl' , -file => 'myfile.dat'); $seqobj = $seqio->next_seq(); # get from database $db = Bio::DB::swiss->new(); $seqobj = $db->get_Seq_by_acc('ROA1_HUMAN'); # make from strings in script $seqobj = Bio::Seq->new( -display_id => 'my_id', -seq => $sequence_as_string);
$seqobj = new_sequence("ATTGGTTTGGGGACCCAATTTGTGTGTTATATGTA", “un nombre", "AL12232"); $seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj); $monomer_ref =$seq_stats->count_monomers(); $codon_ref = $seq_stats->count_codons(); $weight = $seq_stats->get_mol_wt($seqobj);
# gets sequence as a string from sequence object $seqstr = $seqobj->seq(); # actual sequence as a string $seqstr = $seqobj->subseq(10,50); # slice in biological coordinates # retrieves information from the sequence # features must implement Bio::SeqFeatureI interface @features = $seqobj->get_SeqFeatures(); # just top level foreach my $feat ( @features ) { print "Feature ",$feat->primary_tag," starts ",$feat->start," ends ", $feat->end," strand ",$feat->strand,"\n"; # features retain link to underlying sequence object print "Feature sequence is ",$feat->seq->seq(),"\n" } # sequences may have a species if( defined $seq->species ) { print "Sequence is from ",$species->binomial_name," [",$species->common_name,"]\n"; } # annotation objects are Bio::AnnotationCollectionI's $ann = $seqobj->annotation(); # annotation object # references is one type of annotations to get. Also get # comment and dblink. Look at Bio::AnnotationCollection for # more information foreach my $ref ( $ann->get_Annotations('reference') ) { print "Reference ",$ref->title,"\n"; } # you can get truncations, translations and reverse complements, these # all give back Bio::Seq objects themselves, though currently with no # features transfered my $trunc = $seqobj->trunc(100,200); my $rev = $seqobj->revcom(); # there are many options to translate - check out the docs my $trans = $seqobj->translate(); # these functions can be chained together my $trans_trunc_rev = $seqobj->trunc(100,200)->revcom->translate();
¿Qué se puede hacer en bioperl? Databases • Bio::DB::GenBank proporciona acceso a GenBank • Bio::Tools::Run::StandAloneBlast corre BLAST localmente. • Bio::Tools::Run::RemoteBlast corre BLAST remotamente. • Bio::Tools::BPlite parsea un reporte BLAST • Bio::Tools::BPpsilite parsea un reporte psiblast • Bio::Tools::HMMER::Results parsea resultados de Cadenas de Markov.
Alignments • Bio::SimpleAlign manipula y despliega alineamientos de múltiples secuencias • Bio::LocatableSeq Son objetos secuencias con puntos de inicio y final para su localización relativa a otras secuencias o alineamientos. • Bio::Tools::pSW alinea dos secuencias con el algoritmo Smith-Waterman. • Bio::AlignIO Alinea dos secuencias con el algoritmo blast • Bio::Clustalw es una interface del paquete Clustalw. • Bio::TCoffee es una interface del paquete Tcoffee • Bio::Variation::Allele maneja conjuntos de allelos. • Bio::Variation::SeqDiff maneja conjuntos de mutaciones y variantes
Features and genes on sequences • Bio::SeqFeature es un objeto con las caracteristicas de una secuencia en Bioperl. • Bio::Tools::RestrictionEnzyme localiza sitios de restriccion sitios en secuencias • Bio::Tools::Sigcleave Encuentra sitios de corte en aminoacidos. • Bio::Tools::OddCodes Rescribe secuencias de aminoacidos con codigos abreviados para especificar analisis estadisticos. (e.g., a hydrophobic/hydrophilic two-letter alphabet). • Bio::Tools::SeqPattern Proporciona soporte para encontrar secuencias de patrones. • Bio::LocationI proporciona una interface para localizar información de una seucencia • Bio::Location::Simple maneja información de la localización de una secuencia, como una simple localización y como un rango. . • Bio::Location::Fuzzy proporciona información de la localización que puede ser inexacta. • Bio::Tools::Genscan es una interface para encontrar genes con el progrma Genscan • Bio::Tools::Sim4::Results (and Exon) es una interface para encontrar exones de genes con el programa Sim4
>more $RHIZO_PUB/RE1PF/CP000138.gbk LOCUS CP000138 642517 bp DNA circular BCT 10-MAR-2006 DEFINITION Rhizobium etli CFN 42 plasmid p42f, complete sequence. ACCESSION CP000138 VERSION CP000138.1 GI:86284836 KEYWORDS . SOURCE Rhizobium etli CFN 42 REFERENCE 1 (bases 1 to 642517) AUTHORS Gonzalez,V., Santamaria,R.I., Bustos,P., Hernandez-Gonzalez,I., Medrano-Soto,A., Moreno-Hagelsieb,G., Janga,S.C., Ramirez,M.A., Jimenez-Jacinto,V., Collado-Vides,J. and Davila,G. TITLE The partitioned Rhizobium etli genome: Genetic and metabolic redundancy in seven interacting replicons … • FEATURES Location/Qualifiers • source 1..642517 • /organism="Rhizobium etli CFN 42" • /mol_type="genomic DNA" • /strain="CFN 42" • /db_xref="taxon:347834" • /plasmid="p42f" • promoter 516..546 • /note="sigma54 panCp promoter; Putative transcription • initiation." • gene 709..1602 • /gene="panC" • /locus_tag="RHE_PF00001" • CDS 709..1602 • /gene="panC" • /locus_tag="RHE_PF00001" • /EC_number="6.3.2.1" • /product="pantoate beta alanine ligase protein"
vjimenez> more /home/genomas/pub/RE1PF/RE1PF_gene_from_GK3.dat • LocusTag GI gene_name product_name position strain gbaccession crossrefs • RHE_PF00001 GI:86284837 panC pantoate beta alanine ligase protein 709..1602 forward CP000138 • CDD:COG0414,CDD:PF02569.4,GI:86284837,InterPro:IPR003721 • RHE_PF00002 GI:86284838 panB ketopantoate hydroximethyltransferase protein 1599..2420 forward • CP000138 CDD:COG0413,CDD:PF02548.4,GI:86284838,InterPro:IPR003700 • RHE_PF00003 GI:86284839 oxyR hydrogen peroxide sensing transcriptional regulator protein, LysR family • complement(2500..3414) reverse CP000138 CDD:COG0583,CDD:PF00126.10,CDD:PF03466.5,GI:86284839,Int • erPro:IPR000847,InterPro:IPR005119 • RHE_PF00004 GI:86284840 katG catalase protein 3559..5745 forward CP000138 CDD:COG0 • 376,CDD:PF00141.9,GI:86284840,InterPro:IPR000763,InterPro:IPR002016
use Bio::Seq; • use Bio::Seq::RichSeq; • use Bio::SeqIO; • use Bio::SeqIO::genbank; • $pathIN=$ARGV[0]; • $filetoRead=$ARGV[1]; #File to Read • $filetoStore=$ARGV[2]; • print "Archivo $filetoRead\n"; • open(OUT,">$filetoStore")|| die "Cannot open output file..$filetoStore\n"; • $in = Bio::SeqIO->new(-file => "$pathIN/$filetoRead", '-format' => 'genbank'); • while ((my $seq = $in->next_seq())){ • $gbaccession = $seq->accession(); foreach my $f ($seq->get_SeqFeatures) { **… • } • } • close(OUT); • print "# finished processing: n_of_seqs=$n_of_seqs\n"; • translate_GBK_to_FileTabs4.pl
if($f->primary_tag() =~ /CDS/){ • $posleft=$f->location->{"_start"}; • $posrigth=$f->location->{"_end"}; • if($f->location->{"_strand"} == 1){ • $strain= "forward"; • }else{ • $strain= "reverse"; } • if($f->has_tag('db_xref')){ • $crossrefs = join(',',sort $f->each_tag_value('db_xref')); • } • if($f->has_tag('locus_tag')){ • $id = join(',',sort $f->each_tag_value('locus_tag')); • } • if($f->has_tag('gene')){ • $gene = join(',',sort $f->each_tag_value('gene')); • } • if($f->has_tag('product')){ • $product = join(',',sort $f->each_tag_value('product')); • } • $header = $id."\t".$gi."\t".$gene."\t".$product."\t".$genepos."\t".$strain."\t".$gbaccession."\t$crossrefs\t"; print OUT "$header\n"; • $n_of_seqs++; }
¿y para Illumina? • Generar una archivo fastq http://www.eead.csic.es/compbio/material/bioperl/node24.html • http://www.lcg.unam.mx/~compu2/cei/
Problemas con Bioperl… • La documentación de Bioperl esta incompleta • Bioperl es grande (mas de 500 módulos) escrito por muchos voluntarios
Referencias: • Curso: Perl en Bioinformática. Autor: Bruno Contreras. url: http://www.eead.csic.es/compbio/material/bioinfoPerl/ • http://www.bioperl.org/ • http://www.bioperl.org/Core/Latest/bptutorial.html • http://www.pasteur.fr/recherche/unites/sis/formation/bioperl/index.html • http://www.bioperl.org/wiki/News