Saturday, January 11, 2014

VCF to calls format

Perl script that reads through a VCF file and prints it out in the "calls" format.

Although i could not find any proper documentation for the calls format, it can be described as diploid genotype calls represented using ambiguous DNA bases for each individual with the first line being the individual id's or names and the following lines being the genotypes at different positions.

 #!/usr/bin/perl  
 use List::Util qw[min max];  
 # Input parameters  
 my $vcf_file = $ARGV[0];  
 #Go through fasta file, extract sequences  
 open(IN, $vcf_file);  
 while($z=<IN>) {  
      if($z=~m/^##/){}  
      elsif($z=~m/^#/){  
           print "scaffold\tposition";  
           chomp $z;  
           @values=split(/\t/,$z);  
           $totind=scalar @values;  
           for($i=9;$i<$totind;$i++){  
           print "\t$values[$i]";  
           }  
      print "\n";  
      }  
      else{  
           chomp $z;  
           @values=split(/\t/,$z);  
           $refbase=$values[3];  
           $altbase=$values[4];  
                if(length $altbase < 2){#skip tri-allelic sites  
                     print "$values[0]\t$values[1]";  
                     $totind=scalar @values;  
                     for($i=9;$i<$totind;$i++){  
                     @parts=split(/\:/,$values[$i]);  
                $parts[0] =~ s/\///g;$parts[0] =~ s/0/$refbase/g;$parts[0] =~ s/1/$altbase/g;  
                $parts[0] =~ s/[A][A]/A/ig;$parts[0] =~ s/[T][T]/T/ig;$parts[0] =~ s/[C][C]/C/ig;$parts[0] =~ s/[G][G]/G/ig;$parts[0] =~ s/[.][.]/N/g;  
                $parts[0] =~ s/[AT][AT]/W/ig;$parts[0] =~ s/[GT][GT]/K/ig;$parts[0] =~ s/[AC][AC]/M/ig;$parts[0] =~ s/[AG][AG]/R/ig;$parts[0] =~ s/[CT][CT]/Y/ig;$parts[0] =~ s/[CG][CG]/S/ig;  
                     print "\t$parts[0]";  
                     }  
           print "\n";  
           }  
      }  
 }#end of file while  
 close IN;  
 #perl vcftocalls.pl genotyped.vcf > genotyped.calls  

Saturday, January 4, 2014

When The machine that won the war

The machine that won the war is a short story written by Issac Asminov in the 1960's. The story revolves around the end of a long war that has decimated an alien race leading to the survival of the human race.  One huge computer system is initially credited with making the decisions that lead to the victory. 

Data about results of various battle results, availability of resources and their locations are all collected and fed into this system. These data are processed by the super computer to provide accurate guidance to the commanders as to what decisions to make. However, over the course of the story we learn that the data being fed to it was too unreliable as most people involved in supplying the data could not be trusted. Moreover, even the data that was collected from these untrustworthy people was being manipulated so that it was "correct" before being sent to the computer. To make matters more interesting the computer system itself was not in a working conditions and could not trusted to interpret the data reliably. So the output was again being manipulated to take into account Murphy's law that says that anything that has to go wrong will go wrong. 

Finally it is revealed that although the data was being fixed at multiple stages, the final person responsible for using the decisions was using a toss of the coin to make the calls. 

Does this point to the fact that very complex systems are very hard to model, interpret and predict? So much so that a probabilistic approach performed as well as a very complicated model. May be the war could have been won much earlier if all the data was perfectly accurate and computer was in perfect working order and its instructions were followed to the letter. On the other hand having such perfect data without reliability issues might be hard to find in many systems that are very complex. Hence, the need to include a factor that takes the unreliability in the data into account. 

Sunday, December 22, 2013

Pseudogene distribution across the human genome

Pseudogene's are those genes which have lost their ability to code proteins or are not expressed due to other changes in the genome. The Ensemble genes 74 annotation of the human genome hg19 has 15,605 annotated pseudogenes. Based on extensive manual curation and automated predictions, the number of known pseudogenes has increased in the human genome over time. 

Pseudogene distribution in human
Distribution of pseudogenes across the human chromosomes (hg19)
Above figure (Ideograph generated by Idiographica ) shows the distribution of the pseudogenes annotated in the latest build of Ensemble. The complete lack of pseudogenes on the small arms of chromsome 13, 14, 15 and 22 is rather striking. A more comprehensive annotation dedicated to the identification and analysis of pseudo genes can be found at pseudogene.org. The latest build consists of 17172 records. The same pattern can be seen even in this more extensive annotation.


While the pattern is striking, it might be due to changes to the chromosome builds affecting the short arms of these chromosomes. However, the possibility of this being biological is of course very interesting. Could it correspond to chromatin type or some other genomic feature? Apparently it does correspond to the hetero-chromatic region of the genome that has not been sequenced.


Monday, November 11, 2013

White hair and assortative mating

 Assortative mating is defined as non-random within genotype or phenotype mating. Such non-random patterns of mating are of significance as they lead to population divergence and promote speciation. Due to its significance to population genetics as well as ecology it has been studied in a wide range of species for a traits such as age, behavior, phenology, genotype, ecotype, visual (color) differences among others. While most assortative mating is considered to be positive, few cases of negative assortative mating or diss-assortative mating have been recorded.

Various phenotypes like bird song, body shape, size or color have been seen to be the trait based on which assortment of mating preferences happen. More than one trait can also define mating preferences, as seen in the Northern Cardinalis. Depending on the nature of the trait that leads to assortative mating, consequences for the evolution of that trait are profound.

While ecologists, have no trouble defining assortative mating around phenotypes, geneticists seem to be very picky. Even the third edition of the book "Evolution" by Ridley defines assortative mating as "Tendency of like to mate with like. It can be for a certain genotype (e.g., individuals with genotype AA tend to mate with other individuals of genotype AA) or phenotype (e.g., tall individuals mate with other tall individuals)."

On the other hand in the book "An Introduction to Population Genetics" by Rasmus Nielsen and Montgomery Slatkin define assortative mating as "A mating structure in which pairs of individuals that are (genetically) similar to each other mate with higher probability than expected under random mating."

Is the "genetic" basis of the non-randomness implied or even necessary? Does the lack of a genetic basis make the phenomenon less interesting? The recent review/meta-analysis by Kirkpatrick's group seems to suggest that these and many more questions could help explain trends in the process of species formation.





Monday, June 10, 2013

Randomise order of lines in a file

Different programming languages are good at different things. R has many powerful statistical functions, while perl is good at data handling.

N random numbers from a certain range of numbers without re-sampling can be easily done in R with the "sample" function. To do the same thing in Perl, looping has (or atleast some form of iteration) to be used along with storing the results and checking to avoid re-sampling. 
 args<-commandArgs(TRUE)  
 totalsnps<-as.integer(args[1])  
 runumber<-as.integer(args[2])  
 sample(1:totalsnps,totalsnps,replace=F)->N  
 write.table(file=paste("rands.out",runumber,sep="."),N,col.names=F,row.names=F,quote=F)  

This Rscript can be run using the below line

 Rscript sampleit.r $linecount $iterationcount  

Once, the file with the new order of lines has been generated, it can be used by the perl script to write the file in the new order. We also keep the first two columns of the file unchanged and just randomise the remaining parts of the file.


 #!/usr/bin/perl  
 open RANDS, $ARGV[1] or die $!;  
 my %rhash = ();  
 my $count=1;  
#read in the file created by the R script in previous step
 while($lines = <RANDS>){  
 chomp $lines;  
 $rhash{$count}=$lines;  
 $count++;  
 }  
 close RANDS;  
#read the file that needs to be randomised and store it in hash with new order
 open STATS, $ARGV[0] or die $!;  
 my $hash = {CHR =>my $genename,POS =>my $pid,RESTATS =>my $pco1};  
 $mycount=1;  
 while($line = <STATS>){  
 chomp $line;  
 @tabs=split(/[ \t]+/,$line);  
 $hash{$mycount}{CHR}=$tabs[0];  
 $hash{$mycount}{POS}=$tabs[1];  
 $line =~ m/\w*\t\w*\t(.*)$/;  
 $hash{$rhash{$mycount}}{RESTATS}=$1;  
 $mycount++;  
 }#end of file while loop  
#check if number of lines match in old and new file
 if($mycount!=$count){print "mismatch in counts\n";}  
#print the file out in new order
 foreach $contigs (sort { $a <=> $b } keys %hash) {   
 print "$hash{$contigs}{CHR}\t$hash{$contigs}{POS}\t$hash{$contigs}{RESTATS}\n";   
 }  

This perl script just reads in the input file, stores it in a hash with new line order and then prints it out.

Friday, June 7, 2013

Distribution of the six most sighted birds

The folks over at the molecular Ecologist have a instructive tutorial on map making using R. The data from the birders in India, provides a novel dataset to try out these new map making skills. To start of with something simple, we look at the distribution of the 6 most sighted birds across India.

Only 6 species of birds have more than 700 recorded observations. The "Rosy Starling" leads the record with 864 observations followed by "Grey Wagtail" (857), "Common Sandpiper"(787),  "Barn Swallow" (752),"Pied Cuckoo" (750) and "Greenish Warbler" with 710 observations.

 read.csv("migrantwatch_reports.csv",row.names=NULL)->M  
 colnames(M)<-c("Species","Location name","City","State","Reporter","Date","Sighting type","Observation frequency","Start date","On behalf of","Latitude","Longitude")  
 numberofclusters<-10  
 N<-data.frame(M$Latitude,M$Longitude)  
 kmeans(N, numberofclusters)->P  
 M$cluster<-P$cluster  
 as.data.frame(P$centers)$M.Latitude->lat  
 as.data.frame(P$centers)$M.Longitude->lon  
 BSwallow<-as.data.frame(table(M$cluster[M$Species=="Barn Swallow"]))$Freq  
 CSandpiper<-as.data.frame(table(M$cluster[M$Species=="Common Sandpiper"]))$Freq  
 GWarbler<-as.data.frame(table(M$cluster[M$Species=="Greenish Warbler"]))$Freq  
 GWagtail<-as.data.frame(table(M$cluster[M$Species=="Grey Wagtail"]))$Freq  
 PCuckoo<-as.data.frame(table(M$cluster[M$Species=="Pied Cuckoo"]))$Freq  
 RStarling<-as.data.frame(table(M$cluster[M$Species=="Rosy Starling"]))$Freq  
 library(maps)  
 library(mapdata)  
 library(mapplots)  
 jpeg("top6.map.jpeg")  
 map("worldHires","India",col="gray90", fill=TRUE)  
 for (i in 1:numberofclusters) {  
 add.pie(z=c(BSwallow[i],CSandpiper[i],GWarbler[i],GWagtail[i],PCuckoo[i],RStarling[i]), x=lon[i], y=lat[i], col=c("blue","brown","Green","Grey","black","yellow"), labels="")  
 }  
 legend("topright",c("Barn Swallow","Common Sandpiper","Greenish Warbler","Grey Wagtail","Pied Cuckoo","Rosy Starling"),col=c("blue","brown","Green","Grey","black","yellow"),pch="*")  

Running the above code gives a nice(take up the politics with the guys who made the R package) map with pie-charts that looks below image:

The pie-charts are generated for the mean locations that are identified by k-means clustering on all observations with a K of 10. So, this also gives us an idea of where most of the observations come from.

Rosy Starling's are definitely observed more often on the west coast while the Grey Wagtail is seen more often to the south. A large number (337) of the Rosy Starling observations come from Shantilal Varu, whose favorite birds are waders. 


Thursday, June 6, 2013

Importance of recording observation frequency


With some feedback from MigrantWatch regarding the post about weekend bias in bird observations, it seems that the "Observation frequency" field is very important indeed. The whole weekend effect  does seem to disappear when the data is filtered for only Daily observations.



While, only 3029 of the 22622 observations ~13% of the observations are Daily, the weekend effect does have a nice control and should affect first and last sighting results much less. Just staring at the graphs might not be appealing to the more statistically minded, so doing a test like analysis of variance might be more productive.


 > data.frame(Day=format(as.Date(M$Date), format="%A"),Year=format(as.Date(M$Date), format="%Y"))->WD  
 > as.data.frame(table(WD))->WDD  
 > WDD[WDD$Year %in% c(2007:2013),]->WDDF  
 > aov(Freq~Day,WDDF)->WDDFA  
 > summary(WDDFA)  
       Df Sum Sq Mean Sq F value  Pr(>F)    
 Day     6 3875954 645992 5.0784 0.0005365 ***  
 Residuals  42 5342556 127204             
 ---  
 Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   
 > M[M$"Observation frequency" == "Daily",]->MD  
 >   
 > data.frame(Day=format(as.Date(MD$Date), format="%A"),Year=format(as.Date(MD$Date), format="%Y"))->DWD  
 > as.data.frame(table(DWD))->DWDD  
 > DWDD[DWDD$Year %in% c(2007:2013),]->DWDDF  
 > aov(Freq~Day,DWDDF)->DWDDFA  
 > summary(DWDDFA)  
       Df Sum Sq Mean Sq F value Pr(>F)  
 Day     6  466  77.75 0.0548 0.9993  
 Residuals  42 59625 1419.64          

While the ANOVA shows a significant effect of Day for the whole dataset, it does not show one for data filtered for only daily observations. A much larger sample size would be needed to ensure that this field is being used in a proper way.

The other point about these observations being affected by Holidays also seems very valid. 
 as.data.frame(table(data.frame(Day=format(as.Date(M$Date), format="%j"))))->DY  
 > DY[DY$Freq>200,]  
   Var1 Freq  
 20  020 246  
 27  027 224  
 337 337 261  
 358 358 213  
 365 365 204  
By tabulating the observations by day of year, annual holidays like Christmas and New years should be captured. Lo and behold! Christmas (358th day of the year) and New years eve (365th day) are in the top five.

Why does 337th day (Lawyers day in India) have the most observations? May be most of the birders are lawyers!!