tag:blogger.com,1999:blog-89277782024-03-13T16:18:00.144-04:00The Plant InformaticsRandom musings on the morphology and genomes of flowering plantsAnonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.comBlogger53125tag:blogger.com,1999:blog-8927778.post-62618986376947458362011-10-16T14:52:00.001-04:002011-10-16T14:52:37.315-04:00Rod Page’s VIZBI 2011 annotated links<p> <p><a href="http://iphylo.blogspot.com/">iPhylo</a> is hands-down my favorite blog, whose author is Rod Page. In the blog post “<a href="http://iphylo.blogspot.com/2011/03/some-vizbi-2011-links.html">Some VIZBI 2011 links</a>”, Rod grabbed a list of visualization software that he noted. I stole the whole list and post it here, but with my own annotations. There are some good ones there that I would like to use in the future, so this post is merely a reminder and click-saver. <ul> <li><a href="http://arena3d.org/">Arena 3D</a> - network in 3D <li><a href="http://www.bioblender.net/">BioBlender</a> - cool, artistic protein structure rendering <li><a href="http://www.mquter.qut.edu.au/bio/blastatlas_default.aspx">BLAST Atlas</a> - COG results in 3D <li><a href="http://www.pathogenomics.ca/cerebral/">Cerebral v.2.0</a> - cytoscape plugin <li><a href="http://www.research.ibm.com/visual/projects/chromogram.html">Chromograms</a> - visualize text editing <li><a href="http://clovr.org/">CloVR</a> - microbial / metagenomics seq analysis using cloud computing <li><a href="http://www.cytoscape.org/">Cytoscape</a> - need-i-say-more <li><a href="http://www.emouseatlas.org/emage/home.php">EMAGE</a> - more of a database than standalone tool <li><a href="http://epmv.scripps.edu/">embedded Python Molecular Viewer (ePMV)</a> - protein structure <li><a href="http://www.bewitched.com/fleshmap.html">Fleshmap</a> - desire and body <li><a href="http://genomeview.org/">GenomeView</a> - like Tablet, nextgen genome viewer <li><a href="http://www.research.ibm.com/visual/projects/history_flow/">History Flow</a> - track history of edits, see also Chromograms above <li><a href="http://mkweb.bcgsc.ca/linnet/">HivePlots</a> - novel visualization method called hive plots, didn't like it though <li><a href="http://www.sci.utah.edu/cibc/software/41-imagevis3d.html">ImageVis3D</a> - image viewer <li><a href="http://www.profilegrid.org/">JProfileGrid</a> - edit MSA <li><a href="http://www.many-eyes.com/">Many Eyes</a> - text cloud, maps, etc. <li><a href="http://www.molecularmovies.com/toolkit/">Molecular Maya</a> <li><a href="http://www.molecularmovies.com/">Molecular Movies</a> <li><a href="http://www.pathblast.org/">PathBLAST</a> - query PPI networks <li><a href="http://portnoy.iplantcollaborative.org/">Phyloviewer</a> - from iPlant, biggish tree viewer <li><a href="http://www.virtualpathology.leeds.ac.uk/research/HCI/Powerwall/virtual_reality_powerwall.php">Powerwall</a> - tumor slides <li><a href="http://dev.gramene.org/mockups/dna.html">Quartz Composition of DNA</a> - on MAC only <li><a href="http://chmille4.github.com/Scribl/">Scribl</a> - lightweight, nice module for building genome browsers <li><a href="http://pages.cs.wisc.edu/~dalbers/">Sequence Surveyor</a> - show sequence compositions along phylogenetic tree <li><a href="http://www.bewitched.com/song.html">Shape of Song</a> - not sure why it is in this list <li><a href="http://sybil.sourceforge.net/">Sybil</a> - from TIGR/JCVI, comparative genomics tool <li><a href="http://topiaryexplorer.sourceforge.net/">Topiary Explorer</a> - useful in microbial ecology studies <li><a href="http://caltech.wormbase.org/virtualworm/">Vrtual Worm</a> - as the name suggests <li><a href="http://hint.fm/seer/">Web Seer</a> - visualize google search results <li><a href="http://wholebraincatalog.org/">Whole Brain Catalog</a></li></ul> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-92138403781646071652010-11-09T23:23:00.001-05:002010-11-09T23:23:54.153-05:00Why 10,000 heads from 20,000 tosses suggests fake data?<p>I have this statistical question in mind for some time now. The reason I am interested in this question is that when we some data in a biology paper, sometimes the data are just “too good to be true”. Well, this by itself is not a sufficient reason to discredit a sound paper, but I was always wondering how we can quantify the confidence when seeing such data. </p> <p>I formulated my question in terms of coin tosses and asked on <a href="http://stats.stackexchange.com/questions/4360/statistical-argument-for-why-10-000-heads-from-20-000-tosses-suggests-invalid-dat">CrossValidated</a>. Let me just repost it here. </p> <blockquote> <p>Let's say we are repeatedly tossing a fair coin, and we know number of heads and tails should be roughly equal. When we see a result like 10 heads and 10 tails for a total of 20 tosses, we believe the results and are inclined to believe the coin is fair. <p>Well when you see a result like 10000 heads and 10000 tails for a total of 20000 tosses, I actually would question the validity of the result (did the experimenter fake the data), as I know this is more unlikely than, say a result of 10093 heads and 9907 tails. <p>What is the statistical argument behind my intuition?</p></blockquote> <p>I got my answer back very quickly and someone suggested I had implicitly invoked <a href="http://en.wikipedia.org/wiki/Bayes'_theorem">Bayes theorem</a>. It was nicely explained. After that answer, I did some thinking and calculation and was able to calculate the probability using an explicit model.</p> <p>From prior experience, I know the data faker would forge a data that is very close to the “expected”, i.e. small variance. I will just model it as a <font face="Courier New">Normal(mean=10000, sd=10)</font>. In the mean time, I know the actual coin tosses would follow a binomial distribution <font face="Courier New">Binom(size=20000, p=0.5)</font>. </p> <p><a href="http://lh5.ggpht.com/_srvRoIok9Xs/TNoeVYKZ3mI/AAAAAAAABJs/jlY2XXp_XKU/s1600-h/two_distributions%5B5%5D.png"><img style="background-image: none; border-right-width: 0px; padding-left: 0px; padding-right: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px; padding-top: 0px" title="two_distributions" border="0" alt="two_distributions" src="http://lh4.ggpht.com/_srvRoIok9Xs/TNoeVuzXsdI/AAAAAAAABJw/WaYAjs5Owv4/two_distributions_thumb%5B3%5D.png?imgmax=800" width="510" height="550"></a></p> <p>Note that the fake data has a much sharper peak (low variance) around 10000 – well that’s what fakers do. Let’s say now I have a prior belief of P(fake)=0.1, i.e. before seeing the data, I think the data is forged 10% of the time. </p> <p>I now have everything I need for the Bayes calculation. The formula is:</p> <p><a href="http://lh6.ggpht.com/_srvRoIok9Xs/TNoeV44QBDI/AAAAAAAABJ0/NfRDWCYtUzo/s1600-h/bayes_formula%5B3%5D.png"><img style="background-image: none; border-right-width: 0px; padding-left: 0px; padding-right: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px; padding-top: 0px" title="bayes_formula" border="0" alt="bayes_formula" src="http://lh6.ggpht.com/_srvRoIok9Xs/TNoeWAVDHvI/AAAAAAAABJ4/-3qEgCjQ82Q/bayes_formula_thumb%5B1%5D.png?imgmax=800" width="531" height="83"></a></p> <p>I left all the calculations to R, in <font face="Courier New">calc_prob()</font> function. </p> <p><font face="Courier New">> calc_prob(10000)<br>[1] 0.4399905<br>> calc_prob(10093)<br>[1] 3.088808e-19<br></font></p> <p>Seeing 10000 heads, the posterior probability of the data being rigged is 0.43 (recall my prior belief is only 0.1). In contrast, seeing 10093 heads, the posterior probability is effectively zero. <strong>After all, who would have rigged a data that’s not pretty!</strong></p> <p>R source code is attached here. </p><script src="https://gist.github.com/670348.js?file=cointoss.R"></script> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-66471061018983242222010-08-01T18:32:00.001-04:002010-08-01T18:32:17.886-04:00Gambler’s ruin<p>There is a common casino strategy when playing <a href="http://en.wikipedia.org/wiki/Roulette">roulette</a>: <blockquote> <p><em>Bet on either odd or even, and stick to your choice. Bet $20 at start, if you lose, double the bet; if you win, start with the $20 bet again. Don’t worry about losing, just keep going.</em></p></blockquote> <p>Let’s assume you have $10000. Well, if you adopt this strategy, you’ll only go bankrupt ~1/500 of the time – this happens when you lose 9 times in a row (2**9=512), then you’ll have no money at all. But if you win at any point, you go back to where you were, plus $20 more. For example, let’s say you lose 3 times followed by a win: that’s –20-40-80+160= $20.</p> <p>This simulation (<a href="http://gist.github.com/503788">python script</a>) suggests that if you follow this strategy for a few rounds (say 50), more than 90% of the people will make money, with few people bankrupt.</p> <div style="border-bottom: silver 1px solid; text-align: left; border-left: silver 1px solid; padding-bottom: 4px; line-height: 12pt; background-color: #f4f4f4; margin: 20px 0px 10px; padding-left: 4px; width: 97.5%; padding-right: 4px; font-family: 'Courier New', courier, monospace; direction: ltr; height: 99px; max-height: 200px; font-size: 8pt; overflow: auto; border-top: silver 1px solid; cursor: text; border-right: silver 1px solid; padding-top: 4px" id="codeSnippetWrapper"> <div style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px" id="codeSnippet"><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum1"> 1:</span> $ python casino_strategy.py</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum2"> 2:</span> Simulating 50 rolls per trial and generate 10000 trials</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum3"> 3:</span> mean: 10006, median: 10480</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum4"> 4:</span> 93% of the samples earned money</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum5"> 5:</span> 3% of the samples got bankrupt </pre><!--CRLF--></div></div><br /><p>The following graph simulates 20 trials, and each line tracks their balance at any given round. In the end (round 50), only 1 guy lost money, while the others all earn money, and no people went bankrupt. Notice the typical pattern is <strong>exponential drop followed by rebound in a single round</strong>. Most people earn little, few guys lose a lot. Now the question is, do you want to be that poor guy who lost a lot? <br /><p><a href="http://lh5.ggpht.com/_srvRoIok9Xs/TFX174kcvBI/AAAAAAAABBM/cWDJbYIKQ10/s1600-h/casino_strategy%5B3%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="casino_strategy" border="0" alt="casino_strategy" src="http://lh3.ggpht.com/_srvRoIok9Xs/TFX18Hvp8VI/AAAAAAAABBQ/uFz-oAGAGzM/casino_strategy_thumb%5B1%5D.png?imgmax=800" width="487" height="374"></a> <br /><p>In addition, if you adopt this strategy long enough, eventually you'll go bankrupt, that’s called <a href="http://en.wikipedia.org/wiki/Gambler's_ruin">Gambler's ruin</a>, since “bankrupt” is an absorption state. <br /> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-20034056565238937872010-07-30T23:45:00.001-04:002010-07-30T23:45:55.634-04:00Promethease report: Hyungyong Kim<p>What can you do with your genome deciphered by <a href="https://www.23andme.com/">23andme</a>? For those who are not into tech startups, this is a personal genomics company that reads your genetic information and sells it back to you. There are a couple of others, like <a href="http://www.decodeme.com/">deCODEme</a>, <a href="http://www.navigenics.com/">Navigenics</a>, and <a href="http://www.bioresolve.com/">BioResolve</a>. Among these, I think 23andme is particularly promising – oh did I mention its cofounder <a href="http://en.wikipedia.org/wiki/Anne_Wojcicki">Anne Wojcicki</a>, is the wife of Google co-founder Sergey Brin?</p> <p>Get back to the genome data topic. This is a typical report you get from company <a href="https://www.23andme.com/">23andme</a>, which is essentially a spreadsheet for all the SNPs that they test on their 2nd generation chip. Each line is a single SNP, with the typical SNP ID (rsid), its location on the reference human genome, and the genotype call. The genotype shows two nucleotides, since we have two copies (from both parents) of each SNP position.</p> <div style="border-bottom: silver 1px solid; text-align: left; border-left: silver 1px solid; padding-bottom: 4px; line-height: 12pt; background-color: #f4f4f4; margin: 20px 0px 10px; padding-left: 4px; width: 97.5%; padding-right: 4px; font-family: 'Courier New', courier, monospace; direction: ltr; max-height: 200px; font-size: 8pt; overflow: auto; border-top: silver 1px solid; cursor: text; border-right: silver 1px solid; padding-top: 4px" id="codeSnippetWrapper"> <div style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px" id="codeSnippet"><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum1"> 1:</span> # rsid chromosome position genotype</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum2"> 2:</span> rs3094315 1 742429 AA</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum3"> 3:</span> rs12562034 1 758311 AG</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum4"> 4:</span> rs3934834 1 995669 CC</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum5"> 5:</span> rs9442372 1 1008567 GG</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum6"> 6:</span> rs3737728 1 1011278 GG</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: white; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum7"> 7:</span> rs11260588 1 1011521 GG</pre><!--CRLF--><pre style="border-bottom-style: none; text-align: left; padding-bottom: 0px; line-height: 12pt; border-right-style: none; background-color: #f4f4f4; margin: 0em; padding-left: 0px; width: 100%; padding-right: 0px; font-family: 'Courier New', courier, monospace; direction: ltr; border-top-style: none; color: black; font-size: 8pt; border-left-style: none; overflow: visible; padding-top: 0px"><span style="color: #606060" id="lnum8"> 8:</span> rs6687776 1 1020428 CC</pre><!--CRLF--></div></div><br /><p>Hyungyong Kim, a bioinformatician working in a Korean biotech firm, had his genome sequenced by 23andme, is open enough to put the <a href="http://biohackers.net/wiki/yong27/Genome">genome on his blog</a> – so that everyone can download it, hack it and know all his weaknesses or strengths. The problem is that you’ll have to read through a ton of literatures to understand what this spreadsheet tells you.</p><br /><p>Enter <a href="http://www.snpedia.com/index.php/Promethease">Promethease</a> – the software that does all the magic – will annotate the SNPs for you. So for Kim’s data, it will report interesting SNPs, drug metabolism, medical conditions, plus more complicated information.</p><br /><p><a href="http://lh3.ggpht.com/_srvRoIok9Xs/TFOccc8-O1I/AAAAAAAABBE/seFmw7i5qV4/s1600-h/Promethease_Kim%5B7%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="Promethease_Kim" border="0" alt="Promethease_Kim" src="http://lh4.ggpht.com/_srvRoIok9Xs/TFOcctO2e3I/AAAAAAAABBI/sK7btGxV5Q4/Promethease_Kim_thumb%5B5%5D.png?imgmax=800" width="583" height="351"></a></p><br /><p>Going down the list. Kim is a male, asian, slow caffeine metabolizer, maybe a few disease susceptibilities here and there. Is this good news or bad news for Kim?</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-30911921820255326382010-07-28T23:50:00.001-04:002010-07-28T23:50:59.135-04:00My dear frond<p>I have been into <a href="http://en.wikipedia.org/wiki/Chaos_game">Chaos game</a> for quite a while. The game proceeds by iteratively generating a sequence of points, using a series of simple linear rewriting rules, or in geek terms – <a href="http://en.wikipedia.org/wiki/Iterated_function_system">iterated function system</a>. Unexpectedly, complex shapes arise with simple formulas, particularly for shapes that are self-similar – trees, leaves, snowflakes, broccoli, coastlines etc. So with that in mind, the following is actually a fractal drawing of fern leaf, versus a real fern leaf. Can you tell the difference?</p> <p><a href="http://lh5.ggpht.com/_srvRoIok9Xs/TFD6nZ8b0nI/AAAAAAAABAw/XRzqU93nE20/s1600-h/ifs%5B13%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="ifs" border="0" alt="ifs" src="http://lh5.ggpht.com/_srvRoIok9Xs/TFD6n3cNbOI/AAAAAAAABA0/oR_2uAkDjYA/ifs_thumb%5B7%5D.png?imgmax=800" width="229" height="269"></a> <a href="http://lh5.ggpht.com/_srvRoIok9Xs/TFD6ocuTsbI/AAAAAAAABA4/Un8F4DY_7d0/s1600-h/fern-leaf-isolated-largethumb1135377%5B1%5D%5B12%5D.jpg"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="fern-leaf-isolated-largethumb1135377[1]" border="0" alt="fern-leaf-isolated-largethumb1135377[1]" src="http://lh6.ggpht.com/_srvRoIok9Xs/TFD6or7XHAI/AAAAAAAABA8/bjDaTRX9BYI/fern-leaf-isolated-largethumb1135377%5B1%5D_thumb%5B10%5D.jpg?imgmax=800" width="220" height="261"></a> </p> <p>For those who are into details, the following is the python code to generate it, using the excellent <a href="http://matplotlib.sourceforge.net/">matplotlib library</a>. Note the re-writing functions and probabilities were taken shamelessly from the <a href="http://en.wikipedia.org/wiki/Iterated_function_system">wiki page</a>.</p><script src="http://gist.github.com/497160.js"> </script> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-70680397092019127282010-07-28T18:32:00.001-04:002010-07-28T18:32:09.520-04:00The world’s first Megatron prototype<p>A <a href="http://www.pnas.org/content/early/2010/06/24/0914069107">recent report</a> on the Proceedings of the National Academy of Sciences by Mit-Geeks et al., reveals a miniaturized machine “shape-shifter” that can achieve two distinct shapes – a boat or a plane. The actual robot is cunningly disguised as a sheet of interconnected triangles. The key innovation is the <strong>actuator – </strong>a device that can bend itself when heated. The <strong>creases</strong> (white stretch on the sheet) are made of silicon thus quite flexible.</p> <p><a href="http://lh6.ggpht.com/_srvRoIok9Xs/TFCv4e3wjrI/AAAAAAAABAc/cjm046dzTnY/s1600-h/megatron%5B1%5D%5B3%5D.jpg"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="megatron[1]" border="0" alt="megatron[1]" src="http://lh5.ggpht.com/_srvRoIok9Xs/TFCv4zyBblI/AAAAAAAABAg/VZS0YlpzxI0/megatron%5B1%5D_thumb%5B1%5D.jpg?imgmax=800" width="288" height="255" /></a> <a href="http://lh6.ggpht.com/_srvRoIok9Xs/TFCv5Nr9qTI/AAAAAAAABAk/mffbH09bCRo/s1600-h/origami%5B8%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="origami" border="0" alt="origami" src="http://lh5.ggpht.com/_srvRoIok9Xs/TFCv5uWS-NI/AAAAAAAABAs/LhWh_XcuuZY/origami_thumb%5B4%5D.png?imgmax=800" width="365" height="254" /></a> </p> <p>The research is the product of a growing field called <strong>computational origami – </strong>which, in addition to fueling my own imagination, have proclaimed applications in machine folding, protein folding, airbag design, among many others. Note the entire operation from the sheet to a targeted geometry can be achieved within 20 seconds in the following video. Un-surprisingly, the research was funded by <a href="http://www.darpa.mil/">DARPA</a>.</p> <object width="480" height="385"><param name="movie" value="http://www.youtube.com/v/Pg8VAVWkz3k&hl=en_US&fs=1"></param><param name="allowFullScreen" value="true"></param><param name="allowscriptaccess" value="always"></param><embed src="http://www.youtube.com/v/Pg8VAVWkz3k&hl=en_US&fs=1" type="application/x-shockwave-flash" allowscriptaccess="always" allowfullscreen="true" width="480" height="385"></embed></object> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-35543505609623603642010-07-13T16:53:00.001-04:002010-07-13T16:53:30.953-04:00Javascript bio-sequence highlighter - I<p>A few weeks back during lunch chat with <a href="http://hackmap.blogspot.com/">brentp</a>, we discussed a simple javascript highlighter for biological sequences, much like the syntax highlighter for programming languages (for example, <a href="http://astyle.sourceforge.net/">Artistic Style</a> for C-like languages or <a href="http://pygments.org/">Pygments</a> for python). The raw sequence files, usually in <a href="http://en.wikipedia.org/wiki/FASTA_format">FASTA format</a>, is just plain text – and very boring to stare at. </p> <p>The biggest issue with biological code is that there are no punctuation marks. Luckily, computational biologists over the last few decades have come up with ways to identify various sequence elements, like genes and repeats. Such <em>annotated</em> sequence features extend over “genomic intervals”, and are typically stored in a GFF file or BED file.</p> <p>The javascript code that I started to write today takes as input a FASTA file and a BED file, and then outputs the HTML codes for the annotated sequence. As a first attempt, it is fairly straightforward, with some preset styles defined in a default CSS file. For example, <em>genes</em> are always plotted with a lavendar background, etc.</p> <p>This code is far from stable, especially stuff that involves nesting of sequence features (e.g. <em>exons</em> should be nested within <em>genes</em>), but right now this is not modeled at all, and resulting in nasty bugs.</p> <p>See a working demo <a href="http://biocon.berkeley.edu/~bao/dna-pygments/">here</a>, and screenshot below. Comments and suggestions are welcome, as always.</p> <p><a href="http://lh3.ggpht.com/_srvRoIok9Xs/TDzSSG6muqI/AAAAAAAAA-I/1_cyO4Ey9EY/s1600-h/js_sequence_highlighter_screenshot%5B4%5D.png"><img style="border-bottom: 0px; border-left: 0px; display: inline; border-top: 0px; border-right: 0px" title="js_sequence_highlighter_screenshot" border="0" alt="js_sequence_highlighter_screenshot" src="http://lh5.ggpht.com/_srvRoIok9Xs/TDzSSsKGiHI/AAAAAAAAA-M/cA7tx_NVePs/js_sequence_highlighter_screenshot_thumb%5B2%5D.png?imgmax=800" width="540" height="343" /></a></p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-15897155399278783382010-07-07T00:14:00.001-04:002010-07-07T00:14:34.594-04:00TimeTree App – evolution in your palm<p>I have rarely looked at <a href="http://www.apple.com/iphone/apps-for-iphone/">iPhone AppStore</a>. The only time I did - I just looked at Top25 or 50 to see what’s hot at the moment. An unexpectedly useful app, <a href="http://itunes.apple.com/us/app/timetree/id372842500?mt=8">TimeTree</a>, turned up when I searched for the keyword <strong>genome </strong>(please don’t ask why I entered this). Here is the description:</p> <blockquote> <p>TimeTree is a public knowledge-base for information on the evolutionary timescale of life. This application allows easy exploration of the thousands of divergence times among organisms in the scientific literature. A tree-based (hierarchical) system is used to identify all published molecular time estimates bearing on the divergence of two chosen organisms, such as species, compute summary statistics, and present the results. Names of two taxa to be compared are entered in the search window and the results are presented on a set of self-explanatory tabs.</p> </blockquote> <p>My own query with <strong>arabidopsis</strong> and <strong>grape</strong> gave me a fast and informative result. </p> <p><a href="http://lh3.ggpht.com/_srvRoIok9Xs/TDP_IFNHmvI/AAAAAAAAA9I/dcTTAn1HSpw/s1600-h/image_2%5B2%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="image_2" border="0" alt="image_2" src="http://lh3.ggpht.com/_srvRoIok9Xs/TDP_IsZpkdI/AAAAAAAAA9M/AbmHhq0yByg/image_2_thumb.png?imgmax=800" width="164" height="244" /></a><a href="http://lh5.ggpht.com/_srvRoIok9Xs/TDP_JFYtWHI/AAAAAAAAA9Q/JeI4QO8DEgU/s1600-h/image%5B2%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="image" border="0" alt="image" src="http://lh5.ggpht.com/_srvRoIok9Xs/TDP_JlHr16I/AAAAAAAAA9U/deIZCTZfhTE/image_thumb.png?imgmax=800" width="164" height="244" /></a> <a href="http://lh3.ggpht.com/_srvRoIok9Xs/TDP_KAIPIHI/AAAAAAAAA9Y/1O_JCuhXSNc/s1600-h/image_1%5B2%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="image_1" border="0" alt="image_1" src="http://lh4.ggpht.com/_srvRoIok9Xs/TDP_Kd-PBdI/AAAAAAAAA9c/WsQX-iQe4uM/image_1_thumb.png?imgmax=800" width="164" height="244" /></a>  </p> <p>The divergence time estimates are mostly based on sequence data, and categorized into nuclear, chloroplast etc, depending on the source. Multiple estimates are weighted (with large variations among different estimates as often expected from this kind of study). References are nicely summarized with  # of genes sampled, data type, table # in the publication… then there is also a reference to <a href="http://www.timetree.org/book.php">the book</a>. I recognized one of the authors - <a href="http://www.kumarlab.net/">Sudmir Kumar</a> wrote the <a href="http://www.megasoftware.net/">MEGA software</a>, a user-friendly GUI tool for exploratory evolutionary analysis.</p> <p>I also played with <strong>escherichia</strong> and <strong>rice</strong> (prokaryote vs. eukaryote), giving me an estimated divergence time of 2622.2Mya. Can you get a longer evolutionary distance than this?</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-23564865141672555452010-04-29T20:30:00.001-04:002010-04-29T20:30:38.627-04:00Get plant gene coordinates from Phytozome<p>For most of my genomics research, all I care is the gene position on the chromosome, without worrying about the exons. As a result, the popular-powerful-intuitive <a href="http://www.sequenceontology.org/gff3.shtml">gff3 format</a> is often over-kill for me. There are common python programming tools for parsing the gff3, including <a href="http://genometools.org/">genometools</a> and Brad Chapman’s <a href="http://www.biopython.org/wiki/GFF_Parsing">BCBio</a>, and they work on well-formed gff3. However, it take some substantial (>10s) time to parse the Arabidopsis TAIR9 gff file (which, by the way, is not standard gff3).</p> <p>The work-around for me, in the last couple of months, along with a colleague in the lab, is to adopt a simpler file format for our calculations and graphics. The format is called <a href="http://genome.ucsc.edu/FAQ/FAQformat.html#format1">“bed” format</a>, and was introduced by UCSC, mainly for graphics purpose. The bed format, like gff, is tab-delimited, but contain no hierarchy among different types. There are less required fields, too. So the bed files I often have contain 4 columns.</p> <pre>Chr2 16848438 16850487 AT2G40340<br />Chr2 14239213 14242365 AT2G33640</pre><br /><br /><p>As it stands, these are all I need for my <a href="http://github.com/tanghaibao/quota-alignment/blob/master/scripts/bed_utils.py">synteny analysis</a>. But is there a easy way to generate bed-formatted files? In the past, I usually generate the bed file by parsing the gff file, either through the library or a custom script. It occurs to me today that I can get all that information from <a href="http://www.phytozome.net/">Phytozome</a>. For a brief introduction, phytozome is a great site that stores many plant gene sequences and their annotations. Their <a href="http://www.phytozome.net/biomart/martview/">biomart site</a> provides an entry point for programmers to get the data out.</p><br /><br /><p>So our script has an xml template. I got the template by going to the <a href="http://www.phytozome.net/biomart/">biomart site</a> and then click on the “XML” button. Then we can send the xml query to the website. Two key things are “filters” and “attributes”. The rest is to send the query through http. You can send the query through “curl” command or through urllib in python.</p><br /><br /><p>Python script for this is attached below. You can specify the species filter to get the data from just a few species (like ‘Arabidopsis thaliana’ in my code). The tab-delimited bed file will be written to stdout.</p><br /><script src="http://gist.github.com/384470.js?file=gistfile1.py"></script> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-24965122669227570862010-03-18T15:33:00.001-04:002010-03-18T15:33:39.000-04:00Review article on Nature Review Genetics retracted<p>I have never seen a review article retracted. Now here is one,</p> <p><strong>Plant genetic engineering for biofuel production: towards affordable cellulosic ethanol : Article : Nature Reviews Genetics (<a href="http://www.nature.com/nrg/journal/v11/n4/full/nrg2777.html">link</a>)</strong></p> <p>So what happened?</p> <blockquote>I am retracting this invited Nature Reviews Genetics article due to a paragraph being paraphrased without attribution. The paragraph in question was from an early version of an article to which I had access as a peer reviewer and which has since been published in Plant Science. </blockquote> <p>Apparently there is some plagiarism going on.. tracking the Plant Science paper. I noticed some Editor’s comment attached to the end. </p> <blockquote>The third paragraph of Section 4 of an early version of this article was previously published, nearly verbatim, in Nature Reviews Genetics, Paragraph 2, Page 441, Volume 9, June 2008 by a reviewer of this manuscript, while this paper was still under review. When this was discovered, the Plant Science Review Editor immediately reported the apparent abuse of the reviewer’s privilege to the academic institution of the author of the paper published in Nature Reviews Genetics. </blockquote> <p>I simply cannot understand this.. my impression for invited review is that they’d have a near-100% acceptance rate. That said, there's not much motivation for plagiarism. I tend to believe that this paragraph gets in the NRG manuscript by some kind of mistake. Although this incident definitely brings disgrace to the NRG author.</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-40684172392010767302010-02-11T22:46:00.001-05:002010-02-16T11:56:42.145-05:00Getting the phylogeny from a list of organisms<p>I started a postdoc job in <a href="http://plantbio.berkeley.edu/~freeling/labweb/welcome.html">this Berkeley lab</a>, and I am quite excited about working with new people. I haven’t posted anything for this year yet, so today I made up something from scratch. </p> <p>Let’s say we have a list of organisms, with their <a href="http://www.ncbi.nlm.nih.gov/Taxonomy/">NCBI taxonomy ids</a>. So how should we know how close they are on the tree of life? Well, you can look it up on the NCBI. For example, below is the lineage for <em><a href="http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=3702">Arabidopsis thaliana</a></em>, the botanical model species. </p> <p><font size="2" face="Courier New">cellular organisms; Eukaryota; Viridiplantae; Streptophyta; Streptophytina; Embryophyta; Tracheophyta; Euphyllophyta; Spermatophyta; Magnoliophyta; eudicotyledons; core eudicotyledons; rosids; malvids; Brassicales; Brassicaceae; Arabidopsis</font> </p> <p>This is very useful indeed. Nowadays people are less devoted to taxonomy than they used to with everything in molecular level, but knowing the lineage can help making phylogenetic inference of when a certain morphological feature arose. </p> <p>How about looking for a handful of species? is there a better way to check one-by-one? You’d need some help from programming. I’ll use Python here and we need a few library dependencies (<a href="http://wwwsearch.sourceforge.net/ClientForm/">ClientForm</a> and <a href="http://www.crummy.com/software/BeautifulSoup/">BeautifulSoup</a>), make sure that you install BeautifulSoup version 3.0.8, because of the problem described <a href="http://www.crummy.com/software/BeautifulSoup/3.1-problems.html">here</a>. </p> <p>Okay I admit. I am lazy. I don’t really want to go to NCBI and grab the long list of lineage and do tons of matches in my code. I found this excellent website “Interactive Tree of Life” (<a href="http://itol.embl.de/other_trees.shtml">iToL</a>) that helps me do this job, so I’ll just deliver my query there. In that, <a href="http://wwwsearch.sourceforge.net/ClientForm/">ClientForm</a> help me do prepare the form and <strong>click</strong> for me, and I’ll use <a href="http://www.crummy.com/software/BeautifulSoup/">BeautifulSoup</a> to retrieve the information I want. The returned data is <a href="http://evolution.genetics.washington.edu/phylip/newicktree.html">Newick-formatted</a>. We’ll use the python library <a href="http://ete.cgenomics.org/">ete</a> to visualize the tree. Full source is below, try to run the example first. If anything is unclear, please see the documentation for the above three python packages (also for installations, more functionalities, etc). </p> <pre style="border-bottom: #cecece 1px solid; border-left: #cecece 1px solid; padding-bottom: 5px; background-color: #fbfbfb; min-height: 40px; padding-left: 5px; width: 650px; padding-right: 5px; overflow: auto; border-top: #cecece 1px solid; border-right: #cecece 1px solid; padding-top: 5px"><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">"<span style="color: #8b0000"></span>""<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">Example:<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">>>> mylist = [3702, 3649, 3694, 3880]<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">>>> t = TaxIDTree(mylist)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">>>> <span style="color: #0000ff">print</span> t<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">(((Carica_papaya,Arabidopsis_thaliana)Brassicales,(Medicago_truncatula,Populus_trichocarpa<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">)fabids)rosids);<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">>>> t.print_tree()<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"><BLANKLINE><br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> /-Carica_papaya<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> /---|<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> | \-Arabidopsis_thaliana<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">---- /---|<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> | /-Medicago_truncatula<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> \---|<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> \-Populus_trichocarpa<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">"<span style="color: #8b0000"></span>""<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">from urllib2 <span style="color: #0000ff">import</span> urlopen<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">from ClientForm <span style="color: #0000ff">import</span> ParseResponse<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">from BeautifulSoup <span style="color: #0000ff">import</span> BeautifulSoup<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px">URL="<span style="color: #8b0000">http://itol.embl.de/other_trees.shtml</span>"<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"><span style="color: #0000ff">class</span> TaxIDTree(object):<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">def</span> __init__(self, list_of_taxids):<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> # the <span style="color: #0000ff">data</span> to send in<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> form_data = "<span style="color: #8b0000">\n</span>".<span style="color: #0000ff">join</span>(str(x) for x in list_of_taxids)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> response = urlopen(URL)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> forms = ParseResponse(response, backwards_compat=False)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> form = forms[0]<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> form["<span style="color: #8b0000">ncbiIDs</span>"] = form_data<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> page = urlopen(form.click()).<span style="color: #0000ff">read</span>()<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> soup = BeautifulSoup(page)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> for element in soup("<span style="color: #8b0000">textarea</span>"):<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">if</span> element["<span style="color: #8b0000">id</span>"]=="<span style="color: #8b0000">nameCol</span>":<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> self.nameCol = str(element.contents[0])<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">def</span> __str__(self):<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">return</span> self.nameCol<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">def</span> print_tree(self):<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> from ete2 <span style="color: #0000ff">import</span> Tree<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> t = Tree(self.nameCol)<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"> <span style="color: #0000ff">print</span> t<br /></pre><pre style="background-color: #fbfbfb; margin: 0em; width: 100%; font-family: consolas,'Courier New',courier,monospace; font-size: 11px"></pre></pre> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-59880279986484693452009-12-07T07:58:00.001-05:002009-12-07T07:58:14.133-05:00gzip.open<p>I read about this <a href="http://blog.revolution-computing.com/2009/12/r-tip-save-time-and-space-by-compressing-data-files.html">blog</a> on how you can save time and space by compressing datafile with gzip and then read in using the zlib library, in R. I did a small test in python, but as it is informal I’ll just state the results. I started with a 730Mb sorghum assembly file (with lots of ACGT), after compressions, the .gz file was around 200Mb. Using normal file read on the raw file took about 5sec on the test server, whereas the .gz file took about 22sec to read. Not impressed.</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-76410774193809391712009-11-14T10:04:00.001-05:002009-11-14T10:04:59.839-05:00The "Gang of Four" in Population Genetics<p>The following are classic theories in the field of population genetics.</p> <ul> <li><a href="http://en.wikipedia.org/wiki/Fisher%27s_fundamental_theorem_of_natural_selection">Fundamental theorem of natural selection</a>, by R. A. Fisher</li> <li><a href="http://en.wikipedia.org/wiki/Shifting_balance_theory">Shifting balance theory</a>, by S. Wright</li> <li><a href="http://www.blackwellpublishing.com/ridley/tutorials/Molecular_evolution_and_neutral_theory7.asp">Genetic load theory</a>, by J. B. S. Haldane</li> <li><a href="http://en.wikipedia.org/wiki/Neutral_theory_of_molecular_evolution">Neutral theory</a>, by M. Kimura</li> </ul> <p>Which one is your favorite?</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-1361911367802511342009-10-29T14:38:00.001-04:002009-10-29T14:38:33.620-04:00Decrease in TAIR funding<p>Although I am not working on Arabidopsis, I think this is still one of the best-managed plant databases out there.</p> <p>Read more from the TAIR homepage and <a href="http://www.arabidopsis.org/doc/about/tair_funding/410">comments from the community</a>. </p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-60223757965541908092009-08-23T22:54:00.002-04:002010-03-01T16:58:43.233-05:00High GC grass genes<p>Several work in the past noted that there are two classes of genes in the grass genome, giving bimodal distribution when plotting the GC content of various genes. The following figure is taken from <a href="http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2270835">Lescot 2008 Musa paper</a>.</p> <p><img src="http://www.pubmedcentral.nih.gov/corecgi/tileshop/tileshop.fcgi?p=PMC3&id=3568&s=14&r=1&c=1" width="424" height="461" /> </p> <p>Please note that the eudicot (Arabidopsis) pattern is quite different from grasses (rice). In some rice genes, in fact there are genes that have the third-codon position virtually all G or C. The third codon positions (also known as synonymous positions) have more freedom to change because the substitutions don’t alter the amino acid composition.From many applications, we need to calculate the <em>Ks</em> distance (substitutions per synonymous positions) as a rough estimate of how divergent two proteins become. Most calculations are based on an <a href="http://en.wikipedia.org/wiki/Substitution_model">evolutionary model</a> of nucleotide bases. With the abnormal gene class noted above, much of the assumptions that simplistic models are no longer valid. In particular, let us assume the <a href="http://en.wikipedia.org/wiki/Models_of_DNA_evolution">JC69 model</a>, which is what Nei and Gojobori’s Ks calculation method is based on. When there are only G<->C changes as allowable transitions, the assumption of equal transition from G to other three nucleotides no longer hold.</p> <p><a href="http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6T39-4JH6C6V-2&_user=655127&_rdoc=1&_fmt=&_orig=search&_sort=d&_docanchor=&view=c&_acct=C000033918&_version=1&_urlVersion=0&_userid=655127&md5=1d929e1ef285b7e9553379b749d89c00">Several published work</a> and I noted in the past that when calculating the <em>Ks</em> for the grass genes, Nei & Gojobori’s method under-estimates the <em>Ks </em>values while PAML over-estimates the <em>Ks</em> values. In practice, we noted that yn00 program in the PAML software, when calculating <em>Ks</em> values between pairs of grass genes, often resulted in values larger than 2. Some further investigations reveal that some results are simply correlated with the protein length! </p> <p>I looked very briefly at the PAML codes and did a bit debugging and found the following code in <em>yn00.c</em>.</p> <pre class="prettyprint">int DistanceF84(double n, double P, double Q, double pi[],<br /> double*k_HKY, double*t, double*SEt)<br />{<br />/* This calculates kappa and d from P (proportion of transitions) & Q<br /> (proportion of transversions) & pi under F84.<br /> When F84 fails, we try to use K80. When K80 fails, we try<br /> to use JC69. When JC69 fails, we set distance t to maxt.<br /> Variance formula under F84 is from Tateno et al. (1994), and briefly<br /> checked against simulated data sets.<br />*/</pre><br /><br /><p>So the code tries to fit the data using several substitution model, and when all fails, switch to the default, so what is the default?</p><br /><br /><pre class="prettyprint">if(failK80) {<br /> if((P+=Q)>=.75) { failJC69=1; P=.75*(n-1.)/n; }<br /> *t = -.75*log(1-P*4/3.);<br /> if(*t>maxt) *t=maxt;<br /> if(SEt) {<br /> *SEt = sqrt(9*P*(1-P)/n) / (3-4*P);<br /> }<br />} <br /></pre><br /><br /><p>Let us look at this code. If transitions (<em>P</em>) plus transversions (<em>Q</em>) are larger than 75%, which means it really is getting completely randomized to the point that it is even more likely to get substituted than not. Under this case, JC69 model fails, fair enough. Then the code enforces the transitions to be slightly less than 75%.</p><br /><br /><p><br /> <br /><a href="http://lh6.ggpht.com/_srvRoIok9Xs/S4w4IGGFwHI/AAAAAAAAA0w/PLySAngrg94/s1600-h/chart%5B1%5D%5B8%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="chart[1]" border="0" alt="chart[1]" src="http://lh4.ggpht.com/_srvRoIok9Xs/S4w4Ic8eRAI/AAAAAAAAA00/uw6zzdOBrY8/chart%5B1%5D_thumb%5B2%5D.png?imgmax=800" width="151" height="24" /></a> ,<a href="http://lh4.ggpht.com/_srvRoIok9Xs/S4w4I1ucQVI/AAAAAAAAA04/gsKjQkqnTcA/s1600-h/chart%5B1%5D%5B5%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="chart[1]" border="0" alt="chart[1]" src="http://lh3.ggpht.com/_srvRoIok9Xs/S4w4JLWWV5I/AAAAAAAAA08/VSM0rLimPOA/chart%5B1%5D_thumb%5B1%5D.png?imgmax=800" width="202" height="24" /></a> </p><br /><br /><p>With a bit of math you can get</p><br /><a href="http://lh3.ggpht.com/_srvRoIok9Xs/S4w4JLZaQ_I/AAAAAAAAA1A/jSp70S4SV3s/s1600-h/chart%5B1%5D%5B2%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="chart[1]" border="0" alt="chart[1]" src="http://lh3.ggpht.com/_srvRoIok9Xs/S4w4JcARf-I/AAAAAAAAA1E/yWGVmt7Z8uo/chart%5B1%5D_thumb.png?imgmax=800" width="105" height="25" /></a> <br /><br /><br /><br /><br /><br /><br /><br /><p>Now you can see that this becomes a very simple function of the sequence length <em>n</em>. I could imagine that for the grass calculation, the three models it tried get thrown off quite often so that this approximation gets called. With <em>n </em>equals to 333 (for a 1kb gene), you have <em>Ks</em> roughly equal 4. That why for a significant portion of the whole dataset, we have a distribution of <em>Ks</em> values that is simply dependent on the sequence length.</p><br />I will emphasize here that this post is not critic to the <a href="http://abacus.gene.ucl.ac.uk/software/paml.html">PAML</a>. I love PAML and use it heavily. But it is important to understand when some unexpected substitutions (such as the grass GC rage) occur, evolutionary models will fail miserably. <br /><br /> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-8374979822480357772009-08-04T14:03:00.001-04:002009-08-04T14:24:58.614-04:00TKF91 model on a tree<div xmlns="http://www.w3.org/1999/xhtml"> <p><object height='350' width='425'><param value="http://youtube.com/v/fmIh-OSt7YU" name="movie" /><embed height="350" width="425" type="application/x-shockwave-flash" src="http://youtube.com/v/fmIh-OSt7YU" /></object></p> <p>This is a bio-informatics related animation on the youtube. It says that the animation was made by "hacking" RASMOL. This is simply amazing and educating at the same time. Naturally, when I looked at the name of the youtube author, it is Ian Holmes.</p> </div> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-56014023007115453382009-07-23T09:32:00.002-04:002009-07-23T09:35:41.040-04:00Facebook Puzzle Master<p>I just discovered <a href="http://www.facebook.com/careers/puzzles.php">this</a> yesterday, and submitted one solution to the test problem. The judge-bot checks the submissions every four hours, so this is unlike some of the online judge (<a href="http://acm.pku.edu.cn/JudgeOnline/">POJ</a>, <a href="http://train.usaco.org/usacogate">USACO)</a> that I worked on before. The problems are more realistic, rather than contrived for competition purposes. One particular problem that I glanced through yesterday -- the <a href="http://www.facebook.com/careers/puzzles.php#/careers/puzzles.php?puzzle_id=1">Facebull</a> problem seems NP-hard (in computer terms, bloody hard), and you cannot solve it in reasonable time for the general case. On top of that, unlike the problems for the online judges, it does not explain the scale of the problem, nor time/memory requirements. The judge-bot accepts popular computer languages. It is always inspiring to see the <a href="http://www.python.org">snake</a>! </p> <p><img src="http://photos-a.ak.fbcdn.net/photos-ak-sf2p/v287/113/16/15325934266/n15325934266_467912_3263.jpg" /> </p> <p>But, the problem I submitted yesterday was judged correct, so I will look at some more problems this evening. I like the online problem listings; and hopefully not limited to the programming. Similar ideas I’ve seen so far include <a href="http://domino.research.ibm.com/Comm/wwwr_ponder.nsf/pages/index.html">mathematics/computing</a> and <a href="https://www.innocentive.com/">chemistry/biology</a>. These are not simple problems but incentives are prepared. </p>Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-67740351685913386562009-07-21T16:05:00.001-04:002009-07-24T22:14:35.433-04:00Plants that I had when I was a young boy<p>There are many indigenous plant that are local varieties of common domesticated plants, but when I was young I didn’t know that. The two plants that I will show below, are perhaps only native to the eastern part of China where I lived. The names are best pronounced in local Shanghai-nese accents. Some pictures shown are taken from random places on the internet.</p> <p>甜卢粟 (sweet sorghum or sugar sorghum; <em>Sorghum bicolor</em>), when I was young, I used to eat that a lot. It is a variety of sorghum, with a high sugar content in the stalk, so a bit like sugarcane, except that the sweet sorghum has a unique aroma (and flavor) to it. The plants grow fast and my mom used to grow these next to the corn plants and the sorghum never needs any fertilizer or other attention. The panicles (the seed bearing part) can be cut and bundled to make a broom. Therefore the seeds must be non-shattering… well I digress. There are some recent interest in developing this into a energy crop, since it has lower water requirements than the sugarcane. See paper <a href="http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6V24-4DWHJ1R-2&_user=655127&_rdoc=1&_fmt=&_orig=search&_sort=d&_docanchor=&view=c&_searchStrId=963969235&_rerunOrigin=google&_acct=C000033918&_version=1&_urlVersion=0&_userid=655127&md5=dae62faf28ff436cdac83e37da159593">here</a>.</p> <p> <a href="http://lh6.ggpht.com/_srvRoIok9Xs/SmpqiAajfMI/AAAAAAAAAcw/KgR8pGR6haY/s1600-h/20070929140127-1986%5B6%5D.jpg"><img style="border-bottom: 0px; border-left: 0px; display: inline; border-top: 0px; border-right: 0px" title="20070929140127-1986" border="0" alt="20070929140127-1986" src="http://lh6.ggpht.com/_srvRoIok9Xs/SmpqivGdpFI/AAAAAAAAAc0/a2Cohoxhnko/20070929140127-1986_thumb%5B4%5D.jpg?imgmax=800" width="305" height="237" /></a> <a href="http://lh5.ggpht.com/_srvRoIok9Xs/SmYfodnZz8I/AAAAAAAAAco/4EDL4h-Tnn0/s1600-h/image%5B3%5D.png"><img style="border-right-width: 0px; display: inline; border-top-width: 0px; border-bottom-width: 0px; border-left-width: 0px" title="image" border="0" alt="image" src="http://lh4.ggpht.com/_srvRoIok9Xs/SmYfoj1n3KI/AAAAAAAAAcs/x3cSFUTj_N8/image_thumb%5B1%5D.png?imgmax=800" width="185" height="240" /></a> </p> <p>癞葡萄 (bitter melon; <em>Momordica charantia</em>), we didn’t grow these. The skin is yellow while the kernels are red and sweet (but not too sugary). When we bought vegetables on the market, some farmers would usually send this to kids (aka me) for no cost. It turned out the in the same genus as the bitter melon, even the same species. I always suspected whether the left is at different developmental stage as the right (which is the bitter melon we eat as vegetables). But the resemblance is obvious.</p> <p><img src="http://img.club.pchome.net/upload/club/other/2008/8/20/pics_mr_child_1219199786.jpg" width="305" height="230" /> <img src="http://upload.wikimedia.org/wikipedia/commons/9/94/Bittermelonfruit.jpg" width="337" height="231" /></p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-58210320366884015722009-07-04T17:38:00.001-04:002009-07-04T17:38:23.477-04:00Mozilla Firefox 3.5 is here<p>I tried it out… noticeably fast (rather informal impression here, considering the network traffic is perhaps lower). In addition, <video> tags are now supported, which means you don’t need a flash plugin to stream videos. The firefox people are now inventing many things; <canvas> for example, has not been supported by IE yet.</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-70692003851578957592009-06-22T16:16:00.001-04:002009-06-22T16:16:03.118-04:00Different millets<p>It has been confusing to me that there are different millet varieties I have encountered when reading literatures (finger millet, foxtail millet and pearl millet). The question again occurred to me when I picked a few weedy grass downstairs and came back to ask Changsoo to ID them. It turns out to be goose grass (<em><a href="http://en.wikipedia.org/wiki/Eleusine_indica">Eleusine indica</a></em>). I looked at some close relatives and found that it is in fact a relative of finger millet (<em><a href="http://en.wikipedia.org/wiki/Finger_millet">Eluesine coracana</a></em>). Now I will try to explain the different millet, copying some images from <a href="http://en.wikipedia.org/wiki/Main_Page">wikipedia</a> to illustrate the morphology. Three major millets (in terms of agricultural production) are listed below.</p> <p><a href="http://en.wikipedia.org/wiki/Pearl_millet">Pearl millet</a> (<em>Pennisetum glaucum, </em>subfamily: <i><a href="http://en.wikipedia.org/wiki/Pennisetum">Pennisetum</a></i>), most widely grown millets, mainly in Africa and India.</p> <p><img src="http://upload.wikimedia.org/wikipedia/commons/thumb/f/f0/Grain_millet,_early_grain_fill,_Tifton,_7-3-02.jpg/180px-Grain_millet,_early_grain_fill,_Tifton,_7-3-02.jpg" width="248" height="189" /> <img src="http://upload.wikimedia.org/wikipedia/commons/thumb/4/43/Pearl_millet_after_combine_harvesting.jpg/400px-Pearl_millet_after_combine_harvesting.jpg" width="245" height="189" /> </p> <p><a href="http://en.wikipedia.org/wiki/Foxtail_millet">Foxtail millet</a> (<em>Setaria italica, </em>subfamily: <a href="http://en.wikipedia.org/wiki/Panicoideae"><em>Panicoideae</em></a>), very important in east asia, in fact been grown in China for more than 8000 years (there is a recent <em>PNAS</em> paper on this). The Chinese name for this is <strong>小米</strong>; it also looks like a local weedy grass <em>Setaria viridus</em>, which we called <strong>狗尾巴草</strong>.</p> <p><img src="http://upload.wikimedia.org/wikipedia/commons/thumb/2/28/Foxtailmillet.jpg/200px-Foxtailmillet.jpg" width="248" height="180" /> <img src="http://www.3nong.cc/admin/Shopimages/2009514114645小米.jpg" width="207" height="180" /> </p> <p><a href="http://en.wikipedia.org/wiki/Finger_millet">Finger millet</a> (<i>Eleusine coracana</i>, subfamily: <a href="http://en.wikipedia.org/wiki/Chloridoideae"><em>Chloridoideae</em></a>), annual plant grown as cereal in Africa and Asia. </p> <p><img src="http://upload.wikimedia.org/wikipedia/en/thumb/6/6c/Finger_millet_3_11-21-02.jpg/280px-Finger_millet_3_11-21-02.jpg" width="242" height="183" /> <img src="http://upload.wikimedia.org/wikipedia/commons/c/c6/Finger_millet_grains_of_mixed_color.jpg" width="220" height="184" /></p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com1tag:blogger.com,1999:blog-8927778.post-7395406602398242002008-08-18T22:06:00.001-04:002008-08-18T22:06:47.012-04:00STL usage<pre>发信人: littlebee (小蜜蜂), 信区: C<br />标 题: STL的技巧<br />发信站: 日月光华 (2008年08月07日16:36:25 星期四)<br /><br />网上看到的,都是一些常用技巧,都不是很难,不过有些确实很有用。。<br />转过来大家看看吧~~<br /><br /><br />toupper,tolower<br />地球人都知道 C++ 的 string 没有 toupper ,好在这不是个大问题,因为我们有 STL 算<br />法:<br /><br />string s("heLLo");<br />transform(s.begin(), s.end(), s.begin(), ::toupper);<br />cout << s << endl;<br />transform(s.begin(), s.end(), s.begin(), ::tolower);<br />cout << s << endl;<br /><br />当然,我知道很多人希望的是 s.to_upper() ,但是对于一个这么通用的 basic_string <br />来说,的确没办法把这些专有的方法放进来。如果你用 boost stringalgo ,那当然不在<br />话下,你也就不需要读这篇文章了。<br /><br />------------------------------------------------------------------------<br />trim<br />我们还知道 string 没有 trim ,不过自力更生也不困难,比 toupper 来的还要简单:<br /><br /><br /> string s(" hello ");<br /> s.erase(0, s.find_first_not_of(" \n"));<br /> cout << s << endl;<br /> s.erase(s.find_last_not_of(' ') + 1);<br /> cout << s << endl;<br /><br />注意由于 find_first_not_of 和 find_last_not_of 都可以接受字符串,这个时候它们寻<br />找该字符串中所有字符的 absence ,所以你可以一次 trim 掉多种字符。<br /><br />-----------------------------------------------------------------------<br />erase<br />string 本身的 erase 还是不错的,但是只能 erase 连续字符,如果要拿掉一个字符串里<br />面所有的某个字符呢?用 STL 的 erase + remove_if 就可以了,注意光 remove_if 是不<br />行的。<br /><br /> string s(" hello, world. say bye ");<br /> s.erase(remove_if(s.begin(),s.end(), <br /> bind2nd(equal_to<char>(), ' ')), <br /> s.end());<br /><br />上面的这段会拿掉所有的空格,于是得到 hello,world.saybye。<br /><br />-----------------------------------------------------------------------<br />replace<br />string 本身提供了 replace ,不过并不是面向字符串的,譬如我们最常用的把一个 sub<br />str 换成另一个 substr 的操作,就要做一点小组合:<br /><br /> string s("hello, world");<br /> string sub("ello, ");<br /> s.replace(s.find(sub), sub.size(), "appy ");<br /> cout << s << endl;<br /><br />输出为 happy world。注意原来的那个 substr 和替换的 substr 并不一定要一样长。<br /><br /><br />-----------------------------------------------------------------------<br />startwith, endwith<br />这两个可真常用,不过如果你仔细看看 string 的接口,就会发现其实没必要专门提供这<br />两个方法,已经有的接口可以干得很好:<br /><br /> string s("hello, world");<br /> string head("hello");<br /> string tail("ld");<br /> bool startwith = s.compare(0, head.size(), head) == 0;<br /> cout << boolalpha << startwith << endl;<br /> bool endwith = s.compare(s.size() - tail.size(), tail.size(), tail) == 0;<br /><br /> cout << boolalpha << endwith << endl;<br /><br />当然了,没有 s.startwith("hello") 这样方便。<br /><br />------------------------------------------------------------------------<br />toint, todouble, tobool...<br />这也是老生常谈了,无论是 C 的方法还是 C++ 的方法都可以,各有特色:<br /><br /> string s("123");<br /> int i = atoi(s.c_str());<br /> cout << i << endl;<br /> <br /> int ii;<br /> stringstream(s) >> ii;<br /> cout << ii << endl;<br /> <br /> string sd("12.3");<br /> double d = atof(sd.c_str());<br /> cout << d << endl;<br /> <br /> double dd;<br /> stringstream(sd) >> dd;<br /> cout << dd << endl;<br /> <br /> string sb("true");<br /> bool b;<br /> stringstream(sb) >> boolalpha >> b;<br /> cout << boolalpha << b << endl;<br /><br />C 的方法很简洁,而且赋值与转换在一句里面完成,而 C++ 的方法很通用。<br /><br />------------------------------------------------------------------------<br />split<br />这可是件麻烦事,我们最希望的是这样一个接口: s.split(vect, ',') 。用 STL 算法来<br />做有一定难度,我们可以从简单的开始,如果分隔符是空格、tab 和回车之类,那么这样<br />就够了:<br /><br /> string s("hello world, bye.");<br /> vector<string> vect;<br /> vect.assign(<br /> istream_iterator<string>(stringstream(s)),<br /> istream_iterator<string>()<br /> );<br /><br />不过要注意,如果 s 很大,那么会有效率上的隐忧,因为 stringstream 会 copy 一份 <br />string 给自己用。<br /><br />------------------------------------------------------------------------<br />concat<br />把一个装有 string 的容器里面所有的 string 连接起来,怎么做?希望你不要说是 han<br />d code 循环,这样做不是更好?<br /><br /> vector<string> vect;<br /> vect.push_back("hello");<br /> vect.push_back(", ");<br /> vect.push_back("world");<br /> <br /> cout << accumulate(vect.begin(), vect.end(), string(""));<br /><br />不过在效率上比较有优化余地。<br /><br />-------------------------------------------------------------------------<br /><br />reverse<br />其实我比较怀疑有什么人需要真的去 reverse 一个 string ,不过做这件事情的确是很容<br />易:<br /><br /> std::reverse(s.begin(), s.end());<br /><br />上面是原地反转的方法,如果需要反转到别的 string 里面,一样简单:<br /><br /> s1.assign(s.rbegin(), s.rend());<br /><br />效率也相当理想。<br /><br />-------------------------------------------------------------------------<br /><br />解析文件扩展名<br />字数多点的写法:<br /><br /> std::string filename("hello.exe");<br /><br /> std::string::size_type pos = filename.rfind('.');<br /> std::string ext = filename.substr(pos == std::string::npos ? filename.leng<br />th() : pos + 1);<br /><br />不过两行,合并成一行呢?也不是不可以:<br /><br /> std::string ext = filename.substr(filename.rfind('.') == std::string::npos<br /> ? filename.length() : filename.rfind('.') + 1);<br /><br />我知道,rfind 执行了两次。不过第一,你可以希望编译器把它优化掉,其次,扩展名一<br />般都很短,即便多执行一次,区别应该是相当微小。 <br />STL 算法 <br />distance <br />很多时候我们希望在一个 vector ,或者 list ,或者什么其他东西里面,找到一个值在<br />哪个位置,这个时候 find 帮不上忙,而有人就转而求助手写循环了,而且是原始的手写<br />循环:<br /><br />for ( int i = 0; i < vect.size(); ++i)<br /> if ( vect[i] == value ) break;<br /><br />如果编译器把 i 看作 for scope 的一部分,你还要把 i 的声明拿出去。真的需要这样么<br />?看看这个:<br /><br /> int dist = <br /> distance(col.begin(), <br /> find(col.begin(), col.end(), 5));<br /><br />其中 col 可以是很多容器,list, vector, deque... 当然这是你确定 5 就在 col 里面<br />的情形,如果你不确定,那就加点判断:<br /><br /> int dist;<br /> list<int>::iterator pos = find(col.begin(), col.end(), 5);<br /> if ( pos != col.end() )<br /> dist = distance(col.begin(), pos);<br /><br />我想这还是比手写循环来的好些吧。<br /><br />--------------------------------------------------------------------------<br />max, min<br />这是有直接的算法支持的,当然复杂度是 O(n),用于未排序容器,如果是排序容器...老<br />兄,那还需要什么算法么?<br /><br />max_element(col.begin(), col.end());<br />min_element(col.begin(), col.end());<br /><br />注意返回的是 iterator ,如果你关心的只是值,那么好:<br /><br />*max_element(col.begin(), col.end());<br />*min_element(col.begin(), col.end());<br /><br />max_element 和 min_element 都默认用 less 来排序,它们也都接受一个 binary predi<br />cate ,如果你足够无聊,甚至可以把 max_element 当成 min_element 来用,或者反之:<br /><br /><br />*max_element(col.begin(), col.end(), greater<int>()); // 返回最小值!<br />*min_element(col.begin(), col.end(), greater<int>()); // 返回最大值<br /><br />当然它们的本意不是这个,而是让你能在比较特殊的情况下使用它们,例如,你要比较的<br />是每个元素的某个成员,或者成员函数的返回值。例如:<br /><br />#include <iostream><br />#include <list><br />#include <algorithm><br />#include <string><br />#include <boost/bind.hpp><br /><br />using namespace boost;<br />using namespace std;<br /><br />struct Person<br />{<br /> Person(const string& _name, int _age)<br /> : name(_name), age(_age)<br /> {}<br /> int age;<br /> string name;<br />};<br /><br />int main()<br />{<br /> list<Person> col;<br /> list<Person>::iterator pos;<br /><br /> col.push_back(Person("Tom", 10));<br /> col.push_back(Person("Jerry", 12));<br /> col.push_back(Person("Mickey", 9));<br /><br /> Person eldest = <br /> *max_element(col.begin(), col.end(), <br /> bind(&Person::age, _1) < bind(&Person::age, _2));//>=1.33<br /> <br /> cout << eldest.name;<br />}<br /><br />输出是 Jerry ,这里用了 boost.bind ,原谅我不知道用 bind2nd, mem_fun 怎么写,我<br />也不想知道...<br /><br />-------------------------------------------------------------------------<br />copy_if<br />没错,STL 里面压根没有 copy_if ,这就是为什么我们需要这个:<br /><br />template<typename InputIterator, typename OutputIterator, typename Predicate><br /><br />OutputIterator copy_if(<br /> InputIterator begin, InputIterator end, OutputIterator destBegin, Predicat<br />e p)<br />{<br /> while (begin != end) <br /> {<br /> if (p(*begin))*destBegin++ = *begin;<br /> ++begin;<br /> }<br /> return destBegin;<br />}<br /><br />把它放在自己的工具箱里,是一个明智的选择。<br /><br />------------------------------------------------------------------------<br />惯用手法:erase(iter++)<br />如果你要去除一个 list 中的某些元素,那可千万小心:(下面的代码是错的!!!)<br /><br /><br />#include <iostream><br />#include <algorithm><br />#include <iterator><br />#include <list><br /><br />int main()<br />{<br /> int arr[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};<br /> std::list<int> lst(arr, arr + 10);<br /><br /> for ( std::list<int>::iterator iter = lst.begin();<br /> iter != lst.end(); ++iter)<br /> if ( *iter % 2 == 0 )<br /> lst.erase(iter);<br /> <br /> std::copy(lst.begin(), lst.end(),<br /> std::ostream_iterator<int>(std::cout, " "));<br />}<br /><br />当 iter 被 erase 掉的时候,它已经失效,而后面却还会做 ++iter ,其行为无可预期!<br />如果你不想动用 remove_if ,那么唯一的选择就是:<br /><br />#include <iostream><br />#include <algorithm><br />#include <iterator><br />#include <list><br /><br />int main()<br />{<br /> int arr[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};<br /> std::list<int> lst(arr, arr + 10);<br /><br /> for ( std::list<int>::iterator iter = lst.begin();<br /> iter != lst.end(); )<br /> if ( *iter % 2 == 0 )<br /> lst.erase(iter++);<br /> else<br /> ++iter;<br /> <br /> std::copy(lst.begin(), lst.end(),<br /> std::ostream_iterator<int>(std::cout, " "));<br />}<br /><br />但是上面的代码不能用于 vector, string 和 deque ,因为对于这些容器, erase 不光<br />令 iter 失效,还令 iter 之后的所有 iterator 失效!<br /><br />-------------------------------------------------------------------------<br />erase(remove...) 惯用手法<br />上面的循环如此难写,如此不通用,如此不容易理解,还是用 STL 算法来的好,但是注意<br />,光 remove_if 是没用的,必须使用 erase(remove...) 惯用手法:<br /><br />#include <iostream><br />#include <algorithm><br />#include <iterator><br />#include <list><br />#include <functional><br />#include <boost/bind.hpp><br /><br />int main()<br />{<br /> int arr[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};<br /> std::list<int> lst(arr, arr + 10);<br /><br /> lst.erase(remove_if(lst.begin(), lst.end(),<br /> boost::bind(std::modulus<int>(), _1, 2) == 0),<br /> lst.end()<br /> );<br /> <br /> std::copy(lst.begin(), lst.end(),<br /> std::ostream_iterator<int>(std::cout, " "));<br />}<br /><br />当然,这里借助了 boost.bind ,让我们不用多写一个没用的 functor 。 <br />简单常识——关于stream <br />从文件中读入一行 <br /><br />简单,这样就行了: <br /><br />ifstream ifs("input.txt");<br />char buf[1000];<br /><br /><br />ifs.getline(buf, sizeof buf); <br /><br />string input(buf); <br /><br />当然,这样没有错,但是包含不必要的繁琐和拷贝,况且,如果一行超过1000个字符,就<br />必须用一个循环和更麻烦的缓冲管理。下面这样岂不是更简单?<br /><br />string input;<br />input.reserve(1000);<br />ifstream ifs("input.txt");<br />getline(ifs, input); <br /><br />不仅简单,而且安全,因为全局函数 getline 会帮你处理缓冲区用完之类的麻烦,如果你<br />不希望空间分配发生的太频繁,只需要多 reserve 一点空间。<br /><br />这就是“简单常识”的含义,很多东西已经在那里,只是我一直没去用。<br /><br />---------------------------------------------------------------------------<br /><br />一次把整个文件读入一个 string <br /><br /><br />我希望你的答案不要是这样:<br /><br />string input;<br />while( !ifs.eof() )<br />{<br /> string line;<br /> getline(ifs, line);<br /> input.append(line).append(1, '\n');<br />} <br /><br />当然了,没有错,它能工作,但是下面的办法是不是更加符合 C++ 的精神呢?<br /><br />string input(<br /> istreambuf_iterator<char>(instream.rdbuf()), <br /> istreambuf_iterator<char>()<br />); <br /><br />同样,事先分配空间对于性能可能有潜在的好处:<br /><br />string input;<br />input.reserve(10000);<br />input.assign(<br /> istreambuf_iterator<char>(ifs.rdbuf()), <br /> istreambuf_iterator<char>()<br />);<br /><br /><br />很简单,不是么?但是这些却是我们经常忽略的事实。<br />补充一下,这样干是有问题的:<br /><br /> string input; <br /> input.assign( <br /> istream_iterator<char>(ifs), <br /> istream_iterator<char>() <br /> ); <br /><br /><br />因为它会忽略所有的分隔符,你会得到一个纯“字符”的字符串。最后,如果你只是想把<br />一个文件的内容读到另一个流,那没有比这更快的了:<br /><br /><br /> fstream fs("temp.txt"); <br /> cout << fs.rdbuf(); <br /><br />因此,如果你要手工 copy 文件,这是最好的(如果不用操作系统的 API):<br /><br /> ifstream ifs("in.txt"); <br /> ofstream ofs("out.txt"); <br /> ofs << in.rdbuf(); <br /><br /><br />-------------------------------------------------------------------------<br /><br />open 一个文件的那些选项 <br /><br />ios::in Open file for reading <br />ios::out Open file for writing <br />ios::ate Initial position: end of file <br />ios::app Every output is appended at the end of file <br />ios::trunc If the file already existed it is erased <br />ios::binary Binary mode <br /><br />-------------------------------------------------------------------------<br /><br />还有 ios 的那些 flag <br /><br />flag effect if set <br />ios_base::boolalpha input/output bool objects as alphabetic names (true, fals<br />e). <br />ios_base::dec input/output integer in decimal base format. <br />ios_base::fixed output floating point values in fixed-point notation. <br />ios_base::hex input/output integer in hexadecimal base format. <br />ios_base::internal the output is filled at an internal point enlarging the ou<br />tput up to the field width. <br />ios_base::left the output is filled at the end enlarging the output up to the<br /> field width. <br />ios_base::oct input/output integer in octal base format. <br />ios_base::right the output is filled at the beginning enlarging the output up<br /> to the field width. <br />ios_base::scientific output floating-point values in scientific notation. <br /><br />ios_base::showbase output integer values preceded by the numeric base. <br />ios_base::showpoint output floating-point values including always the decimal<br /> point. <br />ios_base::showpos output non-negative numeric preceded by a plus sign (+). <br /><br />ios_base::skipws skip leading whitespaces on certain input operations. <br />ios_base::unitbuf flush output after each inserting operation. <br />ios_base::uppercase output uppercase letters replacing certain lowercase lett<br />ers. <br /><br />There are also defined three other constants that can be used as masks: <br /><br />constant value <br />ios_base::adjustfield left | right | internal <br />ios_base::basefield dec | oct | hex <br />ios_base::floatfield scientific | fixed <br /><br />--------------------------------------------------------------------------<br /><br /><br />用我想要的分隔符来解析一个字符串,以及从流中读取数据 <br /><br /><br />这曾经是一个需要不少麻烦的话题,由于其常用而显得尤其麻烦,但是其实 getline 可以<br />做得不错:<br /><br /><br /> getline(cin, s, ';'); <br /> while ( s != "quit" ) <br /> { <br /> cout << s << endl; <br /> getline(cin, s, ';'); <br /> } <br /><br /><br />简单吧?不过注意,由于这个时候 getline 只把 ; 作为分隔符,所以你需要用 ;quit; <br />来结束输入,否则 getline 会把前后的空格和回车都读入 s ,当然,这个问题可以在代<br />码里面解决。<br /><br /><br />同样,对于简单的字符串解析,我们是不大需要动用什么 Tokenizer 之类的东西了:<br /><br /><br />#include <iostream> <br />#include <sstream> <br />#include <string> <br /><br />using namespace std; <br /><br />int main() <br />{ <br /> string s("hello,world, this is a sentence; and a word, end."); <br /> stringstream ss(s); <br /> <br /> for ( ; ; ) <br /> { <br /> string token; <br /> getline(ss, token, ','); <br /> if ( ss.fail() ) break; <br /> <br /> cout << token << endl; <br /> } <br />} <br /><br /><br />输出:<br /><br /><br />hello <br />world <br /> this is a sentence; and a word <br /> end. <br /><br /><br />很漂亮不是么?不过这么干的缺陷在于,只有一个字符可以作为分隔符。<br /><br /><br />--------------------------------------------------------------------------<br /><br /><br />把原本输出到屏幕的东西输出到文件,不用到处去把 cout 改成 fs<br /><br /><br />#include <iostream><br />#include <fstream> <br />using namespace std; <br />int main()<br />{ <br /> ofstream outf("out.txt"); <br /> streambuf *strm_buf=cout.rdbuf(); <br /> cout.rdbuf(outf.rdbuf()); <br /> cout<<"write something to file"<<endl; <br /> cout.rdbuf(strm_buf); //recover <br /> cout<<"display something on screen"<<endl; <br /> system("PAUSE");<br /> return 0;<br />} <br /> <br />输出到屏幕的是:<br /><br /><br />display something on screen <br /><br /><br />输出到文件的是:<br /><br /><br />write something to file <br /><br /><br />也就是说,只要改变 ostream 的 rdbuf ,就可以重定向了,但是这招对 fstream 和 st<br />ringstream 都没用。<br /><br /><br />--------------------------------------------------------------------------<br /><br /><br />关于 istream_iterator 和 ostream_iterator<br /><br /><br />经典的 ostream_iterator 例子,就是用 copy 来输出:<br /><br /><br />#include <iostream> <br />#include <fstream> <br />#include <sstream> <br />#include <algorithm> <br />#include <vector> <br />#include <iterator> <br /><br />using namespace std; <br /><br />int main() <br />{ <br /> vector<int> vect; <br /> for ( int i = 1; i <= 9; ++i ) <br /> vect.push_back(i); <br /> <br /> copy(vect.begin(), vect.end(), <br /> ostream_iterator<int>(cout, " ") <br /> ); <br /> cout << endl; <br /> <br /> ostream_iterator<double> os_iter(cout, " ~ "); <br /> *os_iter = 1.0; <br /> os_iter++; <br /> *os_iter = 2.0; <br /> *os_iter = 3.0; <br />} <br /><br />输出:<br /><br /><br />1 2 3 4 5 6 7 8 9 <br />1 ~ 2 ~ 3 ~ <br /><br /><br />很明显,ostream_iterator 的作用就是允许对 stream 做 iterator 的操作,从而让算法<br />可以施加于 stream 之上,这也是 STL 的精华。与前面的“读取文件”相结合,我们得到<br />了显示一个文件最方便的办法:<br /><br /><br /> copy(istreambuf_iterator<char>(ifs.rdbuf()), <br /> istreambuf_iterator<char>(), <br /> ostreambuf_iterator<char>(cout) <br /> ); <br /><br /><br />同样,如果你用下面的语句,得到的会是没有分隔符的输出:<br /><br /><br /> copy(istream_iterator<char>(ifs), <br /> istream_iterator<char>(), <br /> ostream_iterator<char>(cout) <br /> ); <br /><br /><br />那多半不是你要的结果。如果你硬是想用 istream_iterator 而不是 istreambuf_iterat<br />or 呢?还是有办法:<br /><br /><br /> copy(istream_iterator<char>(ifs >> noskipws), <br /> istream_iterator<char>(), <br /> ostream_iterator<char>(cout) <br /> ); <br /><br /><br />但是这样不是推荐方法,它的效率比第一种低不少。<br />如果一个文件 temp.txt 的内容是下面这样,那么我的这个从文件中把数据读入 vector <br />的方法应该会让你印象深刻。<br /><br /><br />12345 234 567<br />89 10<br /><br /><br />程序:<br /><br /><br />#include <iostream> <br />#include <fstream> <br />#include <algorithm> <br />#include <vector> <br />#include <iterator> <br /><br />using namespace std; <br /><br />int main() <br />{ <br /> ifstream ifs("temp.txt"); <br /> <br /> vector<int> vect; <br /> vect.assign(istream_iterator<int>(ifs),<br /> istream_iterator<int>()<br /> ); <br /><br /> copy(vect.begin(), vect.end(), ostream_iterator<int>(cout, " ")); <br />} <br /><br />输出:<br /><br /><br />12345 234 567 89 10 <br /><br /><br />很酷不是么?判断文件结束、移动文件指针之类的苦工都有 istream_iterator 代劳了。<br /><br /><br /><br />-----------------------------------------------------------------------<br /><br /><br />其它算法配合 iterator <br /><br /><br />计算文件行数:<br /><br /><br /> int line_count = <br /> count(istreambuf_iterator<char>(ifs.rdbuf()), <br /> istreambuf_iterator<char>(), <br /> '\n'); <br /><br /><br />当然确切地说,这是在计算文件中回车符的数量,同理,你也可以计算文件中任何字符的<br />数量,或者某个 token 的数量:<br /><br /><br /> int token_count = <br /> count(istream_iterator<string>(ifs), <br /> istream_iterator<string>(), <br /> "#include"); <br /><br /><br />注意上面计算的是 “#include” 作为一个 token 的数量,如果它和其他的字符连起来,<br />是不算数的。<br /><br /><br />------------------------------------------------------------------------<br />Manipulator<br /><br /><br />Manipulator 是什么?简单的说,就是一个接受一个 stream 作为参数,并且返回一个 s<br />tream 的函数,比如上面的 unskipws ,它的定义是这样的:<br /><br /><br /> inline ios_base& <br /> noskipws(ios_base& __base) <br /> { <br /> __base.unsetf(ios_base::skipws); <br /> return __base; <br /> } <br /><br /><br />这里它用了更通用的 ios_base 。知道了这一点,你大概不会对自己写一个 manipulator<br /> 有什么恐惧感了,下面这个无聊的 manipulator 会忽略 stream 遇到第一个分号之前所<br />有的输入(包括那个分号):<br /><br /><br />template <class charT, class traits><br />inline std::basic_istream<charT, traits>&<br />ignoreToSemicolon (std::basic_istream<charT, traits>& s)<br />{<br /> s.ignore(std::numeric_limits<int>::max(), s.widen(';'));<br /> return s;<br />}<br /><br />不过注意,它不会忽略以后的分号,因为 ignore 只执行了一次。更通用一点,manipula<br />tor 也可以接受参数的,下面这个就是 ignoreToSemicolon 的通用版本,它接受一个参数<br />, stream 会忽略遇到第一个该参数之前的所有输入,写起来稍微麻烦一点:<br /><br /><br />struct IgnoreTo {<br /> char ignoreTo;<br /> IgnoreTo(char c) : ignoreTo(c) <br /> {}<br />};<br /> <br />std::istream& operator >> (std::istream& s, const IgnoreTo& manip)<br />{<br /> s.ignore(std::numeric_limits<int>::max(), s.widen(manip.ignoreTo)); <br /> return s;<br />}<br /><br />但是用法差不多:<br /><br /><br /> copy(istream_iterator<char>(ifs >> noskipws >> IgnoreTo(';')), <br /> istream_iterator<char>(), <br /> ostream_iterator<char>(cout) <br /> ); <br /><br /><br />其效果跟 IgnoreToSemicolon 一样。</pre> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-64931124631625023022008-08-10T17:05:00.001-04:002008-08-10T17:05:30.362-04:00Stanford machine learning series<p>This course provides a broad introduction to machine learning and statistical pattern recognition. Topics include supervised learning, unsupervised learning, learning theory, reinforcement learning and adaptive control. Recent applications of machine learning, such as to robotic control, data mining, autonomous navigation, bioinformatics, speech recognition, and text and web data processing are also discussed. <br />Complete Playlist for the Course: <br /><a href="http://www.youtube.com/view_play_list?p=A89DCFA6ADACE599">http://www.youtube.com/view_play_list...</a> <br />CS 229 Course Website: <br /><a href="http://www.stanford.edu/class/cs229/">http://www.stanford.edu/class/cs229/</a></p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-59114473113896477752008-05-27T16:08:00.001-04:002009-07-25T01:02:22.820-04:00Briefings in bioinformatics -- how impact factor should be calculated<p>The review-oriented journal has received its first impact factor measurement from Journal citation reports for 2006, an astounding 24.37; and this would rank as the first in the field of bioinformatics. I have looked at some articles in this journal, including a recent review on operon prediction -- but this has turned out to be quite a surprise to me. A closer look at the announcement both on the <a href="http://bib.oxfordjournals.org/">journal website</a> and <a href="http://bib.oxfordjournals.org/cgi/content/full/8/4/207">an editorial</a> reveals that <a href="http://bib.oxfordjournals.org/cgi/content/abstract/5/2/150">one article</a> on phylogenetic reconstruction software MEGA3 has contributed significantly to the high score. </p> <p>The editors of BIB acknowledge this outlier and explicitly point out that once the above article is removed, impact factor drops to a 4 -- still considered a relatively high score for a new magazine.</p> <p>While appreciating the editors' honesty, one should look at the formula of IF calculation</p> <p><img src="http://latex.codecogs.com/gif.latex?2006\ \textup{IMPACT FACTOR}=\frac{\textup{Cites in 2006 to articles published in 2004 and 2005}}{\textup{Total number of articles published in 2004 and 2005}}" /> </p> <p>For a journal that contains quite few articles, the large number of citations of one or two articles could affect a lot on the score. I imagine the degree (citations) distribution follows power law, in that a few papers actually attracted most citations.</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-80934728983863088242008-05-01T10:26:00.001-04:002008-05-01T10:26:33.875-04:00Lenovo T61 woes<p>I have waited for the new computer for about half a month now, only to discover that the Intel 3945A/B/G is not compatible with my Linksys router. Internet connection drops to about one-fifth of what I used to get. All the symptoms seem to go away as soon as I plug in the cable directly. Played around with it for more than three hour last night, tried a dozen solutions suggested but still no luck. </p> <p> </p> <p>Finally I give up and plug in my old smart card Linksys NIC, everything seems to be OK now. I wish I had sticked to the Macintosh ...</p> Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0tag:blogger.com,1999:blog-8927778.post-49758099412876011942008-03-27T10:16:00.005-04:002008-04-14T03:40:36.740-04:00Knuth shuffle (Fisher-Yates shuffle)<p>The shuffling algorithm occurred to me several years ago when I first learned programming, but recently re-surfaced as I tried to digest a piece of Perl code from lab note for PBIO class.</p> <p>The following is the perl code, </p><pre class="prettyprint">#!/usr/bin/perl<br />use warnings;<br />use strict;<br />open IFILE, 'sequence.txt';<br />chomp (my $seq = <IFILE>);<br />close IFILE;<br />my @oldseq = split('', $seq);<br />srand;<br />my @newseq = ();<br />for( @oldseq ){<br />my $r = rand(@newseq+1);<br />push(@newseq,$newseq[$r]);<br />$newseq[$r] = $_;<br />}<br />print @oldseq, "\n";<br />print @newseq, "\n";<br /></pre><br /><p>It tries to loop through the original array and then exchange the element to an element already passed through (including itself), this would generate n! execution paths, and the possible outcomes are n!, and it is easy to prove that there is no collision in different execution path.</p> <p>Note that this is a slight variant of Knuth shuffle, where the elements are swapped with an element that has NOT passed through (including itself), the same number of execution paths. Optimal solution. O(n log n) time complexity.</p> <p>A more intuitive way (to me) is to think about it in Fisher-Yates original method, where a random item is taken each time out of a hat. Note that there are potential pitfalls as well when implementing this simple permutation algorithm.</p> <p>Reference:</p> <p><a title="http://en.wikipedia.org/wiki/Fisher-Yates_shuffle" href="http://en.wikipedia.org/wiki/Fisher-Yates_shuffle">http://en.wikipedia.org/wiki/Fisher-Yates_shuffle</a></p> <p><a title="http://en.wikipedia.org/wiki/Shuffle" href="http://en.wikipedia.org/wiki/Shuffle">http://en.wikipedia.org/wiki/Shuffle</a></p>Anonymoushttp://www.blogger.com/profile/13607944791843532178noreply@blogger.com0