Extract sequences from multi fasta file in a particular order

 #!/usr/bin/perl  
 use warnings;  
 #perl extractInOrder.pl test.fasta test.list  
 # Input parameters  
 open FASTA, $ARGV[0] or die $!;  
 my $seqst_temp="";  
 my $seqs = ();  
 while($line = <FASTA>){  
 if($line=~ /^>/){  
 if($header){  
 $header=~s/>//;  
 $seqs{$header}=$seqst_temp;  
 }  
 chomp $line;  
 $header="";  
 $header=$line;  
 $seqst_temp="";  
 }  
 else{  
 $line =~ s/[\n\t\f\r_0-9\s]//g;  
 $seqst_temp .= $line;  
 }  
 }#end of while loop  
 if($header){  
 $header=~s/>//;  
 $seqs{$header}=$seqst_temp;  
 }  
 close FASTA;  
 open FASTA, $ARGV[1] or die $!;  
 while($line = <FASTA>){  
 chomp $line;  
 $line=~s/>//;  
 if(exists $seqs{$line}){print ">$line\n$seqs{$line}\n";}  
 else {  
 print "Sequence header $line does not exist in fasta\n";  
 exit();  
 }  
 }  
 close FASTA;  

Merge multiple kmer count hashes into one

A previous attempt at merging two kmer count hashes was neither memory efficient nor capable of merging multiple kmer count hashes. Here, we use the age old trick of sorting to write a more memory efficient script that can handle as N number of hashes.

 cat *_"$kmer"_counts.fa|sort > sorted_"$kmer"_all.fa  

The above command will concatenate all the hashes and sort it. This sorted file can then be used by the below perl script to merge the hashes. Since, all kmers that need to be merged are in adjacent lines, the memory needed for merging is drastically reduced compared to the previous script.

 #!/usr/bin/perl  
 use warnings;  
 # Input parameters  
 open FASTA1, $ARGV[0] or die $!;  
 my $previous="Kmer";  
 my $previousCount="Kmercount"; 

 my @jelly;

  while($line = <FASTA1>){  
 chomp $line;  
 @jelly=split(/\s+/,$line);  
      if($previous=~/$jelly[0]/){  
      $previousCount=$previousCount+$jelly[1];  
      }  
      else{  
      print "$previous\t$previousCount\n";  
      $previous=$jelly[0];$previousCount=$jelly[1];  
      }  
 }  
      #printing last line if it needed merging  
      if($previous=~/$jelly[0]/){  
      print "$previous\t$previousCount\n";  
      }  
 close FASTA1;  

Merge two kmer count hashes into one

 #!/usr/bin/perl  
 use warnings;  
 # Input parameters  
 open FASTA1, $ARGV[0] or die $!;  
 open FASTA2, $ARGV[1] or die $!;  
 my %seqs = ();  
 while($line = <FASTA1>){  
 chomp $line;  
 my @jelly=split(/\s+/,$line);  
 $seqs{$jelly[0]}=$jelly[1];  
 }  
 close FASTA1;  
 while($line = <FASTA2>){  
 chomp $line;  
 my @jelly=split(/\s+/,$line);  
 if(exists $seqs{$jelly[0]}){$seqs{$jelly[0]}=$seqs{$jelly[0]}+$jelly[1];}  
 else{$seqs{$jelly[0]}=$jelly[1];}  
 }  
 close FASTA2;  
 foreach $iso (sort keys %seqs) {  
 print "$iso\t$seqs{$iso}\n";  
 }  

Will take two lists of kmer counts and merge them into one. A kmer count list consists of two columns. The first column being the kmer itself and the second being its count.

This might be useful in merging the output of a program like jellyfish after running it on each chromosome separately. While jellyfish has a merge function, it requires the hashes to be of equal size.

May be this can be re-written as a perl one liner...

Ka ka ka ka ka ka ka ka

Previous crow stuff has been covered by Jerry Coyne himself. This in turn stimulates lot of debate on what constitutes a species. Fortunately, most people agree that understanding the genetics is more important and relevant than fighting over the definition of species. 

Like the title of this post, the call of a crow can be written as "ka ka ka ka ka" or "kaw kaw kaw kaw" or "Caw caw caw caw". They are all trying to capture the same thing. As long as the underlying genetics is well understood, it should be ok to not worry about the definition for now. Hopefully, over time our understanding of genetics will be so good that we may even come up with a single definition that most of us can agree on.

Get overlapping sequences from multifasta file

 #!/usr/bin/perl  
 #use strict;  
 use warnings;  
 # Input parameters  
 open FASTA, $ARGV[0] or die $!;  
 my $seqst_temp="";  
 my %seqs = ();  
 my $iso="";  
 my $maxlen=0;  
 my $maxval="";  
 while($line = <FASTA>){  
 if($line=~ /^>/){  
 if($header){  
 $seqs{$header}=$seqst_temp;  
 }  
 chomp $line;  
 $header="";  
 $header=$line;  
 $seqst_temp="";  
 }  
 else{  
 $line =~ s/[\n\t\f\r_0-9\s]//g;  
 $seqst_temp .= $line;  
 }  
 }#end of while loop  
 if($header){  
 $seqs{$header}=$seqst_temp;  
 }  
 close FASTA;  
 $maxlen=0;  
 foreach $iso (sort keys %seqs) {  
 my $line1=$iso;  
 my $line2=$seqs{$iso};  
 my $flag=0,$overlap=500,$length=1000;  
 $seqlen=length $line2;  
 while(($seqlen-$flag)>$length){  
 if(($seqlen-($flag+$overlap))<$length){  
 $length=$seqlen-$flag;  
 }  
 $nextseq=substr $line2,$flag,$length;  
 print $line1.":".$flag."-".($flag+$length)."\n";  
 print $nextseq."\n";  
 $flag=$flag+$overlap;  
 }#end of seqlen while loop  
 }  

As the crow forks

We had our work on crows (The genomic landscape underlying phenotypic integrity in the face of gene flow in crows) published in the journal Science. Very interesting interpretations of the work ranging from it being proof for evolution to reasoning for human assortative mating in a political context have been presented. It has also been covered by the Swedish Radio, Der Spiegel (a popular German weekly) and few other offline sources.

A well written description of the science that puts the work in a broader context is also available.[It should also be noted that a description of another paper about the human polymorphic inversion 17q21 from the same author was a useful resource while looking at the ancestral state of our fixed differences.]

The whole genome was scanned for differences in DNA sequence content between carrion and hooded crows using 60 crows. While most of the genome is the same, small regions of the genome do show differences in DNA sequence and gene expression. These differences in DNA sequence are localized to specific regions of the genome. Infact all 81 of 82 fixed differences that we find are within a 2 megabase region. 

Snapshots of a region of the genome with fixed differences seen between hooded and carrion crow sequencing reads provides a very visual idea of how these fixed differences "look". Each of the boxes in the figure is a sequencing read (from the bam file) mapped to the reference genome. Positions which differ from the reference genome have the actual base that differs from the reference in each read.

Below figure shows the sequencing reads from one crow from Poland. Note that all the positions are similar to the reference genome. All 30 hooded crow individuals (from Sweden and Poland) look exactly like this at this base.

The next figure shows the sequencing reads from one crow from Germany at the same position as the hooded crow above. However, the blue box has "C" written in all its reads. At this position the Carrion crow has a base that is different from the Hooded crow reference genome. All the 30 carrion crows (from Spain and Germany) have this "C" base at this position.

Similar to the above position, we identify 81 other sites that show such a "fixed difference". The rest of the genome seems to be extremely similar between the carrion and hooded crows.

The Elephant Catchers

It has been 6 years since i read and blogged about Subroto Bagchi's book Go kiss the world. As with the previous book, finished reading the book "The Elephant Catchers" in one sitting. The writing style is very captivating and keeps you reading. 

The Elephant Catchers is all about growing a company that is hunting rabbits (small game) into one that now catches elephants. While many people start a startup, very few can make it into a self-sustaining and growing company. The book deals with the mistakes that happen and how to go about "fixing" those mistakes. 

Many anecdotes from different points along the way makes for an interesting read. The use of personification at different points to explain very different stories is very appealing. For example, at one point the company wins a big client who could be persuaded to hire a dog that is loyal instead of cat that is pleasing. 

Dividing the book into 6 parts with sub-chapters keeps it well organized. All of them have one thing in common, it is about scaling things. Be it scaling ideas or infrastructure the wisdom of having been there and having done that shines through brilliantly.