Mash

Introduction

MinHash is a general dimensionality-reduction technique and it is used by Mash to reduce large sequences and sequence sets to small, representative sketches with the result that global mutation distances (Mash distances) can be rapidly estimated.

Other aspects

• terms itself as an alignment-free method

Usage

Typical analysis

Mash is run on genomes. These will usually be de-novo assembled genomes from tools such as Velvet or SPAdes.

Parallel Usage on gridengine

We'll go through a process here of running Mash on a set of samples, using the DRMAA library to launch Gridengine job arrays.

The scripts will take as argument a file listing of the sample names, and it is assumed there are two pair-ended FASTQ reads per sample. It is also assumed that the paired-ended samples appeared in ordered fashoin in the file-listing: i.e. each consecutive set of two lines represent one sample.

De novo assembly

We use SPAdes with the --meta option here as we are dealing with metagenomes. First we need the python DRMAA script which will control the where and how the job script will be run. The keys issues for it will be:

• what job script to run (this need only be a bash script, BUT it must have executable permissions)
• what is the name of the paired up FASTQ filelisting.
• With how many threads/processes which EACH sample paired will run.
• The queue name is not mentioned, so it will launch on all.q.

#!/usr/bin/env python2.7
import os, sys, drmaa

def main():
    """Submit an array job."""
    argquan=len(sys.argv)
    if argquan != 4:
        print "This script requires two arguments: 1) the script to run in ja mode  2) filelist of absolute paths and filenames 3) Number of threads/CPU for *each* job array"
        sys.exit(2)

    s = drmaa.Session()
    s.initialize()
    print 'Creating job template'
    jt = s.createJobTemplate()
    jt.workingDirectory=os.getcwd() # means sge job output will be deposited here.
    jt.remoteCommand = jt.workingDirectory + '/' +sys.argv[1]

    with open(sys.argv[2]) as x: fl = x.read().splitlines()
    eflsz=len(fl)
    PL=[]
    for i in xrange(eflsz):
        PL.append(fl[i])
    pld2=len(PL)/2
    jt.args =PL 
    nm='-N jadrm0'
    # this is an intensive IO job, don't want to whack the FS too much
    jt.joinFiles=True
    jobid = s.runBulkJobs(jt, 1, pld2, 1)
    print 'Your job has been submitted with id ' + str(jobid)
    print 'Cleaning up'
    s.deleteJobTemplate(jt)
    s.exit()
       
if __name__=='__main__':
    main()

Note the lines:

pld2=len(PL)/2
jt.args =PL 

This takes the file listing and sends out the names in pairs as argument to the jobscript. The number of jobs is half the total number of lines, obviously, as a pair represents one sample.

Next we require the job script itself.

#!/bin/bash
ARGA=( $@ )
TRUEARRIDX1=$((2*$SGE_TASK_ID-2)) # be careful, easy to be confused with this.
TRUEARRIDX2=$(($TRUEARRIDX1+1))
SAMP1=${ARGA[$TRUEARRIDX1]}
SAMP2=${ARGA[$TRUEARRIDX2]}
STRPLDSLASHES=${SAMP1##*/} # gets rid of leading ./ or / however, there could be trouble if there isn't
PRFX=${STRPLDSLASHES%%_*}
OUTDIR=$PRFX
if [ -d $OUTDIR ]; then
        echo "Good $OUTDIR exists, no need to create it"
    else
       echo "Ooops $OUTDIR does not exist. Need to run mkdir $OUTDIR"
       mkdir $OUTDIR
fi
module load SPAdes/3.10.0
spades.py --meta -o $OUTDIR -t $NSLOTS --pe1-1 $SAMP1 --pe1-2 $SAMP2

This loads the right module and tells SPAdes to output its results into a directory representing the sample. Beware the filename parsing that goes on here, it is highly dependent on the way the FASTQ filenames are formulated. It will probably need to be edited for your own filename schema.

The output assembly files will all be given the same name scaffolds.fasta, in a subfolder corresponding (in this case) to a unique sample root name.

The following is a helper script which wraps the two above. Obviously this is entirely specific to the filenames chosen:

#!/bin/bash
./jadrm_spades_e2.py jas_spades_e2.sh fq.lst 4

• the DRMAA python script comes first.
• the shell jobscript comes second
• the filelisting of paired-up FASTQ files comes third.
• the individual parallelisation parameter comes last.

This will launch the jobs in parallel and the gridengine command qstat can be used to monitor them.

Mash step on the assemblies

Bear in mind that Mash's most common usage is in pair-wise comparisons, so each sample must be compared with all the others. For n samples therefore, there will be n*(n-1)/2 comparisons. These can all be run in parallel, but unfortunately the number of operations goes up very quickly, so one needs to be wary, For example 45 samples which can be considered only a medium sized payload, is 990 pairwise comparisons which is five times the number of processors on marvin.

Mash has preprocessing step called "sketch" which can be done first. Note that this is a simple one-for-one operation, as it just converts each assembly into a sketch, or a .msh' file, What is now necessary is a new filelisting with the absolute paths to the sample scaffolds.fasta files.

With that filelisting in hand, we need the DRMAA python script. For a one-to-one operation like this one, all we need is a DRMAA script in its simplest form:

#!/usr/bin/env python2.7
import os, sys, drmaa
def main():
    """Submit an array job """
    argquan=len(sys.argv)
    if argquan != 3:
        print "This script requires two arguments: 1) the script to run in ja mode  2) filelist of absolute paths and filenames"
        sys.exit(2)
    s = drmaa.Session()
    s.initialize()
    print 'Creating job template'
    jt = s.createJobTemplate()
    jt.workingDirectory=os.getcwd() # means sge job output will be deposited here.
    jt.remoteCommand = jt.workingDirectory + '/' +sys.argv[1]
    with open(sys.argv[2]) as x: fl = x.read().splitlines()
    eflsz=len(fl)
    jt.args =fl
    nm='-N ja_mashske0'
    jt.nativeSpecification='-V -q all.q,highmemory.q ' +nm
    jt.joinFiles=True
    jobid = s.runBulkJobs(jt, 1, eflsz, 1)
    print 'Your job has been submitted with id ' + str(jobid)
    print 'Cleaning up'
    s.deleteJobTemplate(jt)
    s.exit()
   
if __name__=='__main__':
    main()

This python script will only need the jobscript and filelisting as arguments. The job script is as follows:

#!/bin/bash
# pairwise treatment: making the sketches, one per genome
ARGA=( $@ )
# ARGA=( $(cat f.l) )
ARGSZ=${#ARGA[@]}
GENOME=${ARGA[$(($SGE_TASK_ID-1))]} # because bash array "ARGA" is zero indexed while SGE TASK IDs are 1-indexed

module load Mash
mash sketch ${GENOME}

The effect of this will be to create a .msh file in the same directory as the scaffolds.fasta file.

Next up is the main Mash process.

• The DRMAA file needs to handle the pairing up of each assembly with every other assembly.
• Mash's output is in fact a mere one line at the end of which is a fraction, this is the Mash distance between the two assemblies. Each of these comparisons will be held in their own file.

Here is the DRMAA script

#!/usr/bin/env python2.7
import os, sys, drmaa
def main():
    """Submit an array job"""
    argquan=len(sys.argv)
    if argquan != 3:
        print "This script requires two arguments: 1) the script to run in ja mode  2) filelist of absolute paths and filenames"
        sys.exit(2)
    s = drmaa.Session()
    s.initialize()
    print 'Creating job template'
    jt = s.createJobTemplate()
    jt.workingDirectory=os.getcwd() # means sge job output will be deposited here.
    jt.remoteCommand = jt.workingDirectory + '/' +sys.argv[1]
    with open(sys.argv[2]) as x: fl = x.read().splitlines()
    eflsz=len(fl)
    PL=[]
    for i in xrange(0,eflsz,2):
        for j in xrange(i+2, eflsz, 2):
            PL.append(fl[i])
            PL.append(fl[j])
    pld2=len(PL)/2
    jt.args =PL
    nm='-N pairing'
    jt.nativeSpecification='-V -q all.q,highmemory.q ' +nm
    jt.joinFiles=True
    jobid = s.runBulkJobs(jt, 1, pld2, 1)
    print 'Your job has been submitted with id ' + str(jobid)
    print 'Cleaning up'
    s.deleteJobTemplate(jt)
    s.exit()
   
if __name__=='__main__':
    main()

Links

• main Mash documentation