Showing posts with label FASTA. Show all posts
Showing posts with label FASTA. Show all posts

Friday, May 09, 2014

A JavaScript Tip for Biohackers

Today I want to give a short bio-hacking code lesson, and I want to keep it as simple as possible in case you're just starting to use JavaScript and want to see what kinds of things it can do. This lesson could prove handy to you even if you're not a gene hacker, because it shows how easy it is to manipulate text in the browser.

For this example, I am going to assume that you have a text file open in the browser. In particular, I'm going to assume the text file consists of a listing of two genes sequences in FASTA format, which is actually very simple: FASTA consists of a header line that starts with the > (greater than) angle bracket, followed by one or more lines of AGTC base sequence data (if it's a nucleotide sequence) or a bunch of letters (KPRMIV etc.) if it's a protein sequence. For example:

>M. leprae MBLr_2224...    <-- this is the header
ATGGCGGTGCTGGATGTC...      <-- this is the data

Your file might have several genes in FASTA format, one after the other. The question is: How can you read this data with JavaScript?

Actually, it's extremely easy to read text data. Open any text file with your browser, then open a JavaScript window: With Firefox, use Shift-F4 to bring up the Scratchpad, or with Chrome use Control-Shift-J. To read the text into a JavaScript variable called text, just type the following line into the script editor:

text = document.body.textContent;

Execute this code with Control-L in Firefox or by simply entering a carriage return in Chrome. (If Firefox fills your Scratchpad with text, bracketed by /* and */,  that's good; it means the code worked. Delete it and proceed.)

Suppose you have several FASTA records in a row and you want to have an array of gene data. The easy thing to do is "split" the records at each header, discarding the header:

r = />[^\n]+\n/g;
genes = text.split( r );

NOTE: To enter more than one line of code in the Chrome console, you have to hold the Shift key down before hitting Enter. Otherwise, Enter executes the code.

The first line defines a regular expression (r) for the pattern: "greater-than symbol followed by one or more non-newline symbols, followed by a newline." (The caret symbol ^ means to negate whatever's in between the square brackets.) When this code executes, genes will be an array of data, but because of the way split() works, the first item (item zero) in the array will be empty, so get rid of it with:

genes.shift();

Now the genes array will contain gene data. The data for gene No. 1 will be in gene[0], the data for gene No. 2 will be in gene[1], etc.

Incidentially, if you want an array of headers, just do:

headers = text.match( r );

No need to do headers.shift(). The match() operation creates an exact array.

If your genes are aligned, you can compare them, base to base, in a loop. (If your genes are not aligned, create an alignment using an online ClustalW alignment tool or using the popular Mega6 program.) The following loop construct compares the first 300 bases in two genes, and tallies the differences according to whether the difference occurred in codon base one, base two, or base three:

snp=[0,0,0]; // array to hold base 1,2,3 results
gene1 = genes[0];
gene2 = genes[1];
 
for (var i=0; i < 300; i++)
   snp[ i % 3 ] += gene1[i] != gene2[i];

You can display the results in the console simply by adding (on its own line) snp; or console.log(snp). Or if you want to see it in a dialog, execute alert(snp).

The final line of code deserves explanation. Results are placed in the snp (single nucleotide polymorphism) array according to whether the "hit" occurred in base 1, base 2, or base 3 of a codon. The i % 3 construct (i modulo 3) means divide i by 3 and throw away the "answer" but keep the remainder. (So for example, 5 % 3 equals 2, 6 % 3 equals zero, 7 % 3 equals 1, etc.) As i increments, i % 3 simply takes on values of 0, 1, 2, 0, 1, 2, etc.

On the right side of the equals sign we have gene1[i] != gene2[i]. This means we want to compare the two genes at an offset of i, and if they are not equal (!=), tally the resulting value (true, or 1) numerically. The numeric result is added to the appropriate slot in snp. Note that JavaScript treats true as having a numeric value of one and false as having a numeric value of zero. It knows to cast the boolean value to a number because of the += symbol (which means numerically add the following value to the existing value of the lefthand variable).

I hope this short tutorial wasn't too painful. If you're new to JavaScript, my advice is: Experiment in the console (Firefox's Scratchpad or Chrome's JS console) a lot, and get familiar with the string functions split() and match(), because they can be incredibly useful, not to mention speedy. Either of those two functions can take a string argument or a regex. But remember to put a 'g' after the regex if you want the operation to occur globally throughout the string. For example:

"abcdefabcdef".match( /abc/ )

will only match the first occurrence of abc, whereas

"abcdefabcdef".match( /abc/g )


will match both occurrences of abc and give you an array of matches.

Sunday, May 05, 2013

More Science on the Desktop

Not to keep harping on the amazing power of desktop omics tools, but I thought I'd share a tip for those of you into genome-mining. The tip in a nutshell is that if you gang-load a bunch of FASTA sequences (DNA sequence data) into the FeatView form at http://genomevolution.org, then click the rather inconspicuous button labeled "Phylogeny.fr" at the bottom left of the FeatView page, you'll be taken automatically to http://www.phylogeny.fr, where you'll get a realtime-generated phylogenetic tree based on the sequence data you provided in FeatView, with no effort on your part (it's truly a one-click operation). Copy and paste DNA sequences into FeatView, click one button, and 30 seconds later a tree shows up on your screen, looking (perhaps) something like this:


The reason I made this tree is that I wasn't satisfied with my knowledge of the relatedness of certain weird microorganisms I've recently run into. Namely:
  • Ralstonia (which I mentioned yesterday), WEIRD BECAUSE: It turns hydrogen gas and CO2 into plastic.
  • Bordetella, a bronchial infection agent; WEIRD BECAUSE: It turns out to be very similar, genetically, to Ralstonia
  • Burkholderia, a soil organism (and human and animal pathogen), WEIRD BECAUSE: It has an unexpectedly large amount of genetic similarity to Ralstonia and Polynucleobacter
  • Polynucleobacter, a ditch-water bacterium, WEIRD BECAUSE: It can live as an intracellular parasite of freshwater ciliates or it can live independently in soil (making it potentially a great study organism for determining the genetic bases of intracellular symbiosis)
  • Thiomicrospira, a very tiny CO2- and sulfur-loving organism, WEIRD BECAUSE: It can only be found near deep-sea thermal vents (see my previous writeup)
  • Polaromonas, a relatively newly discovered and still poorly understood bacterium, WEIRD BECAUSE: It is abundant in glacier ice on multiple continents. Plus it has an amazing (and totally unexpected) amount of genetic overlap with our good friend Bordetella, the whooping-cough bug.
If you're not familiar with how bacterial classification works, let's just say it's a mess. There's a long historical tradition of classifying microorganisms based on a hodgepodge of ad hoc methods involving everything from physical appearance under the microscope (especially after staining with crystal violet), to the habitat of the organism, to its ability to metabolize various substances, its ability to make spores, adaptation to oxygen or lack of oxygen, serological characteristics, etc. It's always been an error-prone system, resulting in many misclassifications and later corrections, owing to its inconsistency and basic irrationality, to put it bluntly. With the advent of molecular genetic techniques, it's now possible to create accurate phylogenies based on little more than DNA sequence differences, usually involving the 16S ribosomal RNA (more here).

Freshwater ciliates (like this Euplotes) are
home for Polynucleobacter endosymbionts.
As big an advance as ribosome-based phylogeny is, it's pretty far from ideal (IMHO), mainly because it ignores phenotypes. In fact it's pretty far removed from anything at all having to do with an organism's ecology, metabolism, mode of living, etc. What are we really measuring when we measure relatedness according to a 16S ribosomal yardstick? Just the rate of random mutation accumulation in a pretty uninteresting cell artifact. I'd rather have a yardstick that's tied to phenotypic reality than to a slow-to-change, "highly conserved" piece of cold dead scaffolding.

So to create my own "family tree" of two dozen or so microbes, I said to hell with 16S ribosomes and decided to use, as my yardstick, genetic variation in the
GroEL gene, which codes for the 60-kiloDalton heat-shock protein. I chose this protein (or rather, the gene for it) as my phylo-yardstick for a number of reasons. First, the DNA sequence is sizable, at about 1643 nucleotides (making it somewhat bigger than the 16S rDNA). It's important to have a large yardstick gene when looking for faint genetic signals. Secondly, this protein is essentially universal in prokaryotes. It's ubiquitous but not necessarily highly conserved, in the same sense that rRNA is highly conserved. ("Highly conserved" is not what you want. Think about it. Taken to the extreme, a "highly conserved" sequence is invariant. It never changes. And is therefore useless for phylogenetics.) Thirdly, the GroEL heat-shock protein has multiple intracellular touchpoints: It's known to interact with GroES, ALDH2, and dihydrofolate reductase, and it's involved in signal tranduction (it's induced not just by heat but by hydrogen peroxide). Not to overlook the obvious, but it is also a touchpoint protein for any enzyme that can be repaired by the 60kDa heat shock protein. That's probably dozens if not hundreds of enzymes. Why is that important? Think about it: A protein that is sensitive to the 3D conformational requirements of other proteins has to evolve in response to the needs of all the proteins it services. A thermophile (Thermomicrospira)  is going to need a different heat-shock repair system than a psychrophile (Polaromonas). A salt-lover needs a different one than a freshwater-lover. GroEL has to reflect, in its own structure, the many shifting requirements of the host proteome. These considerations make GroEL a highly appropriate basis gene for phylogenetic analysis.

And frankly, I think the GroEL-based phylo-tree phylogeny.fr spit out for me (see illustration further above) speaks for itself. It's a remarkably informative (and accurate) tree. GroEL evolutionary differences not only accurately grouped endosymbionts together, soil organisms together, aquatic organisms, etc., it also correctly grouped the "enteric-alike" Erwinia with E. coli and Shigella, and it cannily put Polaromonas with soil organisms (rather than aquatics), which I think is correct, based on recent Polaromonas isolates being found in soil rather than snow. Likewise, it's good to see Bdellovibrio (a freshwater bug) clustered with Polynucleobacter (which is symbiotic with a ciliate protozoan), with Thiomicrospira (the saltwater hydro-vent organism) a very nearby out-node.

If you get an infection while in a hospital, pray
it's not Clostridium difficile, which is often deadly.
A harder call to make is Clostridium difficile, which is present in 1% to 5% of non-ill people's intestines. Is it an enteric (a la E. coli)? Definitely not. The Clostridia (botulism, tetanus, etc.) are spore-forming soil bacteria. Their placement in the tree not far from the soil-dwelling spore-former, Bacillus thuringensis, is thus eminently correct. Bacillus is a proximal out-node relative to Clostridium, which is understandable in that Bacillus is aerobic whereas Clostridia are strict anaerobes.

Buchnera
(an aphid symbiont) comes at an odd location, much further away from the insect-dwelling Wolbachia than I would have predicted, but then again Buchnera's host feeds on cold sap where Wolbachia's hosts typically feed on warm blood. All the organisms around Wolbachia in the tree are hemophiles.

Our good friend
Bordetella (of pertussis fame) is placed firmly in the soil group. I think that's real and significant. When you start to look at Bordetella's high DNA sequence similarity with Ralstonia and Burkholderia, it would be surprising, actually, if it fell anywhere else in the tree.

Honestly, when I took Bacterial Ecology 201 in college, many years ago, it was under duress and I hated the experience. But now, decades later, I'm starting to like it. With tools like those available for free at
http://genomevolution.org and http://www.phylogeny.fr, what's not to like?