Mailund on the Internet

I’ve spent the last two hours debugging scripts only to find that the error wasn’t in the scripts but in my analysis of the result…

I’m scanning a genome alignment for informative indels (indels where exactly two of four share a gap that starts and stops at the same position).  My scripts find the position of each of those and outputs it together with the two species having the gap: HC for human and chimp sharing the gap, HO for human and orangutan sharing the gap, HM for human and macaque sharing the gap, etc.

Now, in most of my analysis of this I do not want to distingush between which pair has the gap and which does not, I am only interested in the quartet topology (HC|OM vs. HO|CM for example). So I want to re-map the pairs CM to HO, CO to HM and OM to HC.

I do my analysis in R, and this is how I did the re-mapping:

data$pair <- factor(sapply(data$pair,
                           switch, HC="HC", HO="HO", HM="HM", CM="HO",CO="HM",OM="HC"))

The result is not what you’d expect:

> table(data$pair)

   CM    CO    HC    HM    HO    OM
  705   336 30377   646   349 13089
> data$pair <- factor(sapply(data$pair,
+                            switch, HC="HC", HO="HO", HM="HM", CM="HO",CO="HM",OM="HC"))
> table(data$pair)

   HC    HM    HO
13794 30726   982

There is clearly something wrong in the mapping.  The total number is correct, but HC now is not the sum of the earlier HC and CM!

This is how I should have done it:

data$pair <- factor(sapply(as.character(data$pair),
                           switch, HC="HC", HO="HO", HM="HM", CM="HO",CO="HM",OM="HC"))

Do you notice the difference?

The type of data$pair is factor so it is encoded as levels (numbers 1 to 6).  The switch function uses these leves as if they were integers, and use them to index into the list HC=”HC”, …, OM=”HC”.

If data$pair contained strings then switch would match them against the names in that list, but when it is integers it doesn’t.

The type really matters here.

> data$pair <- factor(sapply(as.character(data$pair),
+                            switch, HC="HC", HO="HO", HM="HM", CM="HO",CO="HM",OM="HC"))
> table(data$pair)

   HC    HM    HO
43466   982  1054

I just read this post at Gene Expression.

The opening paragraph is:

You might have heard that Richard Dawkins’ website has been blocked in Turkey because of that moron Harun Yahya (H/T Thabet). Here’s the justification:

His press assistant, Seda Aral, said: “We are not against freedom of speech or expression but you cannot insult people. We found the comments hurtful. It was not a scientific discussion. There was a line and the limit has been passed. We have used all the legal means to stop this site. We asked them to remove the comments but they did not.”

(emphasis from the link).

I don’t know what to say to that, really. You have freedom of expression, but you cannot say anything that is insulting?

Bayblab posted on this news story:

The body of an aristocrat who died nearly 90 years ago has been exhumed in the hope that it will help scientists combat a future flu pandemic.

Didn’t they already do roughly the same thing back in 2005?

Characterization of the Reconstructed 1918 Spanish Influenza Pandemic Virus
Tumpey et al. Science 310(5745) pp. 77 – 80, 2005

The pandemic influenza virus of 1918–1919 killed an estimated 20 to 50 million people worldwide. With the recent availability of the complete 1918 influenza virus coding sequence, we used reverse genetics to generate an influenza virus bearing all eight gene segments of the pandemic virus to study the properties associated with its extraordinary virulence. In stark contrast to contemporary human influenza H1N1 viruses, the 1918 pandemic virus had the ability to replicate in the absence of trypsin, caused death in mice and embryonated chicken eggs, and displayed a high-growth phenotype in human bronchial epithelial cells. Moreover, the coordinated expression of the 1918 virus genes most certainly confers the unique high-virulence phenotype observed with this pandemic virus.

The Spanish flu

The Spanish flu is a nasty virus. It’s an influenza that from 1918 to 1920 killed between 20 and 100 million people worldwide.  In comparison, only about 20 million were killed in the entire First World War, so if you look at deaths per year in the early 20th century, it is the Spanish flu that will show a spike, more than the war.

It was a truly global epidemic.  It is called the Spanish flu because the Spanish press wrote about it (most other countries had the press censored during the war, but Spain wasn’t involved in the war), but really it was found everywhere across the globe, with as much as 20% of the world population infected and killing 2.5-5% of the worldwide population.

With the danger of sounding like a mad scientist, I am really fascinated by this epidemic.

Digging up victims, to try to recreate the flu, might not be the brightest idea ever, but the reason people want to do this is to learn more about influenza epidemics.  With all the scare about the bird flu, we might need all the knowledge about influenza we can get.

There is absolutely nothing preventing a new worldwide epidemic like the Spanish flu, and we really want to be prepared if that happens.

Reconstructing the Spanish flu

I’m not sure exactly what the plan is for Sir Mark Sykes — the guy they are digging up now — but what they did in the Science paper was to re-create the Spanish flu virus from the samples they got there.

They managed to get the sequence of the Spanish flu virus genome and then combine it with (the modern) H1N1 virus, to get viruses with more or less of the Spanish flu in them.  These they then grew to learn about the virus.

Knowing the sequence of the Spanish flu will let us analyse the evolution of it and maybe tell us where it came from (jumping from another species; recombining or mutating to that particular variant; etc.)

Actually growing the virus lets us experiment with it, learn how deadly it really is and how it works (see e.g. this story), maybe develop vaccines (not that that would be much help against a new epidemic variant, but we might learn something from it).

There’s a very interesting post on Google’s blog: The Future of Mobile.

Where is the mobile phone heading? How will it change the future?

It is a quiet Saturday morning.  I am slightly hung over.  My scripts are scanning through a genome and I am just sitting here waiting for them to finish with nothing much to do.

So I started thinking.  How much does it actually cost to get a new genome, these days? If I wanted my own genome sequenced, how much would I have to pay and how long would it take to get it?

The (first) human genome project cost about $3 billion (about $300 million for Celera) and took about 10 years (1990 to 2000 for the first assembly, then three more years for completion, but let’s just say 10 here).

Now they want to sequence 1000 humans in three years for $30-50 million. Next generation sequencing techniques lowered the cost of that project by a factor of 10. Of course, it helps a lot to have the original genome to assemble up against as well.

I’ve asked Roald about the price for the first “arab genome”, but I haven’t gotten an answer yet.  I guess he doesn’t work Saturday mornings ;-)

The genomics age

Some people say we are in the “post genomic” age, but really we are just in the middle of the genomics age if anything.  We are seeing an explosion in new genomes sequenced.

From GOLD you can download some statistics on genome projects. Plotting the total number of genomes published against years, you clearly see the explosive increase in data:

It is even more impressive when you consider all genome projects and not just the published genomes so far:

Statistics at NCBI says we have 22 complete Eukaryote genomes, 161 with a draft assembly and 176 in progress. For Prokaryotes the numbers are 749 complete, 540 draft and 676 in progress.

It doesn’t say anything about the cost of sequencing genomes, though, so I don’t know how much the price has dropped over time.

I was a bit surprised to learn that the only mammals considered complete are mouse and man. There are 22 mammals with draft assemblies and another 26 in progress.  Will the draft genomes be completed any time soon?