CMSC423: Bioinformatic Algorithms, Databases and ToolsLecture 4Writing bioinformatics softwareBiological databasesCMSC423 Fall 2008 2Writing bioinformatics softwareCMSC423 Fall 2008 3Libraries/utilities• Bio::Perl (Perl)• BioJava (Java)•BioPython (Python)•BioRuby (Ruby)• seqAn (C++)•Bioconductor (R)•Chado (SQL)CMSC423 Fall 2008 4Bio::Perl• http://www.bioperl.orguse Bio::Perl;my $seq = read_sequence(“mytest.fa”, “fasta”);my $gbseq = read_sequence(“mytest.gb”, “genbank”);write_sequence(“>test.fasta”, 'fasta', $gbseq);' vs “ ?CMSC423 Fall 2008 5Bio::Perl• Homework question #5use Bio:Perl;while ($seq = read_sequence(“test.fa”, 'fasta')) {if ($seq ->length() > 500) {print $seq->primary_id(), “\n”;}}Note: you still need to write your own version...CMSC423 Fall 2008 6Bio::Perl• Other useful stuff$seqio = new Bio::SeqIO(-format => 'largefasta', -file => 't/data/genomic-seq.fasta');$pseq = $seqio->next_seq();$gb = new Bio::DB::GenBank;$seq1 = $gb->get_Seq_by_id('MUSIGHBA1');etc...CMSC423 Fall 2008 7BioJava• http://www.biojava.orgimport org.biojava.bio.*;String filename = args[0]; BufferedInputStream is = new BufferedInputStream(new FileInputStream(filename)); //get the appropriate Alphabet Alphabet alpha = AlphabetManager.alphabetForName(args[1]); //get a SequenceDB of all sequences in the file SequenceDB db = SeqIOTools.readFasta(is, alpha);CMSC423 Fall 2008 8BioJava• Question 5BufferedReader br = new BufferedReader(new FileReader(args[0])); String format = args[1]; String alphabet = args[2]; SequenceIterator iter = (SequenceIterator)SeqIOTools.fileToBiojava(format,alphabet, br);while (iter.hasNext()){Sequence seq = iter.nextSequence();if (seq.length() > 500) {System.out.println(seq.getName());}}CMSC423 Fall 2008 9BioJava...more• Same as Bio::Perl:– can directly connect to databases–various sequence manipulations (reverse complement, translate, etc.)–basic bioinformatics algorithms– etc.CMSC423 Fall 2008 10BioPython• http://www.biopython.orgfrom Bio import SeqIOhandle = open(“file.fasta”)seq_record = SeqIO.parse(handle, “fasta”)SeqIO.write(my_records, handle2, "fasta")CMSC423 Fall 2008 11BioPython• Question 5from Bio import SeqIOhandle = open("test.fasta")for seq_record in SeqIO.parse(handle, "fasta") :if len(seq_record) > 500 :print seq_record.idhandle.close()CMSC423 Fall 2008 12BioPython...more• Same as Bio::Perl:– can directly connect to databases–various sequence manipulations (reverse complement, translate, etc.)–basic bioinformatics algorithms– etc.CMSC423 Fall 2008 13BioRuby• http://www.bioruby.orgrequire 'bio'input_seq = ARGF.read # reads all files inargumentsmy_naseq = Bio::Sequence::NA.new(input_seq)CMSC423 Fall 2008 14BioRuby• Question 5#!/usr/bin/env rubyrequire 'bio'ff = Bio::FlatFile.new(Bio::FastaFormat, ARGF)ff.each_entry do |f| if f.length > 500 puts f.entry_id endendCMSC423 Fall 2008 15BioRuby...more• Same as Bio::Perl:– can directly connect to databases–various sequence manipulations (reverse complement, translate, etc.)–basic bioinformatics algorithms– etc.CMSC423 Fall 2008 16SeqAn• http://www.seqan.de#include <seqan/sequence.h>#include <seqan/file.h>using namespace seqan;using namespace std;String <Dna> seq;String<char> name;fstream f;f.open(“test.fasta”);readMeta(f, name, Fasta());readMeta(f, seq, Fasta());CMSC423 Fall 2008 17SeqAn• Question 5String <Dna> seq;String<char> name;fstream f;f.open(“test.fasta”);while (! f.eof()){readMeta(f, name, Fasta());readMeta(f, seq, Fasta());if (length(seq)){cout << name << endl;}}CMSC423 Fall 2008 18SeqAn...more• Not quite as much as Perl/Java/Python, but still lots of utilities (including graph algorithms)CMSC423 Fall 2008 19R/BioConductor• http://www.bioconductor.org• Mainly for statistical applications, e.g. microarray analysislibrary("affy")library("geneplotter")library("gplots")data <- ReadAffy()eset <- rma(data)e <- exprs(eset)heatmap.2(e, margin=c(15,15), trace="none", col=redgreen(25), cexRow=0.5)CMSC423 Fall 2008 20R/BioConductor• Book has lots of examples• Worth learning more about it – easy to do various cool things•example... if timeCMSC423 Fall 2008 21Chado• http://www.gmod.org• Relational schema for storing biological data types in a relational database (e.g. MySQL, Oracle, Sybase, ...)SELECT o.organism_id,o.abbreviation,o.genus,o.species, o.common_name, count(f.feature_id) as n_features, o.comment FROM organism o LEFT JOIN feature f USING (organism_id) GROUP by o.organism_id,o.abbreviation,o.genus,o.species, o.common_name,o.comment ORDER BY o.genus,o.speciesCMSC423 Fall 2008 22CMSC423 Fall 2008 23Chado...more• Bio... generally provide ability to interface with relational database.• Understanding SQL and Chado is useful irrespective of language used.•Relational DB particularly useful for web services• Gbrowse example....if timeCMSC423 Fall 2008 24Programming for bioinformatics• Details of specialized libraries beyond scope of course• Good software engineering practices are essential•Often, “correct” is undefined – output of program defines correctness•Pitfalls – e.g. papers retracted from Science due to software bugs•Key – be proactive and learn by yourselves from online
View Full Document
Unlocking...