Drosophila transcriptome reveals interesting patterns in splicing diversity

The Drosophila transcriptome has been analysed in different conditions to understand its diversity. Here, we re-analyse the supplementary materials(only from supplementary table-3) and find "interesting" patterns in splicing diversity.

 read.table(file="supp3",header=T)->S  
 #calculate percentage of genes with more than 1 isoform  
 length(S$transcripts[S$transcripts>1])*100/ length(S$transcripts)  
 #calculate percentage of genes with more than 1 protein  
 length(S$proteins[S$proteins>1])*100/ length(S$proteins)  
 summary(S)  

Over half of the 17,564 annotated genes (10,136; 57%) have more than one transcript isoform. However, only 37% (6,584) of the genes have more than one protein, suggesting that most of the transcript diversity is outside the protein coding regions. While the genes have on average 17.39 transcripts, only 2.59 proteins are found on average per gene. 

This contrast between the number of transcripts and proteins is very pronounced in genes with very large number of transcripts. Infact, the gene "gish" which has 18,972 transcripts has only 142 proteins. See below plot that shows the number of records with the words "exon" and "CDS" from the "gish" gene in the new annotation.

The above figure clearly shows that most "exon" records occur at the beginning and end of the gene. The "CDS" records are restricted to the core of the gene. The same pattern is better captured by the authors of the study using the new metric `per cent spliced in’ (Ψ) index.

Genome base pair correction tool

High rates of sequencing errors in Next-generation sequencing methods can percolate down into genome assemblies. Presence of incorrect bases in genome assemblies can have an impact on analysis that use the genome assembly for functional analysis.While programs like REAPR are able to correct genome assemblies based on mapped reads, it is focused on identifying large scale structural errors in genome assemblies. 

Draft genome assemblies are being generated for a large number of species using Next-generation sequencing methods. Here we present BCD (Base Correction Don), a tool that conservatively identifies single base-pair level errors in genome assembly sequences based on paired-end reads mapped to the genome.

 #!/usr/bin/perl  
 #use strict;  
 use warnings;  
 use List::Util qw(min max);  
 #perl fixbases.pl example.fasta example.coverage > example.corrected.fasta  
 # Input parameters  
 my $genome = $ARGV[0];  
 my $coverage = $ARGV[1];  
 my $totalcount=0;  
 my $minreadcoverage=20;# minimum number of reads required at site to consider it  
 my $maxreadcoverage=200;  
 my $maxindelcount=5;  
 open FASTA1, $genome or die $!;  
 open LOGS, ">test.bed" or die $!;  
 my $seqst_temp="";  
 my %seqs = ();  
 while($line = <FASTA1>){  
     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 FASTA1;  
 open FASTA2, $coverage or die $!;  
 while($line = <FASTA2>){  
 chomp $line;$line=$line . "\t}";  
 my @jelly=split(/\s+/,$line);  
 splice @jelly, 4, 0, 'SnowWhite', 'Humbert';  
 my $scaf=$jelly[0];  
 my $pos=$jelly[1];  
 my $oldbase=$jelly[2];  
 my $totalreads=$jelly[3];  
 my @jellyA=split(/\:/,$jelly[7]);#A  
 my @jellyC=split(/\:/,$jelly[8]);#C  
 my @jellyG=split(/\:/,$jelly[9]);#G  
 my @jellyT=split(/\:/,$jelly[10]);#T  
 my @jellyN=split(/\:/,$jelly[11]);#N  
 my $indelcount=0;  
 if(!($jelly[12]=~ /\}/)){  
 my @jellyindel=split(/\:/,$jelly[12]);#indel  
 $indelcount=$jellyindel[1];  
 }  
 my @counts;  
 push @counts,$jellyA[1];  
 push @counts,$jellyC[1];  
 push @counts,$jellyG[1];  
 push @counts,$jellyT[1];  
 push @counts,$jellyN[1];  
 push @counts,$indelcount;  
     if(($totalreads>$minreadcoverage)&&($indelcount<$maxindelcount)&&($totalreads<$maxreadcoverage)){$totalcount++;  
 my $maxbase= max @counts;  
 my $maxbasepercent=$maxbase*0.1;  
 my $maxbaseval=$oldbase;  
 if($jellyA[1]==$maxbase){$maxbaseval="A";}  
 elsif($jellyC[1]==$maxbase){$maxbaseval="C";}  
 elsif($jellyG[1]==$maxbase){$maxbaseval="G";}  
 elsif($jellyT[1]==$maxbase){$maxbaseval="T";}  
 elsif($jellyN[1]==$maxbase){$maxbaseval="N";}  
 else{$maxbaseval=$oldbase;}  
 my $maxbaseval2=lc $maxbaseval;  
 if((($oldbase=~ /[Aa]/)&&($jellyA[1]==0))||(($oldbase=~ /[Cc]/)&&($jellyC[1]==0))||(($oldbase=~ /[Gg]/)&&($jellyG[1]==0))||(($oldbase=~ /[Tt]/)&&($jellyT[1]==0))){  
 print LOGS "$scaf\t$pos\t$pos\t$maxbaseval\t$oldbase.\t$jellyA[1]\t$jellyC[1]\t$jellyG[1]\t$jellyT[1]\t$jellyN[1]\tzero\t$totalreads\n";  
 my $key=">" . $scaf;  
 my $prepos=$pos-1;#1 based positions  
 substr($seqs{$key},$prepos,1)= $maxbaseval;  
 }  
 elsif((($oldbase=~ /[Aa]/)&&($jellyA[1]<$maxbasepercent))||(($oldbase=~ /[Cc]/)&&($jellyC[1]<$maxbasepercent))||(($oldbase=~ /[Gg]/)&&($jellyG[1]<$maxbasepercent))||(($oldbase=~ /[Tt]/)&&($jellyT[1]<$maxbasepercent))){  
 print LOGS "$scaf\t$pos\t$pos\t$maxbaseval\t$oldbase.\t$jellyA[1]\t$jellyC[1]\t$jellyG[1]\t$jellyT[1]\t$jellyN[1]\tpercent\t$totalreads\n";  
 my $key=">" . $scaf;  
 my $prepos=$pos-1;#1 based positions  
 substr($seqs{$key},$prepos,1)= $maxbaseval;  
 }  
     }  
  }#end of while loop  
 close FASTA2;  
 my $iso;  
 foreach $iso (sort keys %seqs) {  
 print "$iso\n$seqs{$iso}\n";  
 }  
 #scaffold_0   50   c    45   GENOMEseq    {    =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 C:1:29.00:44.00:29.00:0:1:0.06:0.00:0.00:0:0.00:100.00:0.03 G:43:28.49:64.21:27.42:20:23:0.59:0.04:234.35:20:0.62:97.95:0.44    T:1:29.00:68.00:29.00:1:0:0.92:0.16:659.00:1:0.53:98.00:0.53  N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00   }  
 #scaffold_0   1    g    43   =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 C:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00   G:40:33.40:59.08:31.00:19:21:0.00:0.03:0.00:0:0.00:0.00:0.00  T:2:31.00:33.00:31.00:0:2:0.00:0.05:0.00:0:0.00:0.00:0.00    N:1:60.00:33.00:25.00:1:0:0.00:0.05:0.00:0:0.00:0.00:0.00  
 print LOGS "$totalcount\n";  
 close LOGS;  

All paired end data used for performing the genome assembly were mapped back to the genome assembly using bwa read mapper with default settings. Number of reads supporting each of the four bases at individual positions was tabulated using the program bam-readcount. The above script was then used to correct the genome for single base pair errors.

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  
 }