Hawaizaada - could the first human controlled powered airplane have taken flight on the beaches of Mumbai?

A recent movie Hawaizaada tries to think of an alternative history. One in which the very first human controlled powered airplane takes flight not on the beaches near Kitty Hawk but the one in Mumbai, India. Would it have been possible for such a character to even exist in pre-independence India? Most definitely, the example of CV Raman comes to mind. Was the environment conducive for innovation? Was it the startup capital :) The very presence of a few stars gives hope for the future.

The ridicule faced by the main character when he tries to convince the Maharaja of Baroda is probably not very different from what was experienced by the people who tried to claim that the book "Vaimānika Shāstra" was a scientific treatise. While it is now established by the leading scientists of today that the book is no more than a bunch of fantasies that have no grounding in reality, one has to at the least appreciate that somebody thought of concocting such brilliant lies. The movie Hawaizaada tries to give color to the fantasies found in the book and follows the life of Subbaraya Shastry who tried to make these models with help from Dr.Talpade (who has been replaced by a school dropout in the movie) of Bombay. In reality none of these models ever worked, obviously due to their mythical nature. 

While the movie is a commendable attempt at trying to create a "new" genre in Indian Cinema, it was a box office failure which received mixed criticism, most of it being negative. This seems more of a comment on the state & taste of the Indian entertainment industry than on the quality of the movie. One can only hope that with increasing literacy and as understanding of the Jaredness of the world increases, such dreamers will flourish. 

Have empires become larger over time?

The list of largest Empires at Wikipedia tries to provide an unbiased list of empires, their size and the approximate time during which these empires existed. Even though the very definition of the term "empire" is fraught with controversies, an operational definition that one might accept would involve large tracts of land, resources and people that have been brought under the control of some sort of central political entity. 

Obviously, over time the world has become more connected and has required larger regions to be under a single dominion to be considered an empire. In the distant past, the control of a collection of tribes might have been enough to be considered an empire. Overtime, this has grown to require the control of more than one continent. Thus, in a way the very way that empires have been defined, requires that they become larger over time. So, one would expect a positive correlation between the size of the empires and the time period in which they have existed. 

We use the list compiled by wikipedia from various scholarly sources to test this hypothesis. The list accessed on November 2015 had 218 empires with a mean maximum size of 1.161 square miles, with a minimum of 0.060 sq.miles (old kingdom Assyria in 1730 BC) and a maximum of 13 sq.miles (British Empire as measured in the year 1922). We find a significant positive correlation (pearson's correlation coefficient r= 0.29, p-value = 1.329e-05) between the size of empires and the "era" in which they happened reach their zenith ( also see figure below). The correlation is not driven by a few outliers and still shows a strong pattern when restricted to empires smaller than 0.5 million square miles (r=0.16) or 0.1 million square miles (r=0.39).

This suggests that empires have obviously become larger over time. It can also be seen that the largest of them all, the British Empire covered ~23% of the world's land area. One can wonder (when modelling the process is beyond..), if this is some sort of a limit on how big an empire can get before running into trouble. Keith Jeffrey in his book "The British Army and the Crisis of Empire, 1918-22: page 160", echoes this very thought expressed by 'Boney' Fuller when he said "we cannot increase the size of the army to fit the Empire; consequently there is only one thing we can do, namely, reduce the size of the Empire to fit our army".

The age of empires has interesting patterns that are potentially driven by complex processes. Would this correlation exist on a planet that completely lacked geographical barriers? How different would this pattern be if technology had not grown as fast as it has on our blue sphere?

How common is positive selection "inference" in different parts of the human genome?

The database of positive selection provides a handy resource for analyzing patterns of population specific positive selection inference trends and their prevalence in different parts of the human genome. After converting the database into a bed file (loosing some records during conversion makes it easier to run a bunch of analysis on the dataset. 

We next divide the genome into 50Kb chunks and count number of distinct populations in which positive selection was inferred in that window. By this method we can identify "hot-spots" of positive selection in the human genome. Below plot shows the number of populations in which positive selection inference has been done in each of these windows. Only one window on chromosome 2 stretching from 72500000 to 72550000 has been implicated in more than 20 populations. While the Y-chromosome and Mt are empty :) the X chromosome has ~2271 windows with atleast one population.

The above plot shows that it some regions of the genome have been implicated in multiple populations more often than others. Such regions could be termed "hot-spots" of positive selection. However, it reflects the bias in our ability to identify/infer such events rather than actual differences in the occurrence of such events.

Are some populations having positive selection inference spread out over the genome than others? Is it possible to make such a comparison with heterogeneity in the dataset? Methods used are different as well as the resolution of the scans. However, we see(in below figure) how the distribution of such inferences look for populations that have spanned more than 1000 of the ~62,000 windows (50 Kb) spanning the human genome.

MKK has the most number of windows followed by CHB, CEU and YRI. Since, not all observations are independent, it is hard to draw any broadscale conclusions from these numbers.

Titin exon specific GERP scores

Coming back to Titin gene, are all exons equally important? GERP scores could tell us something about how the exons vary in conservation level at base-pair resolution.

In the above figure, the y-axis is the exon rank and the x-axis has boxplots of GERP scores for each exon. Exons with mean GERP scores above 0 are colored red and those below zero are colored blue. One can clearly see the GERP score reduce considerably between exons 160 and 200.

HGT candidate genes - GERP score distribution

HGT candidate genes have been identified in vertebrate species based on gene absence/presence in metazoans. The actual workflow can be found in figure S4.

Interestingly, they identify 145 genes that are candidates for HGT into all primates from various taxon and another 51 from viruses. They provide Ensemble gene id's for these genes in the Human genome.

Of the 145 genes, only 135 could be found in the Ensemble biomart today (Ensemble 79). The following ten genes have been "retired" and are not used in below analysis. These genes seem to be randomly drawn from the different "confidence classes" used in the paper.

ENSG00000107618    class A HGT    Validated
ENSG00000157358    class A HGT    Validated
ENSG00000196333    class C HGT    Not-validated
ENSG00000204486    class C HGT    Validated
ENSG00000204502    class C HGT    Validated
ENSG00000204513    class C HGT    Validated
ENSG00000212857    class C HGT    Not-validated
ENSG00000256062    class A HGT    Validated
ENSG00000260383    class B HGT    Not-validated
ENSG00000263074    class B HGT    Validated

ENSG00000107618 has been merged with ENSG00000265203 which is also a HGT candidate. ENSG00000256062 has been merged into ENSG00000175164, this is the much talked about ABO blood group gene. ENSG00000260383 has been merged into ENSG00000179832 which also happens to be a HGT candidate. ENSG00000263074 is merged into ENSG00000141337 (arylsulfatase G).

A total of 4018 exons could be found for these 135 genes. GERP scores could not obtained for 8 of these genes [ENSG00000066813,ENSG00000183248,ENSG00000183549,ENSG00000183747,
ENSG00000204510,ENSG00000229571,ENSG00000232423,ENSG00000243073]. Here, we obtain the GERP scores for each base of all these exons. Example code from the web was updated to use a bed file as input and print the GERP score for each base within the bed intervals.

 #!/usr/bin/env perl  
 use strict;  
 use warnings;  
 use lib '~/biomart-perl/lib';  
 use Bio::EnsEMBL::Registry;  
 use Bio::EnsEMBL::Utils::Exception qw(throw);  
 my $reg = "Bio::EnsEMBL::Registry";  
 my $species = "Homo sapiens";  
 my $line = "";  
 open FILE1, $ARGV[0] or die $!;  
 $reg->load_registry_from_db(  
    -host => "ensembldb.ensembl.org",  
    -user => "anonymous",  
 );  
 #get method_link_species_set adaptor  
 my $mlss_adaptor = $reg->get_adaptor("Multi", "compara", "MethodLinkSpeciesSet");  
 my $mlss = $mlss_adaptor->fetch_by_method_link_type_species_set_name("GERP_CONSERVATION_SCORE", "mammals");  
 throw("Unable to find method_link_species_set") if (!defined($mlss));  
 my $slice_adaptor = $reg->get_adaptor($species, 'core', 'Slice');  
 throw("Registry configuration file has no data for connecting to <$species>") if (!$slice_adaptor);  
 my $cs_adaptor = $reg->get_adaptor("Multi", 'compara', 'ConservationScore');                 
 while($line=<FILE1>){  
 #$line=$ARGV[0];  
 chomp $line;  
 my @parts=split('\t',$line);  
 my $seq_region = $parts[0];  
 my $seq_region_start = $parts[1];  
 my $seq_region_end =  $parts[2];  
 my $slice = $slice_adaptor->fetch_by_region('toplevel', $seq_region, $seq_region_start, $seq_region_end);  
 throw("No Slice can be created with coordinates $seq_region:$seq_region_start-$seq_region_end") if (!$slice);  
 my $display_size = $slice->end - $slice->start + 1;   
 my $scores = $cs_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice($mlss, $slice, $display_size);  
 #print "number of scores " . @$scores . "\n";  
 foreach my $score (@$scores) {  
   if (defined $score->diff_score) {  
 #     printf("position %d observed %.4f expected %.4f difference %.4f\n", $score->position, $score->observed_score, $score->expected_score, $score->diff_score);  
      printf("$parts[3]\t$parts[4]\t%.4f\n",$score->diff_score);  
   }  
 }  
 }  
 close FILE1;  

Negative GERP scores suggest that the sites are probably evolving neutrally.

 read.table(file="GERP.scores",header=FALSE,stringsAsFactors=FALSE)->M  
 as.data.frame(aggregate(M,list(M$V2),mean))->N  
 as.data.frame(aggregate(M,list(M$V2),max))->P  
 N$Group.1[N$V3<0]->negative.mean.GERP  
 jpeg("GERP_HGT_distribution.jpeg")  
 boxplot(as.numeric(M$V3)~M$V2,outline=FALSE,col=ifelse(unique(M$V2) %in% negative.mean.GERP, "red", "blue"),xaxt="n",ylab="GERP score",main="GERP score by gene")  
 legend("bottomright", c("Negative mean GERP", "Positive mean GERP"), fill = c("red", "blue"))  
 dev.off()  

Despite 41 genes having a mean GERP score less than 0, each of these genes have atleast a few sites that are having positive GERP scores.

Based on above arguments, it could be suggested that 10 of the 145 genes are probably annotation artifacts. The 127 genes for which GERP scores could be obtained from whole genome alignments show moderate levels of conservation, suggesting these genes are probably not artifacts.

Lack of Isolation By Distance pattern in Tiger mitochondrial data

A recent paper, that used ancient DNA sequences from museum samples concluded that Bali, Javan and Sumatran tigers had a closer genetic relationship by comparing it with diagnostic mtDNA sequences from various other tiger subspecies.

Due to lack of precise geographic information about the sampling locations, despite having a matrix of genetic distances they could not look at isolation by distance patterns. The historical geographic ranges of these subspecies can be found on wikipedia. One could, just assume sampling locations within the geographic range and check if a pattern of IBD can be seen. While performing such assumptions in a paper might not be possible, a blog provides you some freedom to do so. However, somebody familiar with Tigers might be able to provide better assumptions. 

ALT or Siberian Tiger is assumed to have been sampled from Sikhote-Alin    45°N 136°E. VIR or Caspian Tiger from Gansu 38°N 102°E. AMO or South China Tiger from Qin mountains 33°N 107°E. COR or Indo-chinese Tiger from Yunnan 25°N 101°E. JAX or Malayan Tiger from Kelantan 5°N 102°E. SUM or Sumatran Tiger from Bukit Barisan Selatan National Park 5°S 104°E. SON or Javan Tiger from Java 7°S 110°E. BAL or Bali Tiger from Bali 8°S 115°E. TIG or Bengal Tiger from Hazaribagh National Park 24°N 85°E.

The R package sp provides many useful functions to deal with spatial data. We use the spDistsN1 function to get the using Euclidean or Great Circle distance between two co-ordinates.

 library(sp)  
 library(reshape)  
 library("ecodist")  
 someCoords <- data.frame(long=c(136,102,107,101,102,104,110,115,85), lat=c(45,38,33,25,5,5,7,8,24))  
 apply(someCoords, 1, function(eachPoint) spDistsN1(as.matrix(someCoords), eachPoint, longlat=TRUE))->X  
 X[upper.tri(X)]->G  
 pop1<-c("ALT","ALT","ALT","ALT","ALT","ALT","ALT","ALT","VIR","VIR","VIR","VIR","VIR","VIR","VIR","AMO","AMO","AMO","AMO","AMO","AMO","COR","COR","COR","COR","COR","JAX","JAX","JAX","JAX","SUM","SUM","SUM","SON","SON","BAL")  
 pop2<-c("VIR","AMO","COR","JAX","SUM","SON","BAL","TIG","AMO","COR","JAX","SUM","SON","BAL","TIG","COR","JAX","SUM","SON","BAL","TIG","JAX","SUM","SON","BAL","TIG","SUM","SON","BAL","TIG","SON","BAL","TIG","BAL","TIG","TIG")  
 data.frame(pop1,pop2,G)->Geodist  
 Geodist2 <- with(Geodist, G)  
 nams <- with(Geodist, unique(c(as.character(pop1), as.character(pop2))))  
 attributes(Geodist2) <- with(Geodist, list(Size = length(nams),Labels = nams,Diag = FALSE,Upper = FALSE,method = "user"))  
 class(Geodist2) <- "dist"  
 read.table(file="Tiger_Fst.txt",header=TRUE,sep="\t",row.names = 1)->M  
 na.omit(melt(M))->N  
 data.frame(pop1,pop2,N$value)->Gendist  
 Gendist2 <- with(Gendist, N$value)  
 nams <- with(Gendist, unique(c(as.character(pop1), as.character(pop2))))  
 attributes(Gendist2) <- with(Gendist, list(Size = length(nams),Labels = nams,Diag = FALSE,Upper = FALSE,method = "user"))  
 class(Gendist2) <- "dist"  
 mantel(Geodist2 ~ Gendist2, nperm=10000)  
   mantelr    pval1    pval2    pval3  llim.2.5% ulim.97.5%   
 -0.06313036 0.67800000 0.32230000 0.65900000 -0.25930292 0.12999925   
 jpeg("Tiger_isolation_by_distance.jpeg")  
 plot(G,N$value/(1-N$value),xlab="Distance in Km",ylab="Fst/(1-Fst) from mtDNA sequences",main="Isolation by Distance pattern",pch=16,col="blue",xlim=c(0,6000))  
 text(G,N$value/(1-N$value),labels=paste(pop1,"  ",pop2))  
 dev.off()  
 jpeg("Tiger_Fst_cladogram.jpeg")  
 hc = hclust(dist(Gendist2))  
 op = par(bg = "#DDE3CA")  
 plot(hc, col = "#487AA1", col.main = "#45ADA8", col.lab = "#7C8071", col.axis = "#F38630", lwd = 3, lty = 3, sub = "", hang = -1, axes = FALSE,xlab="",main="Tiger Fst cladogram")  
 axis(side = 2, at = seq(0, 400, 100), col = "#F38630", labels = FALSE,   lwd = 2)  
 mtext(seq(0, 400, 100), side = 2, at = seq(0, 400, 100), line = 1, col = "#A38630", las = 2)  
 dev.off()  

The mantel's test shows that the IBD pattern has a very weak negative correlation. This is not entirely unexpected, given the high values of Fst that are almost saturated. Would using more markers with greater resolution capture a IBD pattern? How strong are the bottlenecks?

 

They have already used the data to build a phylogenetic tree and a network. As i need to do something different, i build a cladogram using the Fst matrix. While this cladogram shows that the SUM and BAL are close to each other, it fails to group SON with them. Various other differences from the phylogenetic tree can be seen. 

Dog genome annotation quality has improved during the past 30 Ensembl releases

Quality of genome annotation changes with each new Ensembl release. How much does the quality actually change? 

We can quantify the amount of change in annotation quality by looking at specific parameters like number of genes, transcripts, exons etc. Here, the annotation of the dog genome from Ensembl release 48 to 78 is compared to understand the effect of change in genome assembly as well as role of annotation curation. 

First, we see the number of annotated genes(y-axis) in the Dog genome vs the Ensembl release(x-axis). It can be seen that the number of genes steadily  keeps increasing with each release. However, the sharp drop in number of genes in release 68 is surprising. The drop in number of genes from release 67 to 68 is less surprising as that release was accompanied by a change in genome assembly from BROADD2 to CanFam3.1

Various other features of the genome annotation such as Number of transcripts per Gene, Median coding length, Average Exons per gene and 3' UTR count change in a very different way (see below figure)

While such large changes in the genome annotation between releases might be scary, the engine of science chugs on. In the broader scheme of things, the fact that the Dog has 23,062 genes or 24,660 genes should hopefully not matter for the Dog breeders who seem to be driven by the discoveries made possible by the Dog genome and its annotation.