10000bp FASTA 1000bp 10000bp 3' i = 1 remainder = seq.window_search(10000, 9000) do subseq puts subseq.to_fasta("segment #{i}", 60) i += 1 puts remain

Size: px

Start display at page:

Download "10000bp FASTA 1000bp 10000bp 3' i = 1 remainder = seq.window_search(10000, 9000) do subseq puts subseq.to_fasta("segment #{i}", 60) i += 1 puts remain"

おきまさひでやま
6 years ago
Views:

1 BioRuby (Bio::Sequence ) atgcatgcaaaa codontable.rb seq = Bio::Sequence::NA.new("atgcatgcaaaa") puts seq puts seq.complement puts seq.subseq(3,8) p seq.gc_percent p seq.composition puts seq.translate puts seq.translate(2) puts seq.translate(1,11) p seq.translate.codes p seq.translate.names p seq.translate.composition p seq.translate.molecular_weight puts seq.complement.translate # # (Sequence::NA ) # 3 8 # GC (Float) # (Hash) # (Sequence::AA ) # # 11 # (Array) # (Array) # (Hash) # (Float) # Bio::Sequence::NA Bio::Sequence::AA Bio::Sequence Bio::Sequence Ruby String String subseq(from,to) String [] Ruby subseq 0 base 1 base from, to 0 nil window_search(window_size, step_size) Bio::Sequence::NA Bio::Sequence::AA GC% seq.window_search(100) do subseq puts subseq.gc 15 seq.window_search(15, 3) do subseq puts subseq.translate

2 10000bp FASTA 1000bp 10000bp 3' i = 1 remainder = seq.window_search(10000, 9000) do subseq puts subseq.to_fasta("segment #{i}", 60) i += 1 puts remainder.to_fasta("segment #{i}", 60) codon_usage = Hash.new(0) seq.window_search(3, 3) do subseq codon_usage[subseq] += 1 10 seq.window_search(10, 10) do subseq puts subseq.molecular_weight Bio::Sequence::NA input_seq = ARGF.read # my_naseq = Bio::Sequence::NA.new(input_seq) my_aaseq = my_naseq.translate puts my_aaseq na2aa.rb gtggcgatctttccgaaagcgatgactggagcgaagaaccaaagcagtgacatttgtctg atgccgcacgtaggcctgataagacgcggacagcgtcgcatcaggcatcttgtgcaaatg tcggatgcggcgtga my_naseq.txt %./na2aa.rb my_naseq.txt VAIFPKAMTGAKNQSSDICLMPHVGLIRRGQRRIRHLVQMSDAA*

3 % ruby -r bio -e 'p Bio::Sequence::NA.new($<.read).translate' my_naseq.txt GenBank (Bio::GenBank ) GenBank ftp://ftp.ncbi.nih.gov/genbank/.seq gb2fasta ID FASTA gets DELIMITER GenBank GenBank // RS (record separator) while entry = gets(bio::genbank::delimiter) gb = Bio::GenBank.new(entry) # GenBank print ">#{gb.accession} " puts gb.definition puts gb.naseq # ACCESSION # DEFINITION # Sequence::NA Bio::FlatFile ff = Bio::FlatFile.new(Bio::GenBank, ARGF) ff.each_entry do gb definition = "#{gb.accession} #{gb.definition}" puts gb.naseq.to_fasta(definition, 60) FASTA ff = Bio::FlatFile.new(Bio::FastaFormat, ARGF) ff.each_entry do f puts "definition : " + f.definition puts "nalen : " + f.nalen.to_s puts "naseq : " + f.naseq Bio::FlatFile.new BioRuby Bio::DB open

4 ff = Bio::GenBank.open("gbvrl1.seq") ff.each_entry do gb definition = "#{gb.accession} #{gb.definition}" puts gb.naseq.to_fasta(definition, 60) GenBank FEATURES ff = Bio::FlatFile.new(Bio::GenBank, ARGF) # GenBank ff.each_entry do gb # ACCESSION puts "# #{gb.accession} - #{gb.organism}" gb.features.each do feature position = feature.position hash = feature.assoc # FEATURES # # /translation= next unless hash['translation'] # gene_info = [ hash['gene'], hash['product'], hash['note'], hash['function'] ].compact.join(', ') # puts ">NA splicing('#{position}') : #{gene_info}" puts gb.naseq.splicing(position) # puts ">AA translated by splicing('#{position}').translate" puts gb.naseq.splicing(position).translate # /translation= puts ">AA original translation" puts hash['translation'] assoc Feature qualifier qualifier 1 feature splicing GenBank position exon BioRuby splicing GenBank position Bio::Locations position Bio::Locations bio/location.rb GenBank feature position

5 naseq.splicing('join( ,complement( ), ') Locations locs = Bio::Locations.new('join(( ) ,1..855)') naseq.splicing(locs) Bio::Sequence::AA splicing aaseq.splicing(' ') GenBank BioRuby GenBank Bio::FlatFile Bio::FlatFile.new BioRuby (Bio::GenBank Bio::KEGG::GENES ) ff = Bio::FlatFile.new(Bio::, ARGF) FlatFile ff = Bio::FlatFile.auto(ARGF) ff = Bio::FlatFile.auto(ARGF) ff.each_entry do entry p entry.entry_id # ID p entry.definition # p entry.seq # entry_id ID definition reference organism seq naseq aaseq bio/db.rb references Bio::Reference Array reference Bio::Reference

6 FASTA Bio::Fasta FASTA query.pep FASTA ssearch FASTA fasta34 FASTA target.pep FASTA query.pep # FASTA ssearch factory = Bio::Fasta.local('fasta34', ARGV.pop) # FastaFormat ff = Bio::FlatFile.new(Bio::FastaFormat, ARGF) # FastaFormat ff.each do entry # '>' $stderr.puts "Searching... " + entry.definition # FASTA Fasta::Report report = factory.query(entry) # report.each do hit # evalue if hit.evalue < # evalue print "#{hit.query_id} : evalue #{hit.evalue} t#{hit.target_id} at " p hit.lap_at f_search.rb %./f_search.rb query.pep target.pep > f_search.out factory FASTA Fasta query seq = ">test seq nyqvleeigrgsfgsvrkvihiptkkllvrkdikyghmnske" seq.fasta(factory) factory fasta FASTA FASTA ktup ktup 1 10 factory = Bio::Fasta.local('fasta34', 'target.pep', '-b 10') factory.ktup = 1

7 Bio::Fasta#query Bio::Fasta::Report Report FASTA report.each do hit puts hit.evalue # E-value puts hit.sw # Smith-Waterman (*) puts hit.identity # % identity puts hit.overlap # puts hit.query_id # ID puts hit.query_def # puts hit.query_len # puts hit.query_seq # puts hit.target_id # ID puts hit.target_def # puts hit.target_len # puts hit.target_seq # puts hit.query_start # puts hit.query_ # puts hit.target_start # puts hit.target_ # puts hit.lap_at # Bio::Blast::Report FASTA Bio::Fasta::Report fasta report = factory.query(entry) puts factory.output query factory output GenomeNet (fasta.genome.jp) Bio::Fasta.remote Bio::Fasta.local GenomeNet nr-aa, genes, vgenes.pep, swissprot, swissprot-upd, pir, prf, pdbstr nr-nt, genbank-nonst, gbnonst-upd, dbest, dbgss, htgs, dbsts, embl-nonst, embnonst-upd, genes-nt, genome, vgenes.nuc program 'fasta' program 'tfasta' program 'fasta' program = 'fasta' database = 'genes'

8 factory = Bio::Fasta.remote(program, database) factory.query BLAST Bio::Blast BLAST GenomeNet (blast.genome.jp) Bio::Fasta API Bio::Blast f_search.rb # BLAST factory = Bio::Blast.local('blastp', ARGV.pop) GenomeNet Bio::Blast.remote FASTA program program 'blastp' program 'tblastn' program 'blastx' program 'blastn' Bio::Blast -m 8 -m 7 XML Ruby XMLParser REXML BLAST XML XMLParser, REXML Bio::Fasta::Report Bio::Blast::Report Hit BLAST bit_score midline report.each do hit puts hit.bit_score # bit (*) puts hit.query_seq # puts hit.midline # midline (*) puts hit.target_seq # puts hit.evalue puts hit.identity puts hit.overlap puts hit.query_id puts hit.query_def puts hit.query_len puts hit.target_id puts hit.target_def puts hit.target_len puts hit.query_start puts hit.query_ puts hit.target_start puts hit.target_ puts hit.lap_at # E-value # % identity # # ID # # # ID # # # # # # #

9 Hsp Hit Hit Hsp Hsp blastpgp Iteration Bio::Blast::Report Bio::Blast::Report::Iteration Array Bio::Blast::Report::Hits Array Bio::Blast::Report::Hsp Array BLAST bio/appl/blast/*.rb BLAST BLAST Bio::Blast Bio::Blast::Report Bio::Blast.reports blastall -m 7 XML # XML Bio::Blast::Report Bio::Blast.reports(ARGF) do report puts "Hits for " + report.query_def + " against " + report.db report.each do hit print hit.target_id, " t", hit.evalue, " n" if hit.evalue < hits_under_0.001.rb %./hits_under_0.001.rb *.xml BLAST *.xml Blast OS XML XML Blast D -m Blast NCBI BioRuby GenomeNet CGI -m 8 BioRuby blast query Bio::Blast::Report.new exec_ Bio::Blast private factory = Bio::Blast.remote(program, db, option, ' ') BioRuby

10 PubMed (Bio::PubMed ) NCBI PubMed ARGV.each do id entry = Bio::PubMed.query(id) # PubMed medline = Bio::MEDLINE.new(entry) # Bio::MEDLINE reference = medline.reference # Bio::Reference puts reference.bibtex # BibTeX pmfetch.rb %./pmfetch.rb PubMed ID (PMID) NCBI MEDLINE BibTeX # keywords = ARGV.join(' ') # PubMed entries = Bio::PubMed.search(keywords) entries.each do entry medline = Bio::MEDLINE.new(entry) # Bio::MEDLINE reference = medline.reference # Bio::Reference puts reference.bibtex # BibTeX pmsearch.rb %./pmsearch.rb genome bioinformatics PubMed BibTeX NCBI E-Utils esearch, efetch keywords = ARGV.join(' ') options = { 'maxdate' => '2003/05/31', 'retmax' => 1000, }

11 entries = Bio::PubMed.esearch(keywords, options) Bio::PubMed.efetch(entries).each do entry medline = Bio::MEDLINE.new(entry) reference = medline.reference puts reference.bibtex pmsearch.rb E-Utils E-Utils bibtex BibTeX bibitem nature nar Bio::Reference REFERENCE BibTeX BibTeX BibTeX TeX %./pmfetch.rb >> genoinfo.bib %./pmsearch.rb genome bioinformatics >> genoinfo.bib genoinfo.bib documentclass{jarticle} begin{document} bibliographystyle{plain} KEGG ~ cite{pmid: } bibliography{genoinfo} {document} hoge.tex % platex hoge % bibtex hoge # genoinfo.bib % platex hoge # % platex hoge # hoge.dvi bibitem.bib Reference#bibitem pmfetch.rb pmsearch.rb puts reference.bibtex puts reference.bibitem

12 documentclass{jarticle} begin{document} KEGG ~ cite{pmid: } begin{thebibliography}{00} bibitem{pmid: } Kanehisa, M., Goto, S. KEGG: kyoto encyclopedia of genes and genomes., { em Nucleic Acids Res}, 28(1):27--30, {thebibliography} {document} begin{thebibliography} hoge.tex % platex hoge # % platex hoge # BioRuby BioRuby samples/ to be written... OBDA OBDA (Open Bio Database Access) 2002 Arizona Cape Town BioHackathon BioPerl, BioJava, BioPython, BioRuby BioRegistry (Directory) BioFlat 2 BDB BioFetch HTTP BioSQL MySQL PostgreSQL schema <URL: spec CVS cvs.open-bio.org obf-common/ BioRegistry ~/.bioinformatics/seqdatabase.ini /etc/bioinformatics/seqdatabase.ini BioRuby /etc/bioinformatics/

13 ~/.bioinformatics/ open-bio.org seqdatabase.ini bioruby stanza [ ] protocol= location= BioRuby location MySQL protocol index-flat index-berkeleydb biofetch biosql bsane-corba xembl BioRuby index-flat, indexberkleydb, biofetch biosql BioRegistry reg = Bio::Registry.new # serv = reg.get_database('genbank') # ID entry = serv.get_by_id('aa2cg') serv [genbank] protocol Bio::SQL Bio::Fetch nil OBDA get_by_id BioFetch BioSQL BioFlat BioFlat index-flat Berkeley DB (bdb) indexberkeleydb bioruby bioflat % bioflat --makeindex [--format ] BioRuby % bioflat ID

14 GenBank gbbct*.seq % bioflat --makeindex my_bctdb --format GenBank gbbct*.seq % bioflat my_bctdb A16STM262 Ruby bdb Berkeley DB % bioflat --makeindex-bdb [--format ] makeindex bdb BioFetch BioFetch CGI CGI HTTP ID BioRuby BioHackathon GenomeNet DBGET BioFetch bioruby.org BioRuby sample/ BioFetch BioRuby EBI BioFetch BioRuby biofetch % biofetch db_name entry_id 3. Bio::Fetch serv = Bio::Fetch.new(server_url) entry = serv.fetch(db_name, entry_id) 4. BioRegistry Bio::Fetch reg = Bio::Registry.new serv = reg.get_database('genbank') entry = serv.get_by_id('aa2cg') BioRegistry seqdatabase.ini [genbank] protocol=biofetch location= biodbname=genbank BioRegistry Bio::Fetch URL BioFetch Bio::KEGG::GENES, Bio::AAindex1 BioFetch KEGG GENES Halobacterium

15 (VNG1467G) AAindex α (BURA740101) 15 entry = Bio::Fetch.query('hal', 'VNG1467G') aaseq = Bio::KEGG::GENES.new(entry).aaseq entry = Bio::Fetch.query('aax1', 'BURA740101') helix = Bio::AAindex1.new(entry).index position = 1 win_size = 15 aaseq.window_search(win_size) do subseq score = subseq.total(helix) puts [ position, score ].join(" t") position += 1 Bio::Fetch.query bioruby.org BioFetch KEGG/GENES hal AAindex aax1 BioFetch query BioSQL to be written... KEGG API KEGG_API.rd.ja <URL: APPENDIX to be written...

BioRuby Ruby Bioinfomatics Blast BioPerl, BioJava, BIoPython Ruby

BioRuby, BioRuby Ruby Bioinfomatics Blast BioPerl, BioJava, BIoPython Ruby Open Bio* O B F -- Open Bio Foundation BioRuby Ensembl BioCaml BioPerl OmniGene BioLisp BioPython GMOD BioConductor BioJava Apollo