Mash perl module

Hey y'all, I made a new perl module to read and write Mash files.  I made it a complete package, adding it to CPAN with documentation and adding unit testing with continuous integration on Travis-CI.  Most of the back end relies on the binary, but I re-implemented the mash distance formula in Perl.

Why?  I am the developer for Mashtree and it is becoming increasingly difficult to parse Mash files, and it became evident that I needed some object-oriented programming.  Because it might be helpful to the community, I made it an actual perl module.

One difficult part of the process was not being able to read the cap'n proto formatted mash files.  This is due to no relevant perl parser that was available for it.  Therefore, I looked for a different way to serialize mash sketches.  Ultimately because mash info displays JSON, I went with JSON and gzipped JSON.

Here is a quick synopsis: 
use strict;
use warnings;
use Mash;

# Sketch all fastq files into one mash file.
system("mash sketch *.fastq.gz > all.msh");
die if $?;
# Read the mash file.
my $msh = Mash->new("all.msh"); 
# All-vs-all distances 
my $distHash = $msh->dist($msh); 
# Or, one liner perl -MMash -e '$msh1 = Mash->new("1.msh"); $msh2 = Mash->new("2.msh"); $dist = $msh1->dist($msh2);

I ran a preliminary all-vs-all test of 303 genomes.  To make the file, I used Mash v2.  To run all-vs-all distances, I ran mash dist vs $msh->dist().  In this test, the binary was 15x faster (1.9s vs 28s).  One more benchmark shows that the resulting json.gz files are approximately the same size as the original msh files.
Let me know what you think!  As of the time of this post, the version is 0.5.
https://metacpan.org/pod/Mash
https://github.com/lskatz/perl-mash
If you would like to contribute, my wish list includes

• Speed optimization
• Mash paste
• Native reading of mash files
• Native writing of mash files

Sampling the taxonomy database

I was a little frustrated that every time I wanted to try out my new Bio::DB::Taxonomy-based script, it would take a few minutes to run.... and then I would find the bug in my script, fix it, and run it again.  I couldn't find a small example database to run my script on, and so I created one.  This is the NCBI taxonomy, spliced for just Firmicutes (taxid: 1239)  Enjoy!

# Filter the nodes file.
# I used a recursive function printChildren to
# print the taxonomy lines.
perl -F'\t\|\t' -MData::Dumper -lane '
  BEGIN{
    sub printChildren{
      my $parent=shift;
      return if(!$child{$parent});
      for my $child(values($child{$parent})){
        print $nodes{$child};
        printChildren($child);
      }
    }
  }
  push(@{$child{$F[1]}},$F[0]);
  $nodes{$F[0]}.=$_;
  END{printChildren(1239);
}' < nodes.dmp > exampleNodes.dmp

# Backtrack and filter the names file.
perl -F'\t\|\t' -Mautodie -lane '
  BEGIN{
    open($fh, "names.dmp");
    while(<$fh>){
      my($nodeID)=split(/\t\|\t/);
      $names{$nodeID}.=$_;
    }
    close $fh;
    print STDERR "Indexed!  Searching and printing.";
  }
  chomp $names{$F[0]};
  print $names{$F[0]};
' < exampleNodes.dmp > exampleNames.dmp

Show that the taxonomy has been significantly filtered to one or two orders of magnitude.

$ wc -l *.dmp
   107700 exampleNames.dmp
    79466 exampleNodes.dmp
  2401017 names.dmp
  1614627 nodes.dmp
  4202810 total


Then, loading the database in BioPerl:
use Bio::DB::Taxonomy;
$db = Bio::DB::Taxonomy->new(-source=>"flatfile", -nodesfile=>"exampleNodes.dmp", -namesfile=>"exampleNames.dmp");

Custom Kraken database

Like many labs around the world, we use Kraken for contamination detection.  This isn't its intended purpose however because it is supposed to be a taxonomic classifier for metagenomics samples.

Despite that, it has been incredibly useful.  We just want to make it more useful.  Many weird species, especially those without genome assemblies on NCBI, show up really weird in Kraken.  So when look for contamination in a weird species, we get random hits from random species or even a majority of unclassified reads.  On the other hand, if we run Kraken on a standard Escherichia coli genome, we get near-perfect results.

Working with our subject matter experts (SMEs), I was able to piece together a preliminary set of completed genome assemblies where the risk of their being contaminated was virtually none.  Then, I tried making a custom Kraken database.  The instructions are straightforward... until they're not.  There were some unsaid things such as the need to erase the description on sequence headers.  Below is my workflow for creaking a totally new Kraken database and not simply adding onto a pre-existing one.

# First, make a directory of fasta files 
# whose extensions are .fasta. The specific
# extension doesn't matter, but it helps with
# the find command later.

# Display the assemblies directory with tree.
# This is a for-example, not necessarily what I have.
# Each subfolder should be formatted Genus_species.
$ tree genomeAssemblies/* 
Escherichia_coli
├──Escherichia_coli_44.fasta
└──Escherichia_coli_45.fasta
Salmonella_bongori
├── NCTC-12419.fasta
└── N268-08.fastaSalmonella_enterica
└── Salmonella_enterica_90.fasta

# Download the taxonomy and initialize the Kraken
# database.
$ kraken-build --download-taxonomy --db custom-kraken

# Alter the Shigella species so that they are subspecies
# to E. coli
cd custom-kraken/taxonomy
mv nodes.dmp nodes.bak
# Anything with Genus Shigella (taxid: 620), 
# change to parent E. coli (taxid: 562)
perl -Mautodie -e 'open($fh, "nodes.bak"); while(<$fh>){@F=split /\t\|\t/; if($F[1]==620){ $F[1]=562; $_=join("\t|\t", @F); } print; } close $fh;' > nodes.dmp
cd -

# Next format the fasta files using a custom script
# which can be found in my misc github repo
# https://github.com/lskatz/lskScripts/blob/master/formatFastaForKraken.pl 
# This script adds a kraken tag and removes the 
# fasta description. Modifications are in-place.
$ ls genomeAssemblies | xargs -P 12 -n 1 bash -c '
  genusSpecies=$(sed "s/_/ /g" <<< $0); 
  taxid=$(grep "|"$'\t'"$genusSpecies"$'\t'"|" custom-kraken/taxonomy/names.dmp | cut -f 1 | head -n 1)
  formatFastaForKraken.pl --taxid=$taxid $0/*.fasta
'

# This automation with find/xargs only works down to 
# the species level. For any deeper diving, you can 
# use the custom script to replace the taxid with 
# something more specific.
$ formatFastaForKraken.pl --taxid=41520 genomeAssemblies/Salmonella_enterica/Salmonella_enterica_90.fasta

# Add genome assemblies to the custom database
find . -name '*.fasta' | \
  xargs -P 12 -n 1 kraken-build --db custom-kraken --add-to-library

# Build the database
kraken-build --build --threads 12 custom-kraken