tag:blogger.com,1999:blog-2662347345150434102024-03-06T03:40:19.075-05:00a frictionless surfaceetalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.comBlogger30125tag:blogger.com,1999:blog-266234734515043410.post-23933331261550172382017-02-02T16:39:00.000-05:002017-08-21T12:57:41.689-04:00What is the UCSF500?<a href="http://precisionmedicine.ucsf.edu/">Precision medicine</a> is more than genome sequencing, but molecular profiling is an essential part of it. Recognizing this, in 2014 UC San Francisco launched the Genomic Medicine Initiative to bring high-throughput DNA and RNA sequencing techniques into routine clinical care. The first product of this effort is the UCSF500, a targeted cancer genome sequencing service that is now available to patients at the UCSF Medical Center.<br />
<!--
Advancing precision medicine is a priority at UCSF and is central to its
advancing health worldwide mission. UCSF faculty lead efforts at the
local, state and national levels to use data-driven tools and analysis to
develop new diagnostics, therapies and insights into disease. The goal is
to make progress in both personal and population health.
--><br />
The UCSF500 service is provided by the Clinical Cancer Genomic Lab (CCGL). This group is directed by <a href="http://cancer.ucsf.edu/people/profiles/bastian_boris.3307">Dr. Boris C. Bastian</a>, and the informatics group, where I work, is led by <a href="http://profiles.ucsf.edu/iwei.yeh">Dr. Iwei Yeh</a>. The pilot program from 2014 to 2016 focused on local patients with metastatic disease, including children and patients with rare or poorly understood cancer types.<br />
<br />
Several published studies from this program have already given us insight into cancer mechanisms and treatment options:<br />
<ul>
<li><a href="http://dx.doi.org/10.1093/neuonc/now254">Childhood gliomas</a> and other brain tumors</li>
<li><a href="http://dx.doi.org/10.1038/modpathol.2016.97">Malignant phyllodes tumors</a> of the breast</li>
<li><a href="http://dx.doi.org/10.1038/modpathol.2016.188">Peritoneal mesothelioma</a>, a rare cancer of the abdomen</li>
<li><a href="http://dx.doi.org/10.1007/s00401-016-1616-3">Anaplastic pleomorphic xanthoastrocytoma (PXA)</a>, a rare neurological tumor</li>
</ul>
<br />
<ul>
</ul>
<h2 id="targeted-sequencing-approach">
Targeted sequencing </h2>
The UCSF500 assay is a targeted panel of approximately <a href="http://cancer.ucsf.edu/intranet/ccgl">500 genes</a> and other genomic regions relevant to cancer diagnosis, prognosis and treatment.<br />
<br />
The sequenced tissue samples are typically a solid tumor biopsy and a matched normal sample, either blood draw or buccal swab. The matched normal is optional, and hematologic (blood) tumors such as leukemia and lymphoma can also be sequenced. In cases of tumor recurrence a previously sequenced normal sample's results will be reused in the new analysis to save time and cost.<br />
<br />
CCGL staff perform DNA extraction, library preparation and hybridization on-site at UCSF. The custom target panel consists of:<br />
<ul>
<li>Exonic regions of about 500 (initially 510, now 480) cancer-associated genes;</li>
<li>Selected introns of about 40 genes;</li>
<li>Microsatellite sequences, for detecting microsatellite instability (MSI);</li>
<li>Scattered SNP sites ("CGH probes") to detect loss of heterozygosity and mutation burden -- these probes are more concentrated near genes where copy number status is known to be actionable.</li>
</ul>
The captured DNA is then sequenced on two Illumina HiSeq 2500 systems at the UCSF Genomics Core.<br />
<br />
The target panel and analysis approach originated in the Bastian lab at UCSF. Elements of this approach can be seen in published papers on <a href="http://dx.doi.org/10.1056/NEJMoa1502583">melanoma progression</a>, the genetic drivers of <a href="http://dx.doi.org/10.1038/ng.3382">desmoplastic melanoma</a>, and our copy number caller <a href="http://dx.doi.org/10.1371/journal.pcbi.1004873">CNVkit</a>.<br />
<!-- .. image: cnvkit diagram, a few chroms (I made an image of this once) --><br />
<h2 id="analysis-on-the-cloud">
Analysis on the cloud</h2>
CCGL's custom-built pipeline for variant detection and analysis runs on the DNAnexus platform. Analysis of a typical sequencing run -- 14 patient samples, sequenced to an on-target coverage depth of 400x, plus 2 control samples -- takes about 4.5 hours.<br />
<br />
The pipeline detects:<br />
<ul>
<li>Small/single nucleotide variants (SNV): GATK HaplotypeCaller and UnifiedGenotyper, FreeBayes, Mutect -- combined into a single "unified" VCF for annotation</li>
<li>Structural variants (SV): Pindel, DELLY -- run independently to detect potential gene fusions; small indels from Pindel are also included in the SNV VCF</li>
<li>Copy number variants (CNV): CNVkit</li>
<li>MSI detection: MSIsensor</li>
<li>Various validity checks and quality metrics.</li>
</ul>
The design of the pipeline was inspired by <a href="https://bcbio-nextgen.readthedocs.io/en/latest/">bcbio-nextgen</a>, but does not share code, due to technical constraints. DNAnexus now provides bcbio-nextgen as an app; this is a recent addition and may be a tempting option for new clinical sequencing services.<br />
<br />
<h2 id="reporting">
Reporting to the oncologist</h2>
The completed analysis results are automatically pulled to an on-site server hosting the signout software. This server and software is used by CCGL staff to review and approve each run's final QC metrics, and by CCGL clinical geneticists (mainly UCSF pathologists) to review the results and generate a PDF report to be entered in the patient's medical record and returned to the ordering oncologist.<br />
<br />
The UCSF500 report highlights clinically relevant genomic features:<br />
<ul>
<li>Somatic SNVs annotated as "Pathogenic" or "Likely Pathogenic" in ClinVar</li>
<li>Copy number alterations</li>
<li>Fusion genes</li>
<li>Microsatellite stability/instability status</li>
<li>Pathogenic germline variants relevant to oncology (incidental)</li>
</ul>
The clinical geneticist also writes a free-text interpretation of the genomic features in light of the patient's disease and medical context, including literature references. The report includes a complete table of detected somatic SNVs of unknown significance as an appendix.<br />
<br />
Finally, the clinical geneticist responsible for the report and the ordering oncologist responsible for the patient discuss therapeutic strategies at the next molecular tumor board meeting.<br />
<br />
This is the UCSF500's key advantage, and the reason this service can only be offered by a medical center, not a startup: CCGL clinical geneticists continue to work with oncologists after delivering the final report so that the patient's medical history and treatment plan can be considered together with genetics and other diagnostic tests. This approach allows us to not only identify relevant therapies and clinical trials, but also to reconsider the initial diagnosis, quickly order follow-up lab tests, and consider germline findings that may affect the patient's family members.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com6tag:blogger.com,1999:blog-266234734515043410.post-90818324841879388162014-11-04T01:30:00.001-05:002014-11-04T01:30:23.285-05:00Preview and preprint: CNVkit, copy number detection for targeted sequencingI've posted a <a href="http://biorxiv.org/content/early/2014/10/29/010876">preprint of the CNVkit manuscript</a> on bioRxiv. If you think this software or method might suit your needs, please take a look and let me know what you think of it!<br />
<br />
<h3>
What is CNVkit?</h3>
CNVkit is a software toolkit for detecting and visualizing germline copy number variants and somatic copy number alterations in targeted or whole-exome DNA sequencing data. (<a href="https://github.com/etal/cnvkit">Source code</a> | <a href="http://cnvkit.readthedocs.org/en/latest/">Documentation</a>)<br />
<br />
The method implemented in CNVkit takes advantage of the sparse, nonspecifically captured off-target reads present in hybrid capture sequencing output to supplement on-target read depths. The program also uses a series of normalizations and bias corrections so it can be used with or without a normal-sample copy number reference to accurately call CNVs. The overall resolution and copy ratio values are very close to those obtained with 180K array CGH.<br />
<br />
We have used CNVkit at UCSF to assess clinical samples for several research projects over the past year.<br />
<br />
<h3>
Putting it in your pipeline</h3>
See the <a href="http://cnvkit.readthedocs.org/en/latest/quickstart.html">Quick Start page</a> for basic usage. The software package is modular so, in addition to the simple "batch" calling style, the underlying commands can be run directly to support your workflow.<br />
<br />
I've attempted to make CNVkit compatible with other software and easy to integrate into sequencing analysis pipelines. The following are currently supported or in development:<br />
<ul>
<li><a href="https://github.com/chapmanb/bcbio-nextgen">bcbio-nextgen</a> -- in progress</li>
<li><a href="https://usegalaxy.org/">Galaxy</a> -- a basic wrapper is in the development <a href="https://testtoolshed.g2.bx.psu.edu/">Tool Shed</a></li>
<li><a href="http://compbio.cs.brown.edu/projects/theta/">THetA2</a> -- CNVkit segmentation output can be used directly as input to THetA</li>
<li><a href="http://www.broadinstitute.org/igv/">Integrative Genomics Viewer</a> -- export segments as SEG, then load in IGV to view tracks as a heatmap</li>
<li>BioDiscovery Nexus Copy Number -- export files to the Nexus "basic" format</li>
<li><a href="http://jtreeview.sourceforge.net/">Java TreeView</a> -- export CDT or .jtv tabular files, then load in JTV for a microarray-like viewing experience</li>
</ul>
If you would like to see CNVkit play nicely with another existing program, and/or support another standard output format, or just want some help getting set up, please let me know on <a href="http://seqanswers.com/forums/showthread.php?t=47910">SeqAnswers</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com3tag:blogger.com,1999:blog-266234734515043410.post-24448485091232143532014-10-03T02:21:00.001-04:002014-10-03T02:51:50.526-04:00On the awesomeness of the BOSC/OpenBio Codefest 2014This summer I was in Boston for a bundle of conferences: <a href="http://www.iscb.org/ismb2014">Intelligent Systems in Molecular Biology (ISMB)</a>, the <a href="http://www.open-bio.org/wiki/BOSC_2014">Bioinformatics Open Source Conference (BOSC)</a> before that, and a very special <a href="http://www.open-bio.org/wiki/Codefest_2014">Open Bioinformatics Codefest</a> before all of it.<br />
<br />
The Codefest was novel, so I'm writing about the highlights here.<br />
<br />
<h3>
A Galaxy Tool for CNVkit</h3>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://wiki.galaxyproject.org/Images/Logos?action=AttachFile&do=get&target=ToolShed.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="87" src="https://wiki.galaxyproject.org/Images/Logos?action=AttachFile&do=get&target=ToolShed.jpg" width="200" /></a></div>
I spent a good portion of the last year developing <a href="https://github.com/etal/cnvkit">CNVkit</a>, a software package for calling copy number variants from targeted DNA sequencing. At the Codefest I wanted to work on making CNVkit compatible with existing bioinformatics pipelines and workflow managements systems, in particular Galaxy and bcbio-nextgen. (I had no prior development experience with either platform.)<br />
<br />
Galaxy is a popular open-source framework for building reusable/repeatable bioinformatic workflows through a web browser interface. In particular, existing software can be wrapped for the Galaxy platform and distributed through the <a href="https://toolshed.g2.bx.psu.edu/%5D">Galaxy Tool Shed</a>. With help from <a href="http://blastedbio.blogspot.com/">Peter Cock</a>, <a href="http://mattshirley.com/about">Matt Shirley</a> and other members of the Galaxy team, I managed to build and successfully run a Galaxy Tool wrapping CNVkit. It's currently visible in the <a href="https://testtoolshed.g2.bx.psu.edu/view/etal/cnvkit">Test Tool Shed</a> and in the main CNVkit source code repository on GitHub. I still need to finalize the Tool, and sometime after that it will hopefully be accepted into the main Tool Shed, making it easily available to all Galaxy users.<br />
<br />
<h3>
bcbio-nextgen</h3>
<a href="http://bcbio.wordpress.com/">Brad Chapman</a>, in addition to his involvement in developing Biopython and Cloud BioLinux and organizing the Codefest itself, is currently leading the development of <a href="http://bcbio-nextgen.readthedocs.org/">bcbio-nextgen</a>, a framework to implement and evaluate best-practice pipelines for high-throughput sequencing analyses. Recent work on this project considered <a href="http://bcbio.wordpress.com/2014/08/12/validated-whole-genome-structural-variation-detection-using-multiple-callers/">structural variants</a>; next steps will consider cancer samples and targeted or whole-exome sequencing, where CNVkit could be a useful component of the analysis pipeline.<br />
<br />
I didn't produce any code for bcbio-nextgen at the Codefest, but I did get a chance to talk to Brad about it a little, and work is now<a href="https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/structural/cnvkit.py"> progressing</a>. A goal of the bcbio-nextgen project is to produce a pipeline that not only works, but works as well as possible. To achieve this, we'll need to develop good benchmarks for evaluating structural variant and copy number variant calls on cancer samples, something of an open problem at the moment.<br />
<br />
<h3>
Arvados and Curoverse</h3>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://arvados.org/images/arvados-logo-scaled-60px.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://arvados.org/images/arvados-logo-scaled-60px.png" /></a></div>
<a href="https://arvados.org/">Arvados</a> is a robust, open-source platform for conducting large-scale genomic analyses. The project originated in George Church's group at Harvard and the Personal Genome Project. <a href="https://curoverse.com/">Curoverse</a> is a startup that has built a user-friendly workflow management system (similar to DNAnexus and Appistry, conceptually) on top of Arvados. The Curoverse front-end can be installed and run locally, and jobs can also be seamlessly dispatched to distributed computing services (like the Amazon cloud); some of Galaxy and bcbio-nextgen run on Curoverse already.<br />
<br />
Curoverse kindly sponsored the Codefest, and a few of the Arvados/Curoverse folks were in attendance and <a href="https://arvados.org/blogs/20">shared some of their work</a> (and stickers, and free trial accounts and compute time) with the rest of us. The Codefest was also blessed with Amazon Web Services bucks, which we could use toward running Cloud BioLinux or Curoverse. Anyway, Curoverse looks cool, and worth keeping an eye on.<br />
<br />
<h3>
Biopythoneers of the world unite</h3>
<div class="separator" style="clear: both; text-align: center;">
<a href="http://biopython.org/w/images/logo/biopython_tiny.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="http://biopython.org/w/images/logo/biopython_tiny.png" /></a></div>
The core developers of Biopython are distributed globally, and BOSC is a fairly unique opportunity for any of us to meet in person. The Codefest provided a nice setting for Peter Cock, <a href="http://bow.web.id/">Wibowo "Bow" Arindrarto</a> and I to get together, stake out a table and hack on Biopython for a couple days.<br />
<br />
We started with a survey of the issue tracker and addressed some long-standing bugs. Bow then moved on to explore an idea for splitting the Biopython distribution into smaller, separately installable modules, while I cleaned some dark corners of Bio.Phylo and enabled automatic testing of the code examples in the Phylo chapter of the Biopython Tutorial and Cookbook. Peter worked on his new BGZIP module and SAM/BAM branch in Biopython, and at some point stated that Biopython will have native (pure-Python) SAM/BAM support soon.<br />
<br />
<h3>
The scene</h3>
We met up at <a href="http://www.hackreduce.org/">hack/reduce</a>, a hackspace next to MIT and Kendall Square -- a fairly unassuming low-rise brick building, converted from an industrial space and retrofitted with good Wi-Fi, coffee urns and other essentials.<br />
<br />
The environment inside was friendly and helpful. Note the distinction between "codefest" and "hackathon": This one was collaborative, not competitive, and welcomed both newcomers and veterans of open-source projects. In addition to the Biopythoneers, Galaxy was well represented, with <a href="https://twitter.com/jmchilton">John Chilton</a> and <a href="https://twitter.com/biocrusoe">Michael Crusoe</a> conveniently within hollering distance of Team Biopython. Groups from Arvados/Curoverse, Cloud BioLinux, and individuals who are involved in a variety of other projects were there, too. Some people just came to meet up and network. <a href="http://scholar.google.com/citations?user=nHYZW2MAAAAJ&hl=en">Chalalai Chaihirunkarn</a> from Carnegie Mellon University was there to study the dynamics of the Codefest itself, and she will report on it at some point.<br />
<br />
At BOSC, kicking off the second day of the conference, Brad summarized our accomplishments at the Codefest: <br />
<ul>
<li><a href="http://video.open-bio.org/video/22/codefest-2014-report">Video</a></li>
<li><a href="https://docs.google.com/presentation/d/114yvrK0Veasc_ns_rg484j2xxRi1h7wNlU2XKONuUqY/edit?usp=sharing">Slides</a></li>
<li><a href="https://docs.google.com/document/d/1yADE2bb0rEU6TASxuSPsvTdHvh_rtCXzJrsL3NWzxXE/edit#heading=h.xx1hyiuk93d3">Collected notes</a> from before, during and after the Codefest</li>
</ul>
I recommend attending the next OpenBio Codefest to anyone who is interested in it. Even if you aren't currently involved in an open-source project, BOSC and the Codefest are unique, useful opportunities for personal education and professional development. In any case it's an interesting and fun experience.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com2tag:blogger.com,1999:blog-266234734515043410.post-715102455613889642014-06-28T13:39:00.002-04:002014-07-02T12:30:14.646-04:00Tomorrowland never dies: A clever rapid transit system sees life in Tel AvivA skyTran demo track <a href="http://www.bbc.com/news/technology-27995437">will be installed in Tel Aviv, Israel</a> during the next year, and a larger commercial installation is scheduled for 2016. If this system works well, will put every other city's mass transit options to shame.<br />
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhNpA0paEIGZLlBDS2xwmq2nPdFv9U9gOj6cf0n8mjsM69vS9FzZa_kDg-4NH_KxPAcBuyXpkMbQnE3FWFmOHEecSI4aK04qQ9GUsNh3YaInCwLq2KPS-z_w8G5e11210oiE8EjTV-u3Cwl/s1600/skyTran_VehicleSketch-Open-003.jpg" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhNpA0paEIGZLlBDS2xwmq2nPdFv9U9gOj6cf0n8mjsM69vS9FzZa_kDg-4NH_KxPAcBuyXpkMbQnE3FWFmOHEecSI4aK04qQ9GUsNh3YaInCwLq2KPS-z_w8G5e11210oiE8EjTV-u3Cwl/s1600/skyTran_VehicleSketch-Open-003.jpg" height="225" width="320" /></a>The <a href="http://skytran.us/intro/">skyTran</a> concept should be easy to grasp if you've been to Disneyland, specifically Tomorrowland:<br />
<ol>
<li>Start with the monorail and split it into small autonomous cars seating two people (like the now-defunct <a href="http://www.yesterland.com/peoplemover.html">PeopleMover</a>, if you remember it).</li>
<li>Flip it upside down, hanging from the track rather than on top of it, so less can go wrong (a hint of <a href="http://www.yesterland.com/skyway.html">SkyWay</a>, now).</li>
<li>Modernize the design: aerodynamic shape, a maglev track, and sensors to guide and space the cars and allow them to brake quickly if there's a problem ahead.</li>
</ol>
The cars are shunted off from the main track to allow passengers to get on or off without interrupting the flow of traffic — like a freeway or express train. The tight spacing, automatic routing, and lack of stops along the way allow the system to carry large numbers of people to their destinations much more efficiently than a standard highway or even a train.<br />
<br />
<div style="text-align: center;">
* * * </div>
<br />
I found <a href="http://www.canosoarus.com/">Doug Malewicki</a>'s website for skyTran in late 2005, fresh out of college. The idea looked solid and I was enthusiastic about it, but the website didn't do justice to the engineering team behind it. I volunteered to revamp the website, and met with Doug a couple of times to go over ideas. I put together a simple static site as a demo, and after much exertion, got the CSS to look all right in both Firefox and IE (because that's what people used at the time, grandkids) on a variety of screen sizes and default font sizes. But for our next meeting, Doug printed out the webpages, and — due to some combination of hardware and software that I will never know — the fonts tripled in size and the layout turned to garbage. Doug was displeased, I was helpless. In conclusion, I don't have what it takes to do front-end web development. It was probably this experience more than any other single event that motivated me to go to grad school.<br />
<br />
Anyway, I'm delighted to see that skyTran is still moving forward.<br />
<br />
When I first saw the skyTran design in 2005, smartphones were not popular yet, so instead of using a smartphone app to summon a car and pay for it, passengers would carry a keychain-sized RFID dongle for payment (e.g. FasTrak), and simply queue up at a raised platform to catch a car. The design I saw also indicated a maximum speed of 150 mph, making it suitable for most medium-distance travel in a metro area, and probably competitive with <a href="http://en.wikipedia.org/wiki/California_High-Speed_Rail">California's long-delayed high-speed rail</a>, but not as fast as airlines for long-distance travel.<br />
<br />
Obviously, this would be great in a spread-out US city like Los Angeles or San Jose, but implementing it would be politically impossible. At the time Doug told me about it, he said he was going to pursue privately funded initiatives, and he had a team in Seattle building a proper prototype. I thought a good candidate to try the technology would be a city-state like Singapore, with a strong centralized authority and a keen interest in efficient, scalable civic development. It seems Tel Aviv has the same will and ability to develop new infrastructure. So, is there any good reason California can't do the same?etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-42302521196328137742013-12-24T19:18:00.001-05:002014-01-13T02:34:31.923-05:00The blinders of peer reviewDoes pre-publication peer-review isolate a finding from the field during the process? Sure, and that's partly the point of it, but it can lead to some inconveniences when two related papers from separate groups undergo peer review at the same time.<br />
<br />
Earlier this year I published a <a href="http://www.biomedcentral.com/1471-2148/13/117">bioinformatic analysis of the rhoptry kinases (ROPK)</a>, a lineage-specific family of signaling proteins involved in the invasion mechanisms of <i>Toxoplasma gondii</i>, <i>Eimeria tenella</i> and related eukaryotic parasites. During this study I found four <i>T. gondii</i> proteins (and their orthologs in other species) that have the hallmarks of ROPKs, including a predicted signal peptide, a protein kinase domain more similar to other ROPKs than to any other known kinases, and mRNA expression patterns matching those of other ROPKs. I named these genes numerically starting after the highest-numbered ROPK <a href="http://www.sciencedirect.com/science/article/pii/S193131281000243X">previously published (ROP46)</a>.<br />
<br />
To informally reserve the names ahead of the publication of my own article, I posted notes on the corresponding ToxoDB gene pages: <a href="http://toxodb.org/toxo/showComment.do?projectId=ToxoDB&stableId=TGME49_252500&commentTargetId=gene#45393">ROP47</a>, <a href="http://toxodb.org/toxo/showComment.do?projectId=ToxoDB&stableId=TGME49_234950&commentTargetId=gene#45413">ROP48</a>, <a href="http://toxodb.org/toxo/showComment.do?projectId=ToxoDB&stableId=TGME49_274170&commentTargetId=gene#45423">ROP49</a> and <a href="http://toxodb.org/toxo/showComment.do?projectId=ToxoDB&stableId=TGME49_249470&commentTargetId=gene#45433">ROP50</a>. My professor and I made some inquiries with other <i>T. gondii</i> researchers to see if it would be possible to confirm the localization of these proteins to the rhoptry organelle, in order to solidify our argument. Without a peer-reviewed publication to point to, though, this seemed to be the most we could do to promote the new gene names.<br />
<br />
In parallel, another well-regarded lab that specializes in <i>T. gondii</i> rhoptry proteins, including but not limited to ROPKs, investigated the <a href="http://www.sciencedirect.com/science/article/pii/S0020751913002221">localization and function of three other proteins</a> that mRNA expression had indicated were associated with other rhoptry proteins. It's great work. However, their paper and ours both passed through peer review at roughly the same time (earlier this year); we both followed the same numerical naming scheme for rhoptry proteins, starting after ROP46; and unfortunately, we ended up assigning the names ROP47 and ROP48 to different <i>T. gondii</i> proteins.<br />
<br />
Crud.<br />
<br />
How could this confusing situation have been avoided? EuPathDB is widely used, but it's not the primary source for gene names and accessions, and a user-submitted comment alone has fairly limited visibility. I presented a poster at the 2012 Molecular Parasitology Meeting, where many of the active Toxo enthusiasts gather each year, but the choice of new gene names was a minor detail on the poster. Heck, I even had breakfast with the other group's PI, but we only talked about <a href="http://www.jbc.org/content/288/48/34968.short">curious features of established rhoptry proteins</a>, not the novel ROPs we were each about to propose.<br />
<br />
The only way to really claim a gene name is with a peer-reviewed publication.<br />
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS9M7RfbJEbxiMwFRzSUCfhIznOEnK5ctS3AGgLvQNAlpzdrdoU-sJ-EXcvCNOcJnJB2OprCqV7KhGVHfBiS57cZUJ4dH_tfjE0LwnAsxf9sS8Eo8TfFqs9YfVczig882hnGgZPkCxrULq/s1600/peer-review-takes-too-damn-long.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="239" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS9M7RfbJEbxiMwFRzSUCfhIznOEnK5ctS3AGgLvQNAlpzdrdoU-sJ-EXcvCNOcJnJB2OprCqV7KhGVHfBiS57cZUJ4dH_tfjE0LwnAsxf9sS8Eo8TfFqs9YfVczig882hnGgZPkCxrULq/s320/peer-review-takes-too-damn-long.png" width="320" /></a>
<br />
<center>
* * *</center>
<br />
Until now I didn't really grasp the importance of public preprint servers like <a href="http://arxiv.org/">arXiv</a>, <a href="http://biorxiv.org/">BioRxiv</a> and <a href="https://peerj.com/preprints/">PeerJ PrePrints</a> — at least in the life sciences where a good article can be published outside a glamor mag within a few months. (In physics and mathematics, peer review and publication typically take much longer.) It was hard enough to get people I knew to review my articles before submitting them to a journal; would anyone really leave useful comments out of the blue if I posted an unreviewed paper on a preprint server? Answer: Maybe, but there's more to preprints than that.<br />
<br />
"Competitors" have their own projects, usually planned around their own grants. They <i>could</i> drop everything and copy your idea if they saw it. More likely, they will do the same thing they'll do when they see your final published paper, which is to take this new information into account as they pursue their own projects. You do want to make an impact on the field, don't you?<br />
<br />
Pre-publication peer-review is a well-established system for gathering detailed suggestions from uninvolved colleagues, a useful stick to force authors to improve their manuscripts, and sometimes a filter for junk. F1000 has an innovative process of publishing submissions first after a cursory screening, then collecting peer reviews and allowing authors to revise the manuscript at their leisure, apparently. Once a manuscript has been reviewed, revised and approved, it receives a tag indicating that it has been properly peer-reviewed. PeerJ takes a more conservative approach, hosting a preprint server alongside but separate from their peer-reviewed articles. Are either of these the way forward?<br />
<br />
F1000 is new on the scene, and it may be too soon to tell if this is going to be a success. For one thing, will authors be motivated enough to correct their manuscripts promptly? PLoS One once fought a mighty battle against the perception that they weren't peer-reviewed. That stigma came out of thin air, and has been overcome — but will F1000 have to fight the same battle again, since their articles really are a mix of various states of peer-review? I hope not, because many scientists could benefit from having a few holes poked in the wall of pre-publication peer review.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-55990948551552434782013-10-08T02:26:00.000-04:002013-10-09T00:08:08.803-04:00Old Cajun wisdom for the young scientist<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEibCh1Y724-M9fQkju4OwumhWwMqcA6dXQDeGa8v_AjKXyKPl70F0qUe5QITMpdMCgwcQxOnkRQOAIRQ2Vylg-IzkKKKzRXz2RIm_ph51XlrrJg_INBpW_bgx52cMWe5O-_3lcxT3Pe6OZP/s1600/BEL+-+Cajun+Wisdom.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="153" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEibCh1Y724-M9fQkju4OwumhWwMqcA6dXQDeGa8v_AjKXyKPl70F0qUe5QITMpdMCgwcQxOnkRQOAIRQ2Vylg-IzkKKKzRXz2RIm_ph51XlrrJg_INBpW_bgx52cMWe5O-_3lcxT3Pe6OZP/s200/BEL+-+Cajun+Wisdom.jpg" width="200" /></a>During the summer after I started grad school, Paul Graham posted an interesting article: "<a href="http://www.paulgraham.com/ramenprofitable.html">Ramen Profitable</a>." I thought it was inspiring in two ways:<br />
<ol>
<li>It made the point that a startup founder is in a much safer, more comfortable and more productive position once enough revenue is coming in to support minimal cost-of-living; raising additional funding is no longer the most important thing. Replace "founder" with "grad student"/"postdoc"/"omg it never ends," and "funding" with "funding" (well, it sounds completely different when you change the context), and you've explained why applying for long-shot grants is such a resource sink, yet we all do it anyway, and why a grim but reliable stipend is an acceptable equilibrium for many academics.</li>
<li>The article included a vague but basically great recipe for beans and rice in the footnotes.</li>
</ol>
I cooked variations of this recipe through the rest of grad school, and now as a postdoc. It hits the sweet spots for flavor, cost, ease of prep, and sheer volume of leftovers. Here's the specific recipe I converged on over the years.<br />
<br />
<h3 style="text-align: center;">
Red and Black Beans and Rice (SEC)</h3>
<div style="text-align: center;">
— or —</div>
<h3 style="text-align: center;">
Life of the Mind Beans and Rice (other conferences)</h3>
In a rice cooker, mix together:<br />
<ul>
<li>1 c. dry <b>rice</b> (white, brown or parboiled)</li>
<li>1/2 c. <b>quinoa</b> </li>
<li>2 1/2 to 3 c. <b>water</b> (see rice packaging and past experience)</li>
<li>(optional: 1 Tbsp. <b>white vinegar</b>)</li>
<li>(optional: 1 tsp. <b>garlic powder</b>)</li>
</ul>
<br />
Start the rice cooker and let it do its thing. Heat in a large pan over medium:<br />
<ul>
<li>1-2 Tbsp. <b>vegetable oil</b></li>
</ul>
<br />
Add:<br />
<ul>
<li>1 medium/large or 2 small <b>yellow onion</b>(s), chopped</li>
<li>4-6 cloves <b>garlic</b>, chopped; or 2 Tbsp. garlic paste</li>
</ul>
<br />
Cook until onions are translucent, about 3 minutes. Add:<br />
<ul>
<li>2-4 oz. <b>andouille sausage</b>; or spicy pork sausage; or other spicy sausage; chopped </li>
<li>2 stalks <b>celery</b>, chopped</li>
<li>1 <b>green bell pepper</b>, chopped</li>
<li>1/4 c. <b>okra</b>, chopped</li>
<li>(optional: 1/2 to 1 <b>jalapeno pepper</b>, chopped)</li>
</ul>
<br />
Stir casually for 4-5 minutes. (This is a good time to grab a beer from the fridge.) Season with:<br />
<ul>
<li>1 tsp. <b>black pepper</b></li>
<li>1/2 to 1 Tbsp. <b>paprika</b></li>
<li>1/2 to 1 Tbsp. <b>cumin</b></li>
<li>1 "serving" <b>condensed chicken broth</b>; or 1 cube chicken boullion, crumbled</li>
<li>(optional: 1/4 to 1/2 tsp. <b>red/cayenne pepper</b>) </li>
</ul>
<br />
Stir for 1 minute to mix thoroughly. Open:<br />
<ul>
<li>1 14-oz can <b>black beans</b></li>
<li>1 14-oz can <b>red beans</b>; or kidney beans; or more black beans</li>
</ul>
<br />
Pour in the liquid and about 1/3 to 1/2 of the beans from each can into the pan, and stir to mix. Put the rest of the beans in a small bowl or sturdy (Pyrex) beaker and mash somewhat. Add the semi-mashed beans to the pan. Stir for another 3-5 minutes to let the stew thicken.<br />
<br />
(When the rice cooker finishes, take off the glass lid and drape a cheesecloth or thin towel over it for a few minutes while the rice cools a little.)<br />
<br />
Serve the beans over the rice/quinoa blend in a shallow bowl.<br />
<br />
Leftovers are even better.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-457761466627825072013-06-09T21:25:00.001-04:002013-06-13T00:31:15.665-04:00Homebrew chroniclesI've started up another blog on home brewing. The first post covers a <a href="http://goodwolfbrew.blogspot.com/2013/06/grains-of-fury-first-brew.html">Pumpkin Spice Porter</a> and some important lessons learned about water.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-61710330993821346982013-01-19T16:56:00.003-05:002014-10-03T02:25:16.757-04:00Overly honest methods / cheap tricks for multiple sequence alignmentWhat's the quickest way to get a so-so sequence alignment? Once you have this,
there are plenty of efficient ways to refine the alignment, infer a tree, build
a profile to search for related sequences, or any combination of these at once.
<a href="http://www.biomedcentral.com/1471-2105/5/113">MUSCLE</a> and <a href="http://mafft.cbrc.jp/alignment/software/algorithms/algorithms.html">MAFFT</a> spend a surprising amount of time computing initial pairwise
distances and computing a guide tree; building a multiple alignment from this
point is quite efficient.
As part of my recent fascination with absurd algorithms and other unpublishable
methods, I looked into this a bit.<br />
<br />
<h3>
Pairwise distances</h3>
Rigorously/naively, we could compute pairwise alignments with a dynamic
programming algorithm like Smith-Waterman and take the pairwise score as the
measure of similarity. Running all-versus-all BLAST and taking the best HSP score
(similarity) or e-value (distance) would give roughly the same result. But we
can do a crappier job of this if we put some thought into it.<br />
<br />
For each sequence, let's count the number of occurrences of each k-mer (k to be
determined). The similarity between sequences is then the sum of absolute
differences in k-mer counts between two sequences.
<br />
<pre class="prettyprint lang-py"># Count distinct tripeptides in each sequence
from collections import Counter
from Bio import SeqIO
K = 3
def get_kmers(seqstr, k):
for i in range(len(seqstr) - k):
yield seqstr[i:i+k]
def kmer_dist(kmers1, kmers2):
"""Sum of absolute differences between k-mer counts."""
keys = set(kmers1.iterkeys()).union(set(kmers2.iterkeys()))
return sum(abs(kmers1.get(key, 0) - kmers2.get(key, 0))
for key in keys)
records = list(SeqIO.parse(fname, 'fasta'))
lookup = dict((rec, Counter(get_kmers(str(rec.seq), K)))
for rec in records)
# Example:
rec1, rec2 = records[:2]
a_distance = kmer_distance(lookup[rec1], lookup[rec2])
</pre>
If some of the sequences have repetitive regions — think of paralogs from
closely related species, or within a species — this approach will
over-emphasize the content and size/number of those repeats at the expense of
less common but more meaningful k-mers. We could down-weight abundant k-mers by
taking the square root of each k-mer count before computing pairwise
differences.
<br />
<pre class="prettyprint lang-py">from math import sqrt
def kmer_dist(kmers1, kmers2):
keys = set(kmers1.iterkeys()).union(set(kmers2.iterkeys()))
return sum(abs(sqrt(kmers1.get(key, 0)) - sqrt(kmers2.get(key, 0)))
for key in keys)
</pre>
But why stop at the square root? <i>Reducto ad awesome</i>, let's drop the k-mer counts and simply track whether the k-mer is present in a sequence or not.
<br />
<pre class="prettyprint lang-py">def kmer_dist(kmers1, kmers2):
"""Proportion of non-overlapping elements between two sets."""
return (len(kmers1.symmetric_difference(kmers2))
/ float(len(s1.intersection(s2)))
lookup = dict((rec, set(get_kmers(str(rec.seq), K)))
for rec in records)
</pre>
This also opens up some bit-twiddling opportunities which I have not pursued. MUSCLE, MAFFT and ClustalW all have options to use this trick to quickly estimate pairwise distances for large sequence sets.<br />
<br />
<h3>
Non-phylogenetic guide tree</h3>
The guide tree determines how the pairs of sequences will be combined.
The tree is typically inferred by neighbor-joining based on the computed
pairwise distances, an approach that's O(N<sup>3</sup>) using the standard algorithm (though rumor has it that quadratic time is possible).
Pairs of sequences are aligned according to the tips of the tree, then
sub-alignments are treated as internal nodes in the tree aligned to each other.
(<a href="http://code.google.com/p/prank-msa/wiki/PRANK">PRANK</a> also uses the guide tree to infer indel states for characters, but that
falls into the category of "good" algorithms which we're not interested in
here.)<br />
<br />
Let's set aside the evolutionary meaning of the alignment for the moment (<i>gasp</i>)
and look for a fast way to determine a reasonable order for performing pairwise
alignments of individual sequences and sub-alignments.
(MUSCLE and MAFFT choose UPGMA over more "robust" NJ methods, despite or perhaps even
because of long-branch attraction, with the same rationale.)
Treat the collection of pairwise distances between sequences as a graph, and
use Kruskal's algorithm to compute the minimum spanning tree in linear time.
<br />
<pre class="prettyprint lang-py">import networkx
G = networkx.Graph()
for i, rec1 in enumerate(records):
for rec2 in records[i+1:len(records)]:
dist = kmer_distance(lookup[rec1], lookup[rec2])
G.add_edge(rec1.id, rec2.id, weight=dist)
H = networkx.minimum_spanning_tree(G)
# Check: Are the sequences grouped together reasonably?
import pylab
nx.draw_graphviz(H, prog='twopi', node_size=0)
pylab.show()
</pre>
The edges of the resulting graph give the order of alignment; start with nodes
of degree 1 as individual sequences to align with their only neighbor, then
treat nodes with degree >1 as sub-alignments to align to each other,
prioritizing shorter branches.<br />
<br />
<i>Application:</i> In the <a href="https://github.com/etal/fammer/blob/master/fammer/tmalign.py">tmalign module</a> of Fammer, I combine pairwise
alignments of 3D protein structures, obtained with <a href="http://zhanglab.ccmb.med.umich.edu/TM-align/">TMalign</a> in a subprocess, using a minimum spanning tree
where the edge weights are the reciprocal of the TM-scores
etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-75770760328207497832012-11-08T22:19:00.004-05:002012-11-08T22:31:45.619-05:00Journal article: Biopython's Bio.PhyloWe're not known for the timeliest reporting here at etalog, but FYI, our article on the Bio.Phylo module for phylogenetics in Biopython is now up in its final form on the BMC Bioinformatics:<br />
<br />
<a href="http://www.biomedcentral.com/1471-2105/13/209/" target="_blank">Bio.Phylo: A unified toolkit for processing, analyzing and visualizing phylogenetic trees in Biopython</a><br />
<br />
It's a quick read. Here's one of the figures we cut from the manuscript which shows some of the key features of the module:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgT9bxenovZqkMp6qZ_qAtFrBex1cR150GblrxWmfb5zwogGjqcfsgyhQk5JzPjHSVMfqZO4B5tPmIh-dzzXuV2xwmrV3qqTfiinvNcwulGsUerWF0fQ6RNrLSdezaSc83ZvtDgtbiNbx75/s1600/fig2-pygments.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgT9bxenovZqkMp6qZ_qAtFrBex1cR150GblrxWmfb5zwogGjqcfsgyhQk5JzPjHSVMfqZO4B5tPmIh-dzzXuV2xwmrV3qqTfiinvNcwulGsUerWF0fQ6RNrLSdezaSc83ZvtDgtbiNbx75/s1600/fig2-pygments.png" /></a></div>
<br />
This manuscript gave us a chance to write about a few things we haven't had a good reason to write about elsewhere -- the design rationale, nice figures, performance, and a couple of real-world use cases.<br />
<br />
To get started, though, I recommend the main documentation:<br />
<ul>
<li><a href="http://biopython.org/DIST/docs/tutorial/Tutorial.html">Biopython Tutorial</a> (Phylo is currently chapter 12)</li>
<li><a href="http://biopython.org/wiki/Phylo">Phylo wiki page</a></li>
<li><a href="http://biopython.org/wiki/Phylo_cookbook">Wiki cookbook</a></li>
</ul>
<h3>
History / a bit of navel-gazing</h3>
This project stemmed from a Google Summer of Code <a href="http://informatics.nescent.org/wiki/PhyloSoC:Biopython_support_for_parsing_and_writing_phyloXML" target="_blank">project</a> in 2009 to implement a phyloXML parser for Biopython, mentored by <a href="http://cmzmasek.blogspot.com/" target="_blank">Christian Zmasek</a> and <a href="http://bcbio.wordpress.com/" target="_blank">Brad Chapman</a>. This was administered by the National Evolutionary Synthesis Center (NESCent); the Open Bioinformatics Foundation (OBF) didn't start administering its own GSoC projects until the following year. It was a fun summer, and I got to become more involved in Biopython as a result.<br />
<br />
After GSoC ended, we decided that rather than plug the phyloXML module into Biopython as-is, we could do something akin to SeqIO and AlignIO -- wrap the format-specific parsers (NEXUS and Newick were already supported under Bio.Nexus) under a common interface, and share the core objects. At first we planned to create "TreeIO" and "Tree" modules like BioPerl, but as this could lead to confusion with other types of trees in bioinformatics (e.g. from clustering or other low-level algorithms), we settled on "Phylo", with due credit to <a href="http://rutgervos.blogspot.com/">Rutger Vos</a>'s Perl package Bio::Phylo.<br />
<br />
The next summer I gave a <a href="http://etalog.blogspot.com/2010/10/biophylo-unified-phylogenetics-toolkit.html">talk at BOSC 2010</a> in Boston about this work. My professor was a bit weary of this open-source stuff by that point, so the travel award really helped. (And Airbnb. Boston is not cheap.) The rest is pretty well covered in the <a href="http://etalog.blogspot.com/2012/10/biopython-project-update-at-bosc-2012.html">BOSC 2012</a> talk -- the hack continued, deftly shepherded by <a href="http://blastedbio.blogspot.com/">Peter Cock</a>; <a href="http://brandoninvergo.com/">Brandon Invergo</a> arrived bearing gifts of pypaml, and we managed to get the module into a fairly stable state in between bouts of "real research", enough to write it up.<br />
<br />etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-37277997795493231302012-10-27T14:42:00.001-04:002012-10-27T14:42:33.902-04:00Biopython project update at BOSC 2012Not so hot off the presses, here are the slides from the talk I gave this summer at the Bioinformatics Open Source Conference (<a href="http://www.open-bio.org/wiki/BOSC_2012" target="_blank">BOSC</a>), a satellite conference of Intelligent Systems for Molecular Biology (<a href="http://www.iscb.org/ismb2012" target="_blank">ISMB</a>). Since <a href="http://blastedbio.blogspot.com/" target="_blank">Peter Cock</a> wasn't able to make it out to California this year, he suggested I fill in.<br />
<br />
In addition to the usual coverage of new features, a big theme this year was the recurring successes we've had bring in new core developers via <a href="http://biopython.org/wiki/Google_Summer_of_Code" target="_blank">Google Summer of Code</a>.<br />
<br />
<iframe allowfullscreen="allowfullscreen" frameborder="0" height="356" marginheight="0" marginwidth="0" scrolling="no" src="http://www.slideshare.net/slideshow/embed_code/14070617" style="border-width: 1px 1px 0; border: 1px solid #CCC; margin-bottom: 5px;" width="427"> </iframe> <br />
<div style="margin-bottom: 5px;">
<b> <a href="http://www.slideshare.net/etalevich/biopython-project-update-bosc-2012" target="_blank" title="Biopython Project Update (BOSC 2012)">Biopython Project Update (BOSC 2012)</a> </b> from <b><a href="http://www.slideshare.net/etalevich" target="_blank">Eric Talevich</a></b> </div>
<br />Jan Aerts has also posted <a href="http://www.slideshare.net/jandot/tag/bosc2012">the rest of the BOSC 2012 slides</a>.
etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-40459155037453186062012-08-15T17:39:00.001-04:002012-10-27T14:43:50.682-04:00Code Harvest: The RefactoringI've been hacking on bioinformatics code for four years now, but until now the only work I've really made available to "the community" is in Biopython, mainly <a href="http://www.biomedcentral.com/1471-2105/13/209/" target="_blank">Bio.Phylo</a>.<br />
<br />
The code I write in the lab is under one big Mercurial repository called <span style="font-family: "Courier New",Courier,monospace;">esgb</span>; there's a shell script to install everything, including a bunch of scripts, sub-projects and a sprawling Python library called <span style="font-family: "Courier New",Courier,monospace;">esbglib</span>. Most of my Python programs depend on some functionality in <span style="font-family: "Courier New",Courier,monospace;">esbglib</span>, and usually Biopython and sometimes SciPy as well.<br />
<br />
Having signed the <a href="http://sciencecodemanifesto.org/">Science Code Manifesto</a>, duty calls for me to bundle some of the programs I've written with the next paper I'm working on, and so I've begun a mighty refactoring of <span style="font-family: "Courier New",Courier,monospace;">esbglib</span> to extract the general-purpose, reusable components into Python packages. At the moment it looks like I'll end up with two: <span style="font-family: "Courier New",Courier,monospace;">biofrills</span> and <span style="font-family: "Courier New",Courier,monospace;">biocma</span>.<br />
<br />
<a name='more'></a><br /><br />
<h3>
BioFrills</h3>
This package contains general-purpose sequence analysis functions that I can't merge into Biopython, for one reason or another:<br />
<ul>
<li>Research-grade approaches (e.g. heuristic solutions for problems with no efficient exact solution, like estimating the effective number of independent sequences in an alignment)</li>
<li>Cython versions of some functions for speed-up (currently just one, summing BLOSUM62 scores for a pairwise alignment)</li>
<li>Python 2.6 or 2.7 only; includes recent PyPy</li>
</ul>
I gave the library a silly name to keep it from being somehow mistaken for a Biopython replacement, or even a well-maintained, well-thought-out supplementary library. But in order to make the packages that depend on it easy_installable, I guess I'll need to upload it to PyPI at some point.<br />
<br />
One function that I had in <span style="font-family: "Courier New",Courier,monospace;">esbglib</span> has already been independently invented and merged into Biopython as <span style="font-family: "Courier New",Courier,monospace;">Bio.File.as_handle</span> (Mine was <span style="font-family: "Courier New",Courier,monospace;">esbglib.sugar.maybe_open</span>).<br />
<br />
<a href="https://github.com/etal/biofrills">Check out BioFrills on GitHub</a>.<br />
<br />
<h3>
BioCMA</h3>
For working with MAPGAPS, CHAIN and mcBPPS alignments. These programs use an undocumented, bespoke multiple alignment format called CMA, which I assume means either "Consensus Multiple Alignment" or "CHAIN Multiple Alignment". There are a few tools by the same author for processing this format, but again, we have neither documentation nor source code for them. For years I've felt uneasy about relying on these quirky binaries to manage our most important sequence data... so, I did something about it. Ta-da.<br />
<br />
<a href="https://github.com/etal/biocma">Check out BioCMA on GitHub</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-75731146049696571832012-08-01T12:47:00.000-04:002012-10-27T14:45:20.598-04:00The well-organized data science project<br />
Someone recently asked me about the basic setup a computational scientist needs to conduct research efficiently. I'm pretty satisfied with my current arrangement, which was inspired by this: "<a href="http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1000424">A Quick Guide to Organizing Computational Biology Projects</a>"<br />
<br />
My work is organized into individual "projects" which are each supposed to become papers at some point. I keep each project in <a href="https://www.dropbox.com/">Dropbox</a> to ensure
everything is synced and backed up remotely all the time -- no file left behind. I also use <a href="https://www.mendeley.com/">Mendeley</a>, with a folder for each project's references. Mendeley can generate a project-specific BibTex file from a folder.<br />
<br />
A well-organized project might look like this:
<br />
<a name='more'></a><br />
<ul>
<li>projects/<i>$project_name</i>/</li>
<ul>
<li><b>README</b> -- notes on current progress, observations, to-do items, vague ideas for the future -- a starting point for "what was I doing here?" A very sloppy lab notebook by itself, but combined with the Mercurial revision history, it captures the essentials.</li>
<li>
<b>data</b>/</li>
<ul>
<li>Files sent from collaborators, downloaded and processed data sets. Big databases are kept on the lab's file server, not here.</li>
</ul>
<li><b>src</b>/</li>
<ul>
<li>The
main code base, if there's a core software component to the project.
Throwaway scripts will also be littered throughout the work/ directory
tree.<br />
</li>
</ul>
<li><b>work</b>/</li>
<ul>
<li>For each new "idea", create a new
directory under here, then hack away, freely generating intermediate
files. Copying and modifying files is common; only a small portion of
this is eventually shared or used in the manuscript.</li>
<li>If things look good and the procedure is worth remembering (i.e. redoing with updated data later), organize the idea-related scripts and Bash history into a Makefile.<br />
</li>
</ul>
<li><b>results</b>/</li>
<ul>
<li>Interesting outputs, preliminary
figures and tables that have been or will be sent to collaborators or
used in a presentation (e.g. lab weekly update meetings)</li>
</ul>
<li><b>manuscript</b>/</li>
<ul>
<li>Top level: the manuscript in progress -- in my case usually
LaTeX and BibTex files and a Makefile; also potentially .doc files sent by collaborators, etc. Eventually, also a cover letter for the first submission.<br />
</li>
<li><b>figures</b>/</li>
<ul>
<li>The "final" images that will be used in
the manuscript. These may be stitched together from several components,
e.g. several plots or diagrams generated from programs.</li>
</ul>
<li><b>response</b>/</li>
<ul>
<li>Reviewer comments and the response(s) to reviewers being drafted</li>
</ul>
<li><b>proofs</b>/</li>
<ul>
<li>PDF copies of the manuscript exactly as it was submitted to the journal, at each stage of the publication process.<br />
</li>
</ul>
</ul>
</ul>
</ul>
<br />
I also start a <a href="http://mercurial.selenic.com/">Mercurial</a> repo at the base of each project
for local snapshotting and logging. For the manuscript, it turns out that diffs of LaTeX documents
are often not that helpful after an editing frenzy because of automatic line wrapping, but it's better than nothing. (Git's index is too much of a hassle for this
kind of work -- the benefit of a "clean" patch series is outweighed by the distraction of explicitly adding files before each commit. If I didn't already have Dropbox to sync everything, Bzr's bound repos would be a decent alternative.)<br />
<br />
<h3>
Zen</h3>
<br />
Academic research is not like engineering, where a team sets their sights on a certain goal (perhaps refined over time), divides up the component tasks, and marches onward. We don't know that a project will "work" (in science: improve on previous knowledge or methods) even if we execute the idea the "right" way.<br />
<br />
It's more like investigative journalism, where you latch onto a potential story, follow every promising lead, and seek out multiple lines of evidence to support your claims.<br />
<br />
Because of this, I've adopted a development "approach" that could be called Results-Driven Development, or maybe Failure-Driven Development. The idea is that I need fast results to prove the value in pursuing a line of inquiry, i.e. whatever is going on in a subdirectory of <span style="font-family: "Courier New",Courier,monospace;">work/</span>. It's only worthwhile to invest in proper engineering practices once it looks like there's a "lead". This is different from typical software engineering projects because much of the time in research the first steps in a new direction quickly lead to a rethinking of the problem, and it's necessary to drop the current work and take a different approach immediately, whereas in software engineering you typically already know what you're trying to achieve from the beginning. In startups, "fail fast" means weeks or months; in research it can mean minutes or hours.<br />
<br />
The takeaway is that technical debt is OK during the first few hours of pursuing an idea. If it turns out I've hacked up something reusable, I copy the important bits into my main utility library or script collection and manage them in a reasonable way from then on. But normally the result is a slightly gnarly Bash history and a small set of one-off Bash and Python scripts that are completely specific to that inquiry, and should just stay there.<br />
<br />
Also, note that my style of research might be completely different from yours. My lab is fairly small, and each project starts off solo or with a few data sets from a remote collaborator. If you're providing bioinformatics services for a group driven by high-throughput experiments, your infrastructure needs are much different from mine. (But I still recommend Mendeley and version control.)<br />
<br />
<h3>
Further reading</h3>
<br />
Michael Barton has written an excellent series called <a href="http://www.bioinformaticszen.com/post/decomplected-workflows-introduction/">Decomplected Workflows</a> on this topic. My favorite so far is on <a href="http://www.bioinformaticszen.com/post/decomplected-workflows-makefiles/">Makefiles</a>, where he shows how a realistic workflow would be managed in Make, something I've wrestled with in the past. I'll note that where he advises switching from Rake to Make, there are two changes happening: (1) switching from Rake to the technically very much inferior but more widespread Make, and (2) writing the workflow logic in standalone scripts rather than Ruby functions embedded in the Rakefile. If you only use your Rakefile to control the workflow and instead call out to standalone scripts from Rake (which is very easy), then I think Rake is still the better choice.<br />
<br />
<br />etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com5tag:blogger.com,1999:blog-266234734515043410.post-8289814676607014212012-01-16T17:26:00.000-05:002012-10-27T14:44:59.825-04:00Building an analysis: How to avoid repeating intermediate tasks in a computational pipelineIn my projects, I tend to start with a simple analysis of a limited dataset,
then incrementally expand on it with more data and deeper analyses. This means
each time I update the data (e.g. add another species' protein sequences) or
add another step to the analysis pipeline, everything must be re-run -- but
only a small part of the pipeline actually needs to be re-run.<br />
<br />
This is a common problem in bioinformatics:<br />
<a class="reference external" href="http://biostar.stackexchange.com/questions/79/how-to-organize-a-pipeline-of-small-scripts-together">http://biostar.stackexchange.com/questions/79/how-to-organize-a-pipeline-of-small-scripts-together</a><br />
<br />
How can we automate a pipeline like this, without running it all from scratch
each time? This is the same problem faced when compiling large programs, and
that particular case has been solved fairly well by build tools.<br />
<br />
<a name='more'></a><br /><br />
<h3>
Make</h3>
The <tt>make</tt> utility is normally used to determine which portions of a program
need to be recompiled, and the dependencies between them; it includes special
features for managing C code, but actions are specified in terms of shell
commands. You can just as well use it to specify dependencies between arbitrary
tasks.<br />
<br />
Giovanni Dall'Olio has posted some helpful instructional material:
<a class="reference external" href="http://bioinfoblog.it/?p=29">http://bioinfoblog.it/?p=29</a><br />
<br />
(One of his links is broken -- here's <a href="http://software-carpentry.org/4_0/make/">Software Carpentry's current course on <tt>make</tt></a>.)<br />
<br />
<h3>
Rake, a Really Awesome Make</h3>
If your pipeline is already organized as a pipeline of stand-alone scripts,
Rake is a more-or-less ideal solution:
<a class="reference external" href="http://rake.rubyforge.org/">http://rake.rubyforge.org/</a><br />
<br />
At some point in your project, you'll likely need to do some string processing;
this is where <tt>make</tt> falls down. Ruby happens to be great for this task. Don't
worry if you don't know Ruby -- the basic string methods in Ruby, Python and
Perl are very similar, and regular expressions work roughly the same way.
You can look up everything you need to know on the fly. (I did.)<br />
<br />
Also, Rake's the mini-language for specifying tasks is both concise and
intuitive. There's a built-in facility for maintaining short descriptions of
each step, which is enticing from a "lab notebook" perspective.<br />
<br />
<h3>
How about a Python solution?</h3>
Inspired by Rake, I wrote a small module called tasks.py. Here's the code, plus
a simple script to demonstrate:<br />
<a class="reference external" href="https://gist.github.com/1623117">https://gist.github.com/1623117</a><br />
<br />
This isn't just for the sake of Python evangelism; I have a bunch of
special-purpose modules written in Python (for my lab) that I would like to
use in an elaborate pipeline of my own. It's silly to put each of these
steps in a separate, throwaway Python script that then gets called in a
separate process by Rake. Instead, I can import <a href="https://gist.github.com/1623117">tasks.py</a> into a Python
script and include the dependency-scheduling features in my own programs.<br />
<br />
Key differences from other build systems (e.g. <a href="http://www.scons.org/">SCons</a>, <a href="http://code.google.com/p/waf/">waf</a>, <a href="http://code.google.com/p/ruffus/">ruffus</a>):<br />
<ul class="simple">
<li>This module is not meant to be run from the command-line -- only for
specific parts of your own code</li>
<li>The implementation of each processing step is separate from the
dependency specification (although single-command tasks can still be
defined inline with a lambda expression).
Separate the algorithm from the data, I always say.</li>
<li>Cleanup is specified within the task, not in a separate area of the
Rakefile/Makefile. This makes more sense for a project with heterogenous
processing steps & intermediate files.</li>
</ul>
To be clear, <tt>tasks.py</tt> is not better than Rake or <tt>make</tt> for arranging a set of
scripts that you've already written. I didn't even implement concurrent jobs
(because most of the CPU-intensive steps I use are calls to programs that are
already multicore-aware; though try adding that feature yourself if you'd
like).etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com1tag:blogger.com,1999:blog-266234734515043410.post-46904448972143198832011-11-02T14:42:00.000-04:002012-10-27T14:45:20.596-04:00Journal article: Comparative kinomics of the malaria pathogen and its relativesHot off the presses!
<br />
<a href="http://www.biomedcentral.com/1471-2148/11/321/abstract">Structural and evolutionary divergence of eukaryotic protein kinases in Apicomplexa</a><br />
<br />
It's a thorough paper, so I'll cover the highlights here.<br />
<br />
<h3>
Why we study apicomplexans</h3>
Apicomplexans are a group of related single-celled organisms which are exclusively parasitic. The best-known member is <a href="https://secure.wikimedia.org/wikipedia/en/wiki/Plasmodium_falciparum_biology"><i>Plasmodium falciparum</i></a>, which causes the most virulent form of malaria. Another well-studied species is <i>Toxoplasma gondii</i>, which primarily lives in cats but can infect most mammals.<br />
<br />
It's a hugely diverse group. But overall, we know very little about them.<br />
<br />
Our main motivation for studying apicomplexan proteins is to find what features make them distinct from human proteins, so we can then design drugs to target those features specifically -- the drug will identify and disable the parasite protein without the risk of affecting the host proteins, too. We study protein kinases, in particular, because a number of drugs have already been designed to <a href="http://www.cancerquest.org/kinase-inhibitors.html">inhibit kinases in cancer</a>. The same or similar compounds could be used to treat parasitic diseases, potentially.<br />
<br />
<a name='more'></a><br /><br />
From an evolutionary biologist's perspective, apicomplexans are also interesting to study because they belong to an <a href="http://tolweb.org/Eukaryotes">evolutionary branch</a> that is quite divergent from the animals, plants and fungi more familiar to us. By learning about apicomplexan biology, and comparing to other model organisms, we can learn more about eukaryotic diversity and the origin of eukaryotes.<br />
<br />
<h3>
Another perspective on the tree of life</h3>
Many people, including
scientists, think of evolution as a ladder, with single-celled organisms
at the bottom and humans at the top. Different lineages, like green
plants and fungi, each branch off the ladder at some intermediate point,
but evolution is nonetheless mistakenly thought of as a directed
progression from bacteria to protists to fish to humans.<br />
<br />
That's
wrong. It leads to mistakes, such as considering all protists
(single-celled eukaryotes) to be closely related to each other. But even
within Apicomplexa, the evolutionary distance between <i>Plasmodium falciparum</i> and <i>Toxoplasma gondii</i> is as great as the distance between humans and mosquitoes.<br />
<br />
I'm
particularly proud of Figure 1 in the paper, which includes a species
tree that inverts the traditional view: The closest human relative,
yeast, is at the bottom, and layers of increasingly strange and
unfamiliar protists build up to the <i>Plasmodium</i> genus.<br />
<br />
<h3>
Interesting features of proteins and genomes</h3>
When apicomplexan parasites invade a host, they secrete a mixture of dozens of different proteins into a <a href="http://www.nature.com/nrmicro/journal/v6/n1/fig_tab/nrmicro1800_F2.html">protective vacuole</a> formed from the host cell membrane. We'd expect that some of these proteins are essential for invasion and virulence, and therefore good targets for inhibition or diagnosis.<br />
<br />
Two apicomplexan-specific protein kinase families are known to be exported. The FIKK family appears in 1 copy in most apicomplexans, but is amplified to 21 copies in <i>P. falciparum</i> and 6 copies in <i>P. reichenowi</i>, and does not appear in any species outside the Apicomplexa. Another family, called rhoptry kinases (ROPK) after the apicomplexan organelle they're localized to, appears in dozens of copies in coccidians (<i>T. gondii</i>, <i>Neospora caninum</i>, <i>Eimeria tenella</i>, <i>Sarcocystis neurona</i>), but not in any other lineage of Apicomplexa. <i>Plasmodium</i> and others still contain rhoptries, but there are no kinases in the protein cocktail those rhoptries contain.<br />
<br />
As obligate parasites, apicomplexans evolve under different evolutionary contraints than free-living organisms like yeast and humans. Many genes are no longer necessary, and some may even be a liability if they interact with the host's own biochemical pathways. Because of this, we see widespread gene loss and overall compaction of apicomplexan genomes.<br />
<br />
One especially curious case is the loss of upstream regulators of the MAPK cascade -- a signaling pathway found in almost all eukaryotes, consisting of 3 or 4 protein kinases each activating the next in a sort of biochemical relay. Apicomplexans contain 2 to 3 copies of the downstream protein kinase, MAPK, but the rest of the pathway components (STE7, STE11, STE20) are generally lost, and none of the surveyed apicomplexans had a complete MAPK cascade. So there's an open question: What other proteins take the place of the STEs in this important pathway, or have MAPK-like features? Is there an Achilles heel to be discovered?<br />
<br />
<h3>
The project</h3>
We:<br />
<ol>
<li>Identified and <b>classified</b> the full set of protein kinases in each of the 17 apicomplexan proteomes available </li>
<li>Devised a pipeline to identify <b>apicomplexan-specific ortholog groups</b> in known protein kinase families</li>
<li>Compared these ortholog groups to the typical members of the kinase family to find specific <b>sequence motifs</b> that distinguish the divergent ortholog group</li>
<li>Mapped these motifs onto protein <b>structures</b>; reviewed the literature to understand possible functions and <b>functional differences</b> related to these motifs</li>
</ol>
Read about what we found <a href="http://www.biomedcentral.com/1471-2148/11/321/abstract">here</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-36910218235260451302011-10-28T10:32:00.000-04:002012-10-27T14:45:20.594-04:00Journal article: Our insights into the structure and activation mechanism of ErbB/EGFR protein kinasesHere's an article my lab published in <i>PLoS One</i>:
<br />
<a href="http://dx.plos.org/10.1371/journal.pone.0014310">Co-Conserved Features Associated with cis Regulation of ErbB Tyrosine Kinases</a>
<br />
<br />
I'll give a quick summary of it here. (Don't worry, this isn't a new direction for this blog.)
<br />
<br />
This is a study of the structural mechanisms of a certain protein family, called <a href="http://en.wikipedia.org/wiki/ErbB">ErbB</a> or EGFR (epidermal growth factor receptor), which is frequently involved in cancer. This family belongs to a protein superfamily called <b><a href="http://kinase.com/wiki/index.php/Introduction_to_Kinases">protein kinases</a></b>.<br />
<br />
<a name='more'></a><br /><br />
<h3>
Biochemistry background</h3>
Kinases are enzymes which perform a type of post-translational modification, <b>phosphorylation</b>: The kinase transfers a phosphate group from adenosine triphosphate (<b>ATP</b>) to another <b>substrate</b> molecule, leaving adenosine diphosphate (ADP) and the phosphorylated substrate.<br />
<br />
Protein kinases are kinases that act on protein substrates, i.e. the phosphorylated molecule is another protein. The substrate could even be another protein kinase, so activation of the first protein kinase causes it to phosphorylate and activate another protein kinase, and so on. This is a type of <b>signal transduction</b>.<br />
<br />
Signal transduction is how the cell senses and reacts to its environment, and also its own internal conditions. In the case of ErbB and other receptor tyrosine kinases, the signal starts at the surface of the cell (e.g. epidermal growth factor binds to the extracellular portion of EGFR) and activates the kinase, which then begins sending these phosphorylation signals. These signals are then relayed throughout the cell to trigger other activities, such as <b>cell division</b> or the <b>transcription</b> of certain genes.<br />
<br />
What happens if a protein kinase gets "locked" into the active state, somehow? In the case of EGFR, it's as if the cell thinks it's constantly receiving the growth factor. If this signal isn't blocked by another "gatekeeper" in the cell, then the cell will grow uncontrollably -- and become <b>cancer</b>.<br />
<br />
<h3>
How the enzyme works</h3>
Protein kinases (PKs) all consist of two large lobes connected by a flexible hinge. Between the lobes is a binding pocket for ATP; this molecule binds inside the smaller lobe (N-terminal lobe, or N-lobe). The larger lobe (C-terminal lobe or C-lobe) provides a binding site for another protein, which will be the kinase's substrate.<br />
<br />
The general mechanism of all protein kinases goes like this:<br />
<ol>
<li>The kinase is initially in an <b>inactive</b> state, with the hinge "open" and the two lobes a bit further apart. Since ATP binds in the N-lobe and the substrate binds to the C-lobe, no phosphate is transferred when the two lobes are apart like this.</li>
<li>By some mechanism (it varies between different kinase families), the two lobes are brought closer together, and the kinase becomes <b>active</b>.</li>
<li>ATP binds to the ATP binding pocket, a substrate binds to the C-lobe, some amino acids shift, and a <b>phosphate</b> group is detached from ATP and reattached to a specific amino acid on the substrate.</li>
<li>The ADP and phosphorylated substrate are released.</li>
</ol>
Step 2 is the part we're interested in. How do some recurring, cancer-associated mutations cause EGFR to become "locked" in the active conformation? And, can we reverse it?<br />
<br />
<h3>
How we think ErbB kinases work</h3>
In the ErbB family, it's not just the two lobes of the kinase domain that are involved in activating the enzyme -- the adjacent sections of the protein, outside the kinase domain, are also involved.<br />
<br />
The long C-terminal tail wraps back around the entire kinase domain and associates with the N-lobe, tethered in place by a few residues in the N-lobe and the other N-terminal flanking region (the juxtamembrane segment, between the kinase domain and the cell membrane). The C-tail is placed so that it can influence the movement and relative positioning of the N- and C-lobes, and therefore regulate the activation of the kinase.<br />
<br />
We also examined the locations of two EGFR mutations (S768I and L861Q) that have been previously identified as occurring frequently in cancers, mapping them onto the structure. These mutations appear in locations that would disrupt the switching mechanism we proposed -- breaking necessary interactions, or forming new interactions that shouldn't be there for proper EGFR function.<br />
<br />
If you'd like to know more, read about it <a href="http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0014310">here</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com1tag:blogger.com,1999:blog-266234734515043410.post-22327110591823127792011-07-07T15:28:00.002-04:002014-06-30T23:09:47.629-04:00The statistics of the "Like" buttonThe launch of <b>Google+</b> reminded me of a question I've had about Facebook and YouTube for a while: What happens when you click the "Like" button?<br />
<br />
Facebook isn't so much about sharing <i>content</i> as <i>sharing</i> content. But YouTube and many other sites like it recommend content based on users' response to the content itself, rather than the shape of the surrounding social network. If you're building an application like this yourself, this article is for you.<br />
<br />
Think of a collection of user-generated pages with a "Like" or "+1" button on each. Users can browse pages at will, or arrive at them randomly from an external site, and after viewing a page will either click the "Like" button or do nothing. I'll refer to the number of times viewers do nothing on a page as "Don't care". I'll also assume you have a site-wide "Top Pages" chart that users can view to see the highest-scoring pages and jump to them.<br />
<br />
<h3>
1. Counting "Likes"<br />
</h3>
The very simplest way to score a page is to count the number of times it's been "Liked":<br />
<pre class="literal-block">score = page_likes</pre>
The main problem with this method is inertia. Old "champion" pages accumulate votes over time, and dominate the rankings. New pages don't have a chance to unseat the champs, even temporarily, to gain visibility for themselves. The site appears stagnant.<br />
<br />
<h3>
2. "Likes" versus views<br />
</h3>
To make good recommendations, you want to measure the quality of a page — the chance that a user will like the recommended page themselves:<br />
<pre class="literal-block">score = page_likes / page_views</pre>
In the long run, this is correct — you estimate the probability based on the frequency of "Likes" in previous views. Old "champion" pages will be unseated in the rankings if a newer page earns a better proportion of likes.<br />
<br />
But for new or little-viewed pages, there's an issue of sample size.<br />
<blockquote>
<ul class="simple">
<li>A page where the first view is "Liked" (probably by the creator/uploader) scores 100% and shoots to the top of the rankings. If a few friends all immediately "Like" the page, it becomes difficult to unseat. There's a lot of noise at the top of the rankings.</li>
<li>A page where the first view is not "Liked" scores 0% and sinks to the bottom of the rankings. If you have some mechanism for purging bad content from your site (i.e. deleting low-scoring pages that are likely spam, trolls or just lame), then this makes that task more difficult.</li>
</ul>
</blockquote>
Intuitively, a page with 80 likes out of 100 views is more likely to be good than a page with 4 likes out of 5 views. A page with zero likes out of 100 views is almost certainly junk, but 5 views without any likes may not mean much at all.<br />
<br />
So, your next goal: Make the best possible estimate of a page's likeability based on the first few views and likes, using some prior knowledge. After that, all reasonable methods should converge on the same score (probability of liking). If a meme catches, people will be able to find that page through other means, and your own rankings will be less crucial to its success.<br />
<br />
<h3>
3. Pseudocounts<br />
</h3>
A <a class="reference external" href="http://en.wikipedia.org/wiki/Pseudocount">pseudocount</a> is a prior estimate of the probability of an event. To make an estimate of actual probabilities based on a small number of samples (the problem at the end of Step 2), add the pseudocounts to the actual counts of each event.<br />
<br />
I'll demonstrate.<br />
<br />
The events here are (a) "Like" and (b) "Don't care". I'm going to use <i>b</i> to represent the pseudocount for "Like". For this section, I choose the probabilities:<br />
<pre class="literal-block">like = b = .1
dontcare = (1 - b) = .9</pre>
The two probabilities should sum to 1.<br />
<br />
How do you get these values? Since <i>b</i> represents the probability that a random user will like an arbitrary page, taking the site-wide average of likes versus views is a good choice:<br />
<pre class="literal-block">b = all_likes / all_views</pre>
To use the pseudocounts, add them to the counts in the formula in step 2:<br />
<pre class="literal-block">score = (page_likes + b) / (page_views + 1)</pre>
(Recall: <tt class="docutils literal">views = likes + dontcares</tt>; after adding pseudocounts, <tt class="docutils literal">(likes + b) + (dontcares + 1 - b) = likes + dontcares + 1 = views + 1</tt>.)<br />
<br />
If the database-wide sums of likes and views are large numbers, this won't significantly affect the "average" score, <tt class="docutils literal">all_likes / all_views</tt>. But it smoothes out the initial scoring for new pages.<br />
<br />
<b>Example:</b> Assume 10% of all views result in a "Like" (b = 0.1).<br />
<br />
A single view without a "Like" places the page slightly below the global average, but not too much. (Odds are, 90% of pages will start out this way.) Additional views without a "Like" slowly sink the page score toward 0.<br />
<br />
<table border="1" class="docutils"><colgroup> <col width="34%"></col> <col width="34%"></col> <col width="31%"></col> </colgroup> <thead valign="bottom">
<tr><th class="head">Likes:Views</th> <th class="head">Calculation</th> <th class="head">Percentage</th> </tr>
</thead> <tbody valign="top">
<tr><td>0:1</td> <td>0.1 / 2</td> <td>5.0%</td> </tr>
<tr><td>0:2</td> <td>0.1 / 3</td> <td>3.3%</td> </tr>
<tr><td>0:3</td> <td>0.1 / 4</td> <td>2.5%</td> </tr>
</tbody> </table>
<br />
A single view with a "Like" gives the page a boost, but not to 100%. This can help it gain traction, but probably won't put it in the top rankings (yet). If subsequent views are also liked, the score continues to rise:<br />
<br />
<table border="1" class="docutils"><colgroup> <col width="34%"></col> <col width="34%"></col> <col width="31%"></col> </colgroup> <thead valign="bottom">
<tr><th class="head">Likes:Views</th> <th class="head">Calculation</th> <th class="head">Percentage</th> </tr>
</thead> <tbody valign="top">
<tr><td>1:1</td> <td>1.1 / 2</td> <td>55.0%</td> </tr>
<tr><td>2:2</td> <td>2.1 / 3</td> <td>70.0%</td> </tr>
<tr><td>3:3</td> <td>3.1 / 4</td> <td>77.5%</td> </tr>
<tr><td>4:4</td> <td>4.1 / 5</td> <td>82.0%</td> </tr>
</tbody> </table>
<br />
Away from the extremes (0% or 100% liked), the effect of pseudocounts is less dramatic, and a mix of "Like" and "Don't Care" (viewed without liking) results in a score closer to what you'd see without pseudocounts — just shifted slightly toward the site-wide average. Notice that a page with two "Likes" out of three views (2:3) is scored almost as well as one "Like" and one view (1:1 above).<br />
<br />
<table border="1" class="docutils"><colgroup> <col width="34%"></col> <col width="34%"></col> <col width="31%"></col> </colgroup> <thead valign="bottom">
<tr><th class="head">Likes:Views</th> <th class="head">Calculation</th> <th class="head">Percentage</th> </tr>
</thead> <tbody valign="top">
<tr><td>1:2</td> <td>1.1 / 3</td> <td>36.7%</td> </tr>
<tr><td>1:3</td> <td>1.1 / 4</td> <td>27.5%</td> </tr>
<tr><td>2:3</td> <td>2.1 / 4</td> <td>52.5%</td> </tr>
<tr><td>2:4</td> <td>2.1 / 5</td> <td>42.0%</td> </tr>
</tbody> </table>
<br />
To increase the effect of pseudocounts, you can put a higher weight on the prior by multiplying the pseudocounts by some constant. If the weighting factor is <i>w</i>, then the calculation is:<br />
<pre class="literal-block">score = (page_likes + (b * w)) / (page_views + w)</pre>
Think of this as the number of "imaginary" users you have rating each page before any real users see it. The calculations above use a weight of 1, equivalent to one user giving a fractional score of .1 to every page before it goes live, and you can see the effect of it. Play with it a bit to see how it affects your rankings.<br />
<br />
<b>Update:</b> Evan Miller has a great writeup on a <a href="http://www.evanmiller.org/bayesian-average-ratings.html">Bayesian approach to modelling this problem</a>, if you'd like to go much further down the rabbit hole.<br />
<h3>
4. Statistical significance<br />
</h3>
You now have a score for each page in your database and a top-to-bottom ranking. Where do you draw the line for "recommending" a page?<br />
<h4>
Quantiles: Best and worst<br />
</h4>
Having sorted all the pages by score, take the top 5% as the "best" and bottom 5% as the "worst". Or choose a fixed number, like 25. It's really up to you.<br />
<br />
The "best" ranking is for users, especially new visitors. Depending on your application, some users might also be interested in the "worst" pages — how else would we find gems like "Friday"?<br />
<h4>
Contingency: Is the score meaningful?<br />
</h4>
Another challenge is to determine when a page's score is statistically meaningful — i.e. the difference between a score of 55% based on 1000 views versus a single view. Using pseudocounts addresses this to some extent at the extremes, but it's still possible for pages with low view counts to score highly. You may also want to purge "junk" content with horribly low rankings — but only once it's been given a fair chance.<br />
<br />
With the <tt class="docutils literal">like</tt> and <tt class="docutils literal">dontcare</tt> counts, site-wide and per-page, set up a 2x2 contingency table:<br />
<br />
<table border="1" class="docutils"><colgroup> <col width="32%"></col> <col width="26%"></col> <col width="42%"></col> </colgroup> <thead valign="bottom">
<tr><th class="head"></th> <th class="head">like </th> <th class="head">dontcare </th> </tr>
</thead> <tbody valign="top">
<tr><td>Page</td> <td>A</td> <td>B</td> </tr>
<tr><td>Global</td> <td>C</td> <td>D</td> </tr>
</tbody> </table>
<br />
To evaluate the significance, use a Chi-square test with one degree of freedom (df=1), or if you're picky, <a class="reference external" href="http://en.wikipedia.org/wiki/Fisher%27s_exact_test">Fisher's exact test</a>.<br />
<br />
The chi-square test, in R:<br />
<pre class="literal-block">> abcd = matrix(c(4, 10, 1000, 10000), nrow=2, byrow=T)
> chisq.test(abcd)
Pearson's Chi-squared test with Yates' continuity correction
data: abcd
X-squared = 4.2691, df = 1, p-value = 0.03881</pre>
With the common p-value cutoff ("alpha") of 0.05, we'd say this page with 4 likes out of 10 views is significant — for that cutoff, at least. And if we applied the same test across all pages in the database, we'd be wrong.<br />
<br />
I'll try to be quick about this, because it matters.<br />
<br />
Remember: A p-value of 0.05 means the given like/view ratio will occur by chance 1 in 20 times. Since the same test is being applied to every page in your database, you need to account for <b>multiple hypothesis testing</b>, or else many pages will meet the cutoff by chance.<br />
<br />
If you only have a few pages — say, less than 40 — then you can divide <i>alpha</i> by the number of pages (e.g. 20) and use that in place of the original cutoff (0.05 / 40 = 0.00125, so the previous p-value of 0.03881 would <i>not</i> be significant).<br />
<br />
More likely, you have many more pages than that — hence the need to use grown-up statistics in the first place. Bonferroni correction (described above) would produce a cutoff that's much too stringent, so you'll need a more powerful method.<br />
<br />
R makes this easy. Starting with a single array of the p-values from Chi-square tests of each page:<br />
<pre class="literal-block">> pvals = sapply(chisq.test, contingencytables)</pre>
Adjust these raw p-values for multiple testing (using the <a class="reference external" href="http://en.wikipedia.org/wiki/Familywise_error_rate">familywise error rate</a>, by default — read the help page for p.adjust for all the details):<br />
<pre class="literal-block">> pvals.adj = p.adjust(pvals)</pre>
What do you after this depends on your own code. You can get a boolean array signifying which adjusted p-values are now smaller than alpha, which is useful for selecting "significant" pages from the original page list:<br />
<pre class="literal-block">> significant = pvals.adj < 0.05</pre>
Note that this selects both significantly liked and significantly disliked pages at the same time. To distinguish between the two, just compare each page's like/view ratio to the global average and select higher or lower.<br />
<br />
Another note about the contingency table: Once your application has counted a very large number of site-wide likes and views (cells C and D), this test will register significance for almost any page. You might have better results by replacing the global view and "Like" counts with a per-month or per-user average. And, you can cache these values and update them only occasionally.<br />
<h4>
Trending<br />
</h4>
Calculating p-values is a lot more work than selecting the top and bottom quantiles. If you've put in the extra effort, here's another feature you can support: a list of newly significant winners and losers.<br />
<br />
Each day (or hour or so), perform the chi-square test (described above) across all pages and note which ones cross the significance threshold. Compare this to the previous run's results to see which pages have crossed over, and add these newly significant hits to a separate chart — "Trending", I'll call it.<br />
<br />
This chart shows the pages that have just recently been determined to be likeable, but (probably) haven't accumulated enough votes to reach the "Top Pages" chart. It's a timelier list than "Top Pages", though the average quality of the "Trending" pages is not as high. This is the place where memes show up first. If they're truly good content, they'll eventually make it onto "Top Pages" — but that's not usually the case with memes.<br />
<br />
I'd treat the "Trending" chart as a queue, adding newly trending pages to the top at the end of each run and dropping pages from the bottom as space permits. Or just keep it rolling by week, like a blog. By adjusting <i>alpha</i> you can tune the number of newly significant pages found in each run, and therefore the turnover rate of your "Trending" queue.<br />
<br />
<h3>
5. But is it general?<br />
</h3>
Under the Google Whatever model (YouTube, Picasa, etc.), the ratio of Likes to total views for any given page is small. The statistics here will work in other cases, though — for example, an "Approve" button which clicked most of the time, or a "Dislike" button in place of the "Like" button. In the case where users have to click either "Like" or "Dislike" (Yes/No, Yay/Nay, or any other two options), this is also fine; just pick one option to count, and count "views" as the sum of likes and dislikes.<br />
<br />
<b>Update:</b> The "Like" event can also be something implicit, like downloads or completed interactions (out of started interactions). This opens up options for sites that don't require users to register.<br />
<br />
What about sites with 5-star ratings, like Amazon? Well, there's an easy way and a hard way. Easy: count the ratings as fractional Likes ([0, .25, .5, .75, 1] if you allow 0 stars, [0, .33, .67, 1] if you don't), and use the pseudocounts just like before. The hard way is to treat each star ranking as a separate event category — but that's going to have to wait for a later post.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com2tag:blogger.com,1999:blog-266234734515043410.post-35054287219424683872010-10-05T11:23:00.000-04:002012-10-27T14:43:12.194-04:00Bio.Phylo: A unified phylogenetics toolkit for BiopythonI presented this at the Bioinformatics Open Source Conference (<a href='http://www.open-bio.org/wiki/BOSC_2010'>BOSC 2010</a>) in early July, but somehow forgot to post it here too. It's an overview of my somewhat new sub-package for working with phylogenetic trees in Biopython, based on my Google Summer of Code 2009 project (a phyloXML parser in Biopython).<br />
<br />
In a nutshell, Bio.Phylo is a library for manipulating finished phylogenetic trees and integrating them into a Biopython-based workflow. It can handle the standard file formats — Newick, Nexus and phyloXML, with the current exception of NeXML — and has particularly good support for <a href="http://phyloxml.org">phyloXML</a>.<br />
<br />
This presentation walks through an example of loading a Newick tree, viewing it a few different ways, adding branch colors, and saving it as a phyloXML file.<br />
<br />
<div style="width:425px" id="__ss_4809399"><strong style="display:block;margin:12px 0 4px"><a href="http://www.slideshare.net/etalevich/biophylo-phylogenetics-in-biopython-bosc-2010" title="Bio.Phylo: Phylogenetics in Biopython (BOSC 2010)">Bio.Phylo: Phylogenetics in Biopython (BOSC 2010)</a></strong><object id="__sse4809399" width="425" height="355"><param name="movie" value="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=bosc2010-biophylo-talevich-100721205437-phpapp01&stripped_title=biophylo-phylogenetics-in-biopython-bosc-2010&userName=etalevich" /><param name="allowFullScreen" value="true"/><param name="allowScriptAccess" value="always"/><embed name="__sse4809399" src="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=bosc2010-biophylo-talevich-100721205437-phpapp01&stripped_title=biophylo-phylogenetics-in-biopython-bosc-2010&userName=etalevich" type="application/x-shockwave-flash" allowscriptaccess="always" allowfullscreen="true" width="425" height="355"></embed></object><div style="padding:5px 0 12px">View more <a href="http://www.slideshare.net/">presentations</a> from <a href="http://www.slideshare.net/etalevich">Eric Talevich</a>.</div></div><br />
The conference abstract is <a href='http://www.open-bio.org/w/images/6/6a/8_BOSC2010.pdf'>here</a>. I also recommend the main documentation in the <a href="http://biopython.org/DIST/docs/tutorial/Tutorial.html">Biopython Tutorial</a> (see chapter 12) and the <a href="http://www.biopython.org/wiki/Phylo">wiki page</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-78797447418977958572010-04-08T17:12:00.003-04:002011-03-01T15:15:49.105-05:00Google Summer of Code 2010: The final draftThe <a href="http://socghop.appspot.com/gsoc/program/home/google/gsoc2010">Google Summer of Code 2010</a> application period is in its final 24 hours.<br />
<br />
I volunteered to mentor with two organizations this year, <a href="http://www.open-bio.org/wiki/Main_Page">OBF</a> and <a href="http://www.nescent.org/index.php">NESCent</a>. Last month I posted a couple of ideas with each org:<br />
<ul><li><a href="https://www.nescent.org/wg_phyloinformatics/Phyloinformatics_Summer_of_Code_2010#Jump-start_MIAPA_protocol_annotation_with_a_user-accessible_demo">Jump-start MIAPA protocol annotation with a user-accessible demo</a></li>
<li><a href="https://www.nescent.org/wg_phyloinformatics/Phyloinformatics_Summer_of_Code_2010#Biopython_and_PyCogent_interoperabilityhttps://www.nescent.org/wg_phyloinformatics/Phyloinformatics_Summer_of_Code_2010#Biopython_and_PyCogent_interoperability">Biopython and PyCogent interoperability</a></li>
<li><a href="http://biopython.org/wiki/Google_Summer_of_Code#PDB-Tidy:_command-line_tools_for_manipulating_PDB_files">PDB-Tidy: command-line tools for manipulating PDB files </a></li>
<li><a href="http://biopython.org/wiki/Google_Summer_of_Code#Integration_with_a_third-party_structural_biology_application">Biopython integration with a third-party structural biology application</a></li>
</ul><br />
The applications that have come in have been pretty good; the only thing I can complain about is that nobody has followed through with my MIAPA project -- we got a nibble from one student, but nothing after that.<br />
<br />
Since we're doing the last round of application reviews now before the deadline, here's some general guidance on what mentors are looking for in a student application.<br />
<br />
First, a couple of outside references:<br />
<ul><li><a href="http://socghop.appspot.com/document/show/program/google/gsoc2009/studentallocations">Google's notes on student allocations</a></li>
<li><a href="http://bcbio.wordpress.com/2010/03/26/biopython-projects-for-google-summer-of-code-2010/">Blue Collar Bioinformatics: Biopython projects for GSoC 2010</a></li>
</ul><br />
<h4>The Zen of GSoC</h4>Google Summer of Code is a program to recruit and foster new long-term open-source contributors.<br />
<br />
Broadly, the mentoring organizations are asking three questions:<br />
<ol><li>Are you motivated enough about this work to continue contributing after the summer?<br />
</li>
<li>Can you write useful code on your own?<br />
</li>
<li>Do you interact well with the community, so that we can work with you to merge your work cleanly into the trunk and rely on you to maintain the codebase?</li>
</ol>You can get a sense of what Google and the mentoring orgs are looking for from the applications the orgs themselves submit to Google. For example: <a href="http://docs.google.com/View?id=dhdjhbvd_10c898hdhc">NESCent's 2010 app</a>.<br />
<br />
Here are some specific tips for demonstrating that you have some committer in you.<br />
<br />
<h4>Put your previous work online</h4>It's remarkable how many ostensible programmers just can't write decent code. They'll have a list of successful past projects they worked on, maybe a legitimate degree in computer science, but their code itself was clearly never fully understood by anyone, original programmer included. (Remember, programming languages exist for humans to understand -- the computer itself runs on machine code.) The only way we can be sure you can write code we can use is if we can look at something you've written previously.<br />
<br />
Biopython uses GitHub for development, so putting a project of your own on GitHub demonstrates two useful things: you can write functioning code, and you're already up to speed with the build tools that Biopython uses.<br />
<br />
If the most relevant code you've written is tied up in some way -- say, it's part of a research project still being prepared for publication -- see if you can use at least a few snippets of it. So far, it seems most professors have been willing to allow that.<br />
<br />
<h4>Subscribe to your mentoring organization's mailing list</h4>I know, e-mail mailing lists seem at least a decade behind the times. But open-source projects like to have a permanent public record of the discussions that happen, and everyone has an e-mail account. We also have IRC channels and Twitter tags (#phylosoc and #obfsoc), but project proposals are generally more than 140 characters so it's best to use e-mail at some point.<br />
<br />
Plus, you'll be able to read all the advice the other students are getting -- mentors get fatigued as the application season wears on, and once we've written the same thing a few times we start skipping details.<br />
<br />
<h4>Write a weekly project schedule</h4>The GSoC application has fields for pointing to external info. Create a Google document or spreadsheet (or README.md on GitHub if you're fancy) detailing your project plan week-by-week.<br />
<br />
Suggested fields:<br />
<ul><li>Date, or week number for referencing later</li>
<li>GSoC events and guidelines (see the <a href="http://socghop.appspot.com/document/show/gsoc_program/google/gsoc2010/faqs#timeline">official timeline</a>)</li>
<li>Deliverables for the week — what's produced, e.g. documentation sections, unit tests, classes, modules</li>
<li>Approach for each of these tasks, in a few words</li>
<li>Potential problems that could occur, specific to the tasks — perhaps a dependency turns out to be inadequate, or an integration step is required</li>
<li>Proposed mitigation for each of the foreseen issues</li>
</ul><br />
(If you want to estimate the number of hours or days each task will take, that's cool too.)<br />
<br />
Here are the examples from previous GSoC projects that we've been sharing on the mailing lists:<br />
<ul><li><a href="http://spreadsheets.google.com/pub?key=puFMq1smOMEo20j0h5Dg9fA&single=true&gid=0&output=html">phyloXML in Biopython</a> (mine)</li>
<li><a href="http://www.nescent.org/wg_phyloinformatics/PhyloSoC:PhyloXML_support_in_BioRuby">phyloXML in BioRuby</a></li>
<li><a href="https://www.nescent.org/wg_phyloinformatics/PhyloSoC:Phylogenetic_XML">Phylogenetic XML</a> (the origin of NeXML)</li>
<li><a href="https://www.nescent.org/wg_phyloinformatics/PhyloSoC:Build_a_Mesquite_Package_to_view_Phenex-generated_Nexml_files">NeXML for Mesquite</a></li>
</ul><br />
<h4>Respect the deadlines</h4>Submit a draft of your application to Google at least a day before the deadline, April 9. There are thousands of applicants each year, and Google has no reason to let the deadline slide — an important function of the application process itself is to screen out students who won't deliver by the stated deadlines. In effect, if your application isn't submitted to Google by noon PST on April 9, then you didn't apply.<br />
<br />
<i>BUT:</i> If you submit something even partially complete, we can contact you later during the review stage and get the remaining information from you. And if you included a link to your weekly plan (as a separate online document), you can edit that after the deadline too.<br />
<br />
Best of luck!etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-23270189340406539322010-02-24T12:29:00.003-05:002010-02-24T13:34:53.671-05:00Python workshop #2: BiopythonAs promised, here are the slides from Monday's Biopython programming workshop:<br /><br /><div style="width: 425px;" id="__ss_3266543"><strong style="margin: 12px 0pt 4px; display: block;"><a href="http://www.slideshare.net/etalevich/biopython-programming-workshop-at-uga" title="Biopython programming workshop at UGA">Biopython programming workshop at UGA</a></strong><object height="355" width="425"><param name="movie" value="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=biopywork-100224112524-phpapp01&stripped_title=biopython-programming-workshop-at-uga"><param name="allowFullScreen" value="true"><param name="allowScriptAccess" value="always"><embed src="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=biopywork-100224112524-phpapp01&stripped_title=biopython-programming-workshop-at-uga" type="application/x-shockwave-flash" allowscriptaccess="always" allowfullscreen="true" height="355" width="425"></embed></object><div style="padding: 5px 0pt 12px;">View more <a href="http://www.slideshare.net/">presentations</a> from <a href="http://www.slideshare.net/etalevich">Eric Talevich</a>.</div></div><br /><br />This was another 2-hour session, with a short snack break in the middle this time -- which was also a nice opportunity to ask everyone about the pacing, and see if who's been following along with the examples in IPython (versus staring at a BSOD or lolcats -- which I didn't notice any of).<br /><br />This went well:<ul><li>Pacing</li><li>Using IPython to inspect objects and display documentation -- this lets some people "read ahead" and perhaps answer their own minor questions, leading to other, better questions</li><li>The general introductory pattern of:<ol><li>Demonstrate how to import a module and instantiate the basic class</li><li>Review, in English, the core features of the module and why they exist</li><li>Walk through a short script that uses real data to accomplish some simple but useful task(s)</li><li>Display the result, completing the mental pipeline of <span style="font-style: italic;">input</span> -> <span style="font-style: italic;">transformation</span> -> <span style="font-style: italic;">output</span></li></ol></li></ul>Room for improvement:<ul><li>I didn't always execute the final draft of each example, so there were a couple typos -- inconvenient for those following along in Python. (I've fixed them in the slides here.)<br /></li><li>Consequently, I didn't have an output file to show at the end of each example -- so I had to describe or draft one on the spot.</li><li>The PDB module was the coolest part of the workshop, and I rushed it a bit. I was afraid the visitors from Genetics and Plant Bio would be bored with it, but I don't think they were, and the Bioinformatics folks were left wanting more.</li></ul>I'm planning to host both Python workshops again in the next academic year, either 1 per semester (as it was this year) or both each semester, maybe 2 weeks apart. The Biopython workshop in particular will be different next time because <a href="http://www.biopython.org/wiki/Phylo">Bio.Phylo</a> will finally be included with the main Biopython distribution -- evolution is cool, and more of the pretty is always a good thing to have in a programming workshop.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com0tag:blogger.com,1999:blog-266234734515043410.post-47801256804024967042010-02-19T15:52:00.004-05:002010-02-24T13:36:54.484-05:00Python workshop #1, now on SlideShareLast November I hosted a workshop on basic Python programming at UGA. The attendees were mostly from the bioinformatics department, but this workshop didn't go into science at all -- just practical Python usage. Today I finally got around to cleaning up the slides and uploading them to SlideShare:<br /><br /><div style="width: 425px; text-align: left;" id="__ss_3228053"><a style="margin: 12px 0pt 3px; font-family: Helvetica,Arial,Sans-serif; font-style: normal; font-variant: normal; font-weight: normal; font-size: 14px; line-height: normal; font-size-adjust: none; font-stretch: normal; display: block; text-decoration: underline;" href="http://www.slideshare.net/etalevich/python-workshop-1-uga-bioinformatics" title="Python workshop #1 at UGA">Python workshop #1 at UGA</a><object style="margin: 0px;" height="355" width="425"><param name="movie" value="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=pywork1-100219144829-phpapp02&rel=0&stripped_title=python-workshop-1-uga-bioinformatics"><param name="allowFullScreen" value="true"><param name="allowScriptAccess" value="always"><embed src="http://static.slidesharecdn.com/swf/ssplayer2.swf?doc=pywork1-100219144829-phpapp02&rel=0&stripped_title=python-workshop-1-uga-bioinformatics" type="application/x-shockwave-flash" allowscriptaccess="always" allowfullscreen="true" height="355" width="425"></embed></object><div style="font-size: 11px; font-family: tahoma,arial; height: 26px; padding-top: 2px;">View more <a style="text-decoration: underline;" href="http://www.slideshare.net/">presentations</a> from <a style="text-decoration: underline;" href="http://www.slideshare.net/etalevich">etalevich</a>.</div></div><br /><br />It looks like LaTeX Beamer and SlideShare's PDF/Flash converter don't play well together. Meh, it's still easy enough to read.<br /><br />I'm working on a Biopython-specific followup right now for a workshop on Monday, 2/22. I'll post that here when it's done, too, with reasonable haste.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com1tag:blogger.com,1999:blog-266234734515043410.post-18063153707289190082009-07-20T20:52:00.001-04:002012-05-23T13:42:49.128-04:00Faster string concatenation in Python<a href="https://www.nescent.org/wg/phyloinformatics/index.php?title=Phyloinformatics_Summer_of_Code_2009#Biogeographical_Phylogenetics_for_BioPython">Nick Matzke</a> pointed me to this discussion of string concatenation approaches in Python:<br />
<br />
<a href="http://www.skymind.com/%7Eocrow/python_string/">Efficient String Concatenation in Python</a><br />
<br />
The issue here is whether adding strings together in a <tt>for</tt> loop is inefficient enough to be worth working around. Python strings are immutable, so this:<br />
<pre>
s = 'om'
for i in xrange(1000):
s += 'nom'</pre>
<br />
means doing this 1000 times:<br />
<ol><br />
<li>Translate the assignment to "<span style="font-family: "Courier New",Courier,monospace;">s = s + 'nom'</span>"</li>
<br />
<li>Allocate another string, <span style="font-family: "Courier New",Courier,monospace;">'nom'</span>. (Or reuse a reference if it's already interred.)</li>
<br />
<li>Call the __add__ method on <span style="font-family: "Courier New",Courier,monospace;">s</span>, with <span style="font-family: "Courier New",Courier,monospace;">'nom'</span> as the argument</li>
<br />
<li>Allocate the new string created by __add__</li>
<br />
<li>Assign a reference to that string back to <span style="font-family: "Courier New",Courier,monospace;">s</span></li>
</ol>
<br />
So, using the + operator 1000 times in a loop has to create 1000 ever-larger string objects, but only the last one gets used outside the loop. There are good reasons Python works this way, but still, there's a trap here in an operation that gets used a lot in ordinary Python code.<br />
<br />
There are a few ways to cope:<br />
<ul><br />
<li>Use a mutable object to build the string as a sequence of bytes (or whatever) and then convert it back to a Python string in one shot at the end. Reasonable intermediate objects are array and StringIO (preferably cStringIO).</li>
<br />
<li>Let the string object's <tt>join</tt> method do the dirty work -- strings are a basic Python type that's been optimized already, so this method probably drops down to a lower level (C/bytecode in the CPython interpreter, not sure about the details) where full allocation of each intermediate string isn't necessary.</li>
<br />
<li>Build a single format string and interpolate with the <tt>%</tt> operator (or the format method, if you're fancy) to fill it in, under the same rationale as with the <tt>join</tt> method. This fits real-world scenarios better — filling in a template of a plain-text table or paragraph with computed values, either all at once with <tt>%</tt> or incrementally with string addition. It could be a performance bottleneck, and it's not obvious which approach would be better.</li>
</ul>
<br />
The original article gives a nice analysis and comes out in favor of intermediate cStringIO objects, with a list comprehension inside the string join method as a strong alternative. But it was written in 2004, and Python has changed since then. Also, it doesn't include interpolation among the tested methods, and that was the one I was the most curious about.<br />
<br />
<h2>
Methods</h2>
<br />
I downloaded and updated the script included with that article, and ran it with Python 2.6 and 2.5 to get some new results. (Source code <a href="http://www.relapsecollapse.com/static/strcat.txt">here</a>.)<br />
<br />
First, a changelog:<br />
<ul><br />
<li>The method numbers are different, and there are a couple more. Method #2 is for the <tt>%</tt> operator, in which I build a gigantic format string and a gigantic tuple out of the number list, then smash them together. It trades memory for CPU time, basically. Method #8 uses <tt>map</tt> instead of a list comprehension or generator expression; no lambda is required and the necessary function (<tt>str()</tt>) is already available, so this is a good candidate.</li>
<br />
<li>I used the standard lib's <tt>time.clock()</tt> to measure CPU time around just the relevant loop for each string concatenation method.</li>
<br />
<li>Fetching the process memory size is similar but uses the subprocess module and different options.</li>
<br />
<li>Docstrings are (ab)used to identify the output.</li>
</ul>
<br />
For example, the string addition method now looks like this:<br />
<pre>
def method1():
"""1. string addition"""
start = clock()
out_str = ''
for num in NUMS:
out_str += str(num)
cpu = clock() - start
return (out_str, cpu, memsize())</pre>
<br />
<br />
<h2>
Results</h2>
<br />
Each method concatenates the string representation of the numbers 0 through 999,999. The methods were run sequentially in separate processes, via a for loop in the shell, for Python versions 2.5 and 2.6. The best of three runs for each method are shown below.<br />
<b>Python 2.6:</b><br />
<pre>
1. string addition CPU (s): 1.99 Mem (K): 11.7
2. %-interpolation CPU (s): 2.42 Mem (K): 23.0
3. array object CPU (s): 3.42 Mem (K): 17.3
4. cStringIO object CPU (s): 3.24 Mem (K): 19.7
5. join + for loop CPU (s): 2.29 Mem (K): 48.0
6. join + list comp CPU (s): 1.93 Mem (K): 11.6
7. join + gen expr CPU (s): 2.08 Mem (K): 11.6
8. join + str map CPU (s): 1.47 Mem (K): 11.6</pre>
<br />
The winner is <tt>map</tt>, with string addition, the list comprehension, and the generator expression also doing well. String addition in a loop did much better than would be expected from reading the original article; the Python developers have put effort into making this less of a trap. Specifically, there's a flag on string objects internally that indicates whether the string is the result of an addition operation. This helps the interpreter identify when a string is being concatenated in a loop, and optimize that case by performing in-place concatenation. Nice. So really, there's no need to worry about the quadratic time behavior that we expected — at least in Python 2.6.<br />
<br />
The array object, a sequence of packed bytes, is supposed to be a low-level but high-performance workhorse. It was embedded in the minds of performance-conscious Python programmers by this essay by Guido van Rossum:<br />
<br />
<a href="http://www.python.org/doc/essays/list2str.html">Python Patterns — An Optimization Anecdote</a><br />
<br />
At a glance, that problem looks similar to this one. However, converting ints to chars is a problem that can be described well in bytes. Converting integers to their string representation is not — we're not even using any features of the array object related to byte representation. Going low-level doesn't help us here; as Guido indicates in his conclusion, if you keep it short and simple, Python will reward you. The StringIO object in method 5 performs similar duties, and the shape of both functions is the same; the only difference in performance seems to be that cStringIO trades some memory space for CPU time.<br />
<br />
The string join method is recommended by the Python standard library documentation for string concatenation with well-behaved performance characteristics. Conveniently, <tt>str.join()</tt> accepts any iterable object, including lists and generator expressions. Method 5 is the dump approach: build a list in a for loop, pass it to <tt>join</tt>. Method 6 pushes the looping operation deeper into the interpreter via list comprehension; it saves some bytecode, variable and function lookups, and a substantial number of memory allocations.<br />
<br />
Using a generator expression in method 7 instead of a list comprehension should have been equivalent or faster, by avoiding the up-front creation of a list object. But memory usage is the same, and the list comprehension runs faster by a small but consistent amount. Maybe <tt>join</tt> isn't able to take advantage of lazy evaluation, or is helped by knowing the size of the list object early on... I'm not sure. Interesting, though. In Python 3, the list comprehension is equivalent to building a list object from a generator expression, so results would probably be different there.<br />
<br />
Finally, in method 8, <tt>map</tt> allows the interpreter to look up the <tt>str</tt> constructor just once, rather than for each item in the given sequence. This is the only approach that gives an impressive speedup over string addition in a loop. So how portable is this result?<br />
<br />
<b>Python 2.5:</b><br />
<pre>
1. string addition CPU (s): 3.77 Mem (K): 10.8
2. %-interpolation CPU (s): 2.43 Mem (K): 22.0
3. array object CPU (s): 5.16 Mem (K): 16.4
4. cStringIO object CPU (s): 4.93 Mem (K): 18.7
5. join + for loop CPU (s): 3.98 Mem (K): 47.1
6. join + list comp CPU (s): 3.30 Mem (K): 10.5
7. join + gen expr CPU (s): 3.59 Mem (K): 10.5
8. join + str map CPU (s): 2.72 Mem (K): 10.5</pre>
<br />
Python 2.6.2 has had the benefit of additional development time, and in notably the <a href="http://code.google.com/p/unladen-swallow/wiki/ProjectPlan#2009_Q1">Unladen Swallow</a> project's first quarter of interpreter optimizations, with impressive improvements across the board. By comparison, Python 2.5 uses generally less memory and more CPU time. String interpolation, however, seems to already have been optimized to the max in Python 2.5, and actually wins the performance shootout here! String addition, on the other hand, is slightly less adept at optimizing in a loop. It still avoids the quadratic-time issue (that enhancement was added in Python 2.4), and memory usage is quite respectable.<br />
<br />
<h2>
Conclusion</h2>
<br />
The recommendations at the end of Guido's essay are still exactly right. In general, Python performs best with code that "looks right", with abstractions that fit the problem and a minimum of branching and explicit looping.<br />
<ul><br />
<li>Adding strings together in simple expressions will be optimized properly in recent Pythons, but could bite you in older ones</li>
<br />
<li>Using string interpolation or templates plays well enough with more complex formatting</li>
<br />
<li>Going too low-level can deprive you of Python's own optimizations</li>
<br />
<li>If built-in functions can do what you need, use them, and basic Haskell-style functional expressions can make your code very concise</li>
</ul>
<br />
There's more discussion on <a href="http://stackoverflow.com/questions/376461/string-concatenation-vs-string-substitution-in-python">Stack Overflow</a>.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com4tag:blogger.com,1999:blog-266234734515043410.post-4696150280173323392009-03-08T18:50:00.000-04:002012-10-27T14:46:02.251-04:00Mnemosyne: Getting Things MemorizedIt had been bothering me since I joined this lab that I couldn't confidently just read a protein sequence and understand what it meant — naming the residues, picturing the side chain structures, and understanding the significance of replacing one residue with another. I expected that I'd just pick it up naturally from working with sequences and structures, and that did happen somewhat. But I wanted it to be as easy as reading English, and that level of completeness doesn't happen without some rote memorization.<br /><br />That brought to mind a <a href="http://www.wired.com/medtech/health/magazine/16-05/ff_wozniak">Wired article</a> about Piotr Wozniak and his spaced-repetition memorization program, <a href="http://www.supermemo.com/">SuperMemo</a>. When I originally read the article I wasn't in grad school and didn't have an urge to memorize any particular list of things. Anyway, SuperMemo appeared to be Windows-only software, and an algorithm like this would be more fun to code from scratch anyway. Enough fun, really, that there <em>had</em> to be one or two open-source implementations floating around.<br /><br /><a href="http://mnemosyne-proj.org/">Mnemosyne</a> popped up as the closest match in an Ubuntu package search, so I'm running with that. Putting the flash cards together was pretty simple; I was able to do it in a few minutes from inside the program and export it in the standard XML format. I zipped it up with a quick plain-text README and uploaded it to the project home page as the <a href="http://mnemosyne-proj.org/node/166">Amino Acids</a> card set.<br /><br />The content came from a slide in a lecture, and I did a quick sanity check on Wikipedia before uploading. The notation for the 20 standard amino acids is complete, and that was the main goal of this. The assignment of amino acid "groups" seems to be a little arbitrary, depending on the source (by structure, functional groups, chemical properties, etc.), and I tried to make the categories complete without too much overlap -- there's a small deviation from my slide here. I also added another category for "side chain properties", pH and polarity. Another enhancement might be the standard codons for each amino acid, though I'm not sure I want to deal with that yet.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com4tag:blogger.com,1999:blog-266234734515043410.post-5467270746305790252009-02-07T14:36:00.000-05:002020-02-22T15:57:35.362-05:00Carrots and sticksWhat's old is new again:<br />
<blockquote class="tr_bq">
<a href="https://www.theglobeandmail.com/news/national/professor-makes-his-mark-but-it-costs-him-his-job/article20440434/">Professor makes his mark, but it costs him his job</a></blockquote>
<div class="graf graf--p" name="c4d8">
In <em class="markup--em markup--p-em">Zen and the Art of Motorcycle Maintenance</em>, Robert Pirsig mentions his own experiment in withholding grades at a university. He didn’t just announce on the first day that everyone would get an A+ — that seems gimmicky. Instead, since it was a class on rhetoric, he spent the course developing an argument for eliminating the grades-and-degrees system and discussing it with his students. </div>
<div class="graf graf--p" name="b35d">
Initially, most students were unenthusiastic or opposed — grades and degrees are what they came for. Nonetheless, Prof. Pirsig assigned, collected and graded papers, but returned them to students with only the comments, not the grade.</div>
<div class="graf graf--p" name="b76c">
At first:</div>
<ul class="postList">
<li class="graf graf--li" name="cd2b">A-students felt annoyed by the uncertainty of the situation, but did the work anyway;</li>
<li class="graf graf--li" name="e35b">B-C students blew off some assignments; and</li>
<li class="graf graf--li" name="680d">C-D students usually skipped class.</li>
</ul>
<div class="graf graf--p" name="ddc6">
He observed this and changed nothing. If students acted up, he let it slide.</div>
<div class="graf graf--p" name="7f66">
Around 3–4 weeks into it:</div>
<ul class="postList">
<li class="graf graf--li" name="42e3">A-students got nervous and pushed themselves harder, in class and in papers;</li>
<li class="graf graf--li" name="bd30">B-C students saw what the A-students were doing and returned to the usual level of effort; and</li>
<li class="graf graf--li" name="8e54">C-D students who had made a routine of skipping class would occasionally show up out of curiosity.</li>
</ul>
<div class="graf graf--p" name="d59a">
And finally:</div>
<ul class="postList">
<li class="graf graf--li" name="4f57">A-students relaxed and began enjoying the class as active participants. In a final essay, still not knowing what their grades were, these students favored eliminating grades by 2–1.</li>
<li class="graf graf--li" name="a8ab">B-C students saw this, panicked, and began putting an unusual amount of effort into their work. Eventually, they joined the A students in engaging class discussions. Ultimately, these were evenly divided over the issue of eliminating grades.</li>
<li class="graf graf--li" name="8761">C-D students — or those who attended — also saw this and began trying to hand in reasonable work. Those who couldn’t hack it freaked even more, and remained in a state of Kafkaesque terror until the quarter mercifully ended. Naturally, in the final essay these students were unanimously opposed to eliminating grades.</li>
</ul>
<div class="graf graf--p" name="a13c">
Interesting as this result was, Pirsig reverted to the regular grading system the next quarter because he couldn’t provide any alternate goal for students — those who can recognize quality in their own work don’t need the university; those who can’t need something to work toward, or they don’t progress.</div>
etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com1tag:blogger.com,1999:blog-266234734515043410.post-1230509700381119662008-10-25T12:29:00.000-04:002012-10-27T14:46:02.249-04:00Psychology research with Mechanical Turk<p><b>Elevator pitch:</b> There's a missing gear in the machine of psychology research. Every significant human study requires weeks or months of data collection, and more time coding that data in a form that can be analyzed statistically. This makes it infeasible to do the sort of fast, iterative refinement of models that biology has seen in recent years.<br /><br />Amazon's <a href="https://www.mturk.com/mturk/welcome">Mechanical Turk</a> provides the missing piece. It provides an accessible interface for building a survey, interactive test, or other psychological measure, pushing it out to thousands of participants, quickly returning the results to the researcher in electronic form, and screening out unusable data. It's flexible enough to allow screening and debriefing, and gives access to a vastly larger pool of participants than Experimetrix. And it's cheap.</p><h2>Background</h2><p>First, take a look at this: <a href="http://www.tenthousandcents.com/index.html">Ten Thousand Cents</a>.<br /><br />When I first heard about bioinformatics I was under the impression that it was the exponentially increasing power of computers that made it irresistible to start using them for biological research. But actually, it was pretty much the reverse -- high-throughput experimental methods like gene sequencing, mass spectrometry and X-ray crystallography generated too much data for humans to process manually. Computers were only barely able to handle this workload in 1986, when the human genome project started -- scientists just did what was needed to move around the mountains of data coming out of their experiments. Similarly, new computational research is coming out of the Large Hadron Collider project now.<br /><br />Psychology researchers (especially in social psychology) currently spend semesters at a time gathering data for their studies and converting it into data that can be quantitatively analyzed. High-throughput experimental methods are scarce and expensive, so there's no "data glut" driving the development of better information-management methods. Progress in the field is slow and lossy -- since there's not much demand for the raw data, conclusions are described qualitatively, which makes it hard to use prior results as a solid foundation for future work.<br /><br />With Mechanical Turk, it's possible to do in one shot a study that would otherwise require a meta-analysis of several studies across particular locations or demographics. With more consistent data and larger populations, data <i>can</i> be reusable.</p><h2>How it could work</h2><p>If it fits behind a web interface, or can be described and completed with plain language or pictures, it can be done with Mechanical Turk. Necessarily, a form of consent can precede the main task, and a blurb of debriefing can finish.<br /><br />To get a feel for how it's done, read this article: <a href="http://waxy.org/2008/11/the_faces_of_mechanical_turk/">The Faces of Mechanical Turk</a><br /><br />Naturally, the first study done this way should be something to determine how the population of Turkers corresponds to the general population and the student populations that have already been characterized in previous studies. A public-domain measure of the Big Five or something like the Narcissistic Personality Inventory would be good candidates. Then, let slip the hounds of statistics. Are Turkers as representative of the general population as psychology undergrads? More so?<br /><br />Some research along these lines has already been blogged here: <a href="http://behind-the-enemy-lines.blogspot.com/2008/03/mechanical-turk-demographics.html">Mechanical Turk Demographics</a><br /><br />Now, let's try some examples.<br /><br /><b>Surveys:</b> You craft a survey, Turkers take it, and you retrieve and filter the results through the Mechanical Turk interface. Pretty straightforward, no?<br /><br /><b>Interactive tasks:</b> This is what Mechanical Turk is designed for; only, the focus was expected to be the task, not the Turker. Anyway, the data's yours. An example of a task like this would be a simple, unbounded game (Flash or JavaScript) that the participant can quit any time (possibly paired with another stimulus). The returned data would be the play duration alongside any personal or demographic information requested.<br /><br /><b>Coding visual or audio data:</b> Following the original intent of Mechanical Turk more closely, this application of the service distributes a repetitive task normally performed over several weeks or months by the researcher or a group of grad students. Rather collect new data about a participant, this simply boils down a vast quantity of data that's already been generated -- this is a problem we want to have. A two-step example: (1) run a Mechanical Turk task in which participants draw or assemble an arbitrary image; (2) run a second task with a different set of participants who look at these images and code (type or select) the relevant traits they see in the images.<br /><br /><b>Measure development:</b> One of the more uncomfortable questions in social psychology research is the validity of personality measures. Devise a series of questions and a method for tabulating the results; run it on some participants; analyze the results to get some answers. But, what's really being tested here -- the population, or the measure? Tragically, there's no time to refine the measure very much; if the results are useful, you run with it. But! With Mechanical Turk, collecting survey results is cheap and quick; and since the general format of the survey isn't changing between revisions, the same set of statistical transformations can be applied programatically to each iteration of the survey.<br /><br />This is a great way to build a psychological measure that you can be confident in: Push an initial draft of the measure out to Turkers, receive some results, perform a statistical analysis and save the operations as an <a href="http://www.r-project.org/">R</a> or SPSS script. Then, manually refine the measure, put it back on Turk, filter the new results through your analysis script, and repeat until it looks good. This can get as advanced as you'd like -- start with several times as many questions as you'd like to see in the final survey, then automatically dispatch random subsets of the question list to Turk, filter through your automatic analysis to get some scores indicating quality, and use a Bayesian classifier to narrow down the best possible subset of questions.</p><br /><p><b>Update:</b> Here's a conference paper on the same topic.<br /><a href="http://www-users.cs.umn.edu/~echi/papers/2008-CHI2008/2008-02-mech-turk-online-experiments-chi1049-kittur.pdf">http://www-users.cs.umn.edu/~echi/papers/2008-CHI2008/2008-02-mech-turk-online-experiments-chi1049-kittur.pdf</a></p>etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com1tag:blogger.com,1999:blog-266234734515043410.post-72755819633880961442008-03-07T15:48:00.002-05:002010-01-01T13:58:38.723-05:00Vimming your way to the topHere's the Vim syntax file I use for highlighting my to-do list. It's based on the syntax file for YAML.<br /><a href="http://www.vim.org/scripts/script.php?script_id=2599">http://www.vim.org/scripts/script.php?script_id=2599</a><br /><br />Benefits:<br /><ul><li>Different colors for lines ending in ':', or starting with '*' or '{'</li><li>Assign keywords to be automatically highlighted, like important locations, coworkers' names, customers, taquerias, etc.</li><li>Start sections with a line of underscores and a heading beginning with the '{' character. The heading stands out (red with GVim's "desert" color scheme), and you can jump between sections just like C blocks using ]] and [[ keystrokes.</li><li>Ordinary text (i.e. not specifically formatted for this syntax) looks sane.</li></ul>Normally I have a line in my .vimrc assigning the filetype "todolist" to the file where I keep my permanent todolist, but another way to add this highlighting to a text file is to add <tt>vim: ft=todolist</tt> to the end of a file. It's harmless.<br /><br /><b>Update (4/2/09):</b> I uploaded the script to vim.org, where it will be easier to track and update.<br /><br /><span style="font-weight: bold;">Update (1/1/10):</span> Here's an example of how to use this color scheme for course notes.<br /><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgq6M-voHuj5i845t61NIDjqu-mAio62nC8PUnOE5wvSws1UUMsyE2MrEwcaJlR9OVIw-gyGdJQg4U_1nwdjVjoSD8MXiMqnh4kE4UeEAD6qfta_Rm3vccI1ozkZ8pRHLZYdLbDjr0-UtGW/s1600-h/todolist-desert.png"><img style="margin: 0px auto 10px; display: block; text-align: center; cursor: pointer; width: 400px; height: 383px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgq6M-voHuj5i845t61NIDjqu-mAio62nC8PUnOE5wvSws1UUMsyE2MrEwcaJlR9OVIw-gyGdJQg4U_1nwdjVjoSD8MXiMqnh4kE4UeEAD6qfta_Rm3vccI1ozkZ8pRHLZYdLbDjr0-UtGW/s400/todolist-desert.png" alt="" id="BLOGGER_PHOTO_ID_5421841357448119970" border="0" /></a><br /><ul><li>60 underscores (my preference) and a curly brace indicate a new section</li><li>Subsection lines end with a colon (generally followed by bullet points)<br /></li><li>Special or out-of-context notes start with an asterisk</li><li>For separation, or to display a different sort of sub-heading, play with asterisks: '* * *' centered, or '** OLD **' for example</li></ul>At school, I run a shell script for each new class that creates a new directory from the course name, copies a skeleton of this example text to a file called lecture-notes.txt, etc., and adds the directory to Mercurial -- so while there's some boilerplate involved with this plugin, it's easy to automate and plays well with Vim's text-munging capabilities.<br /><br />I've also picked up the habit of putting @contexts above unsorted items at the top of my main to-do list, inspired by the GTD approach. The syntax plugin doesn't take advantage of this yet; I'll post another update when that's done.etalhttp://www.blogger.com/profile/10168388850793209768noreply@blogger.com2