updates to v0.3.0

Author: Jon Palmer

bin/augustus_parallel.py

@@ -0,0 +1,105 @@
+#!/usr/bin/env python
+
+import sys, multiprocessing, subprocess, os, shutil, argparse, time
+from Bio import SeqIO
+
+#setup menu with argparse
+class MyFormatter(argparse.ArgumentDefaultsHelpFormatter):
+    def __init__(self,prog):
+        super(MyFormatter,self).__init__(prog,max_help_position=48)
+parser=argparse.ArgumentParser(prog='augustus_parallel.py', usage="%(prog)s [options] -i genome.fasta -s botrytis_cinera -o new_genome",
+    description='''Script that does it all...''',
+    epilog="""Written by Jon Palmer (2016) [email protected]""",
+    formatter_class = MyFormatter)
+parser.add_argument('-i','--input', required=True, help='Genome in FASTA format')
+parser.add_argument('-o','--out', required=True, help='Basename of output files')
+parser.add_argument('-s','--species', required=True, help='Augustus species name')
+parser.add_argument('--hints', help='Hints file (PE)')
+parser.add_argument('--cpus', default=2, type=int, help='Number of CPUs to run')
+args=parser.parse_args()
+
+#check for augustus installation
+try:
+    AUGUSTUS = os.environ["AUGUSTUS_CONFIG_PATH"]
+except KeyError:
+    if not args.AUGUSTUS_CONFIG_PATH:
+        print("$AUGUSTUS_CONFIG_PATH environmental variable not found, Augustus is not properly configured")
+        os._exit(1)
+if AUGUSTUS.endswith('config'):
+    AUGUSTUS_BASE = AUGUSTUS.replace('config', '')
+elif AUGUSTUS.endswith('config'+os.sep):
+    AUGUSTUS_BASE = AUGUSTUS.replace('config'+os.sep, '')
+
+#setup hints and extrinic input, hard coded for protein and transcript alignments from funannotate
+extrinsic = '--extrinsicCfgFile='+os.path.join(AUGUSTUS_BASE, 'config', 'extrinsic', 'extrinsic.E.XNT.cfg')
+
+def countGFFgenes(input):
+    count = 0
+    with open(input, 'rU') as f:
+        for line in f:
+            if "\tgene\t" in line:
+                count += 1
+    return count
+
+def runAugustus(Input):
+    FNULL = open(os.devnull, 'w')
+    species='--species='+args.species
+    hints_input = '--hintsfile='+os.path.join(tmpdir, Input+'.hints.gff')
+    aug_out = os.path.join(tmpdir, Input+'.augustus.gff3')   
+    with open(aug_out, 'w') as output:
+        if args.hints:
+            subprocess.call(['augustus', species, hints_input, extrinsic, '--gff3=on', os.path.join(tmpdir, Input+'.fa')], stdout = output, stderr= FNULL)
+        else:
+            subprocess.call(['augustus', species, '--gff3=on', os.path.join(tmpdir, Input+'.fa')], stdout = output, stderr = FNULL)
+
+
+#first step is to split input fasta file into individual files in tmp folder
+print("Splitting contigs and hints files")
+tmpdir = 'augustus_tmp_'+str(os.getpid())
+os.makedirs(tmpdir)
+scaffolds = []
+with open(args.input, 'rU') as InputFasta:
+    for record in SeqIO.parse(InputFasta, 'fasta'):
+        name = str(record.id)
+        scaffolds.append(name)
+        outputfile = os.path.join(tmpdir, name+'.fa')
+        with open(outputfile, 'w') as output:
+            SeqIO.write(record, output, 'fasta')
+
+#if hints file passed, split it up by scaffold
+if args.hints:
+    for i in scaffolds:
+        with open(os.path.join(tmpdir, i+'.hints.gff'), 'w') as output:
+            with open(args.hints, 'rU') as hintsfile:
+                for line in hintsfile:
+                    cols = line.split('\t')
+                    if cols[0] == i:
+                        output.write(line)
+
+#now loop through each scaffold running augustus
+if args.cpus > len(scaffolds):
+    num = len(scaffolds)
+else:
+    num = args.cpus
+print("Running augustus on %i scaffolds, using %i CPUs" % (len(scaffolds), num))
+p = multiprocessing.Pool(num)
+rs = p.map_async(runAugustus, scaffolds)
+p.close()
+while (True):
+    if (rs.ready()): break
+    remaining = rs._number_left
+    print "Waiting for", remaining, "augustus jobs to complete..."
+    time.sleep(30)
+print("Augustus prediction is finished, now concatenating results")
+with open(os.path.join(tmpdir, 'augustus_all.gff3'), 'w') as output:
+    for file in scaffolds:
+        file = os.path.join(tmpdir, file+'.augustus.gff3')
+        with open(file) as input:
+            output.write(input.read())
+
+join_script = os.path.join(AUGUSTUS_BASE, 'scripts', 'join_aug_pred.pl')
+with open(args.out, 'w') as finalout:
+    with open(os.path.join(tmpdir, 'augustus_all.gff3'), 'rU') as input:
+        subprocess.call([join_script],stdin = input, stdout = finalout)
+shutil.rmtree(tmpdir)
+print("Found %i total gene models" % countGFFgenes(args.out))

bin/funannotate-compare.py

@@ -20,7 +20,7 @@ def __init__(self,prog):
         super(MyFormatter,self).__init__(prog,max_help_position=48)
 parser=argparse.ArgumentParser(prog='funannotate-compare.py', usage="%(prog)s [options] genome1.gbk genome2.gbk",
     description='''Funannotate comparative genomics.''',
-    epilog="""Written by Jon Palmer (2015) [email protected]""",
+    epilog="""Written by Jon Palmer (2016) [email protected]""",
     formatter_class = MyFormatter)
 parser.add_argument('-i','--input', nargs='+', help='List of funannotate genome folders')
 parser.add_argument('-o','--out', default='funannotate_compare', help='Name of output folder')
@@ -30,6 +30,7 @@ def __init__(self,prog):
 parser.add_argument('--bootstrap', default=100, type=int, help='Number of bootstraps to run with RAxML')
 parser.add_argument('--num_orthos', default=500, type=int, help='Number of Single-copy orthologs to run with RAxML')
 parser.add_argument('--outgroup', help='Name of species for RAxML outgroup')
+parser.add_argument('--eggnog_db', default='fuNOG', help='EggNog database')
 args=parser.parse_args()
 
 
@@ -756,7 +757,7 @@ def __init__(self,prog):
             goList.append(description)
         goDict[col[0]] = goList
 
-EggNog = lib.eggnog2dict()
+EggNog = lib.eggnog2dict(os.path.join(parentdir, 'DB', args.eggnog_db+'.annotations.tsv'))
 iprDict = lib.dictFlipLookup(ipr, INTERPRO)
 pfamDict = lib.dictFlipLookup(pfam, PFAM)
 meropsDict = lib.dictFlip(merops)  

bin/funannotate-functional.py

@@ -17,7 +17,7 @@ def __init__(self,prog):
         super(MyFormatter,self).__init__(prog,max_help_position=48)
 parser=argparse.ArgumentParser(prog='funannotate-functional.py', usage="%(prog)s [options] -i genome.fasta -g genome.gff -o test -e [email protected]",
     description='''Script that adds functional annotation to a genome.''',
-    epilog="""Written by Jon Palmer (2015) [email protected]""",
+    epilog="""Written by Jon Palmer (2016) [email protected]""",
     formatter_class = MyFormatter)
 parser.add_argument('-i','--input', help='Folder from funannotate predict.')
 parser.add_argument('--genbank', help='Annotated genome in GenBank format')
@@ -36,6 +36,8 @@ def __init__(self,prog):
 parser.add_argument('--skip_iprscan', action='store_true', help='skip InterProScan remote query')
 parser.add_argument('--force', action='store_true', help='Over-write output folder')
 parser.add_argument('--AUGUSTUS_CONFIG_PATH', help='Path to Augustus config directory, $AUGUSTUS_CONFIG_PATH')
+parser.add_argument('--eggnog_db', default='fuNOG', help='EggNog database')
+parser.add_argument('--busco_db', default='fungi', choices=['fungi', 'metazoa', 'eukaryota', 'arthropoda', 'vertebrata'], help='BUSCO model database')
 args=parser.parse_args()
 
 def runIPRpython(Input):
@@ -113,6 +115,23 @@ def runIPRpython(Input):
         if args.skip_iprscan:
             lib.log.error("To run InterProScan you need to specify an email address to identify yourself to the online service")
             os._exit(1)
+            
+#check EggNog database, download if necessary.
+if not args.eggnog_db in lib.Nogs:
+    lib.log.error("%s is not a valid EggNog group, options are:\n%s" % (args.eggnog_db, ', '.join(lib.Nogs)))
+    os._exit(1)
+if not os.path.isfile(os.path.join(parentdir, 'DB', args.eggnog_db+'_4.5.hmm')):
+    lib.log.error("%s EggNog DB not found, trying to download and format..." % args.eggnog_db)
+    subprocess.call([os.path.join(parentdir, 'util', 'getEggNog.sh'), args.eggnog_db, os.path.join(parentdir, 'DB')], stdout=FNULL, stderr=FNULL)
+    if not os.path.isfile(os.path.join(parentdir, 'DB', args.eggnog_db+'_4.5.hmm')):
+        lib.log.error("Downloading failed, exiting")
+        os._exit(1)
+    else:
+        lib.log.error("%s downloaded and formatted, moving on." % args.eggnog_db)
+
+#check buscos, download if necessary
+if not os.path.isdir(os.path.join(parentdir, 'DB', args.busco_db)):
+    lib.download_buscos(args.busco_db)
 
 #need to do some checks here of the input
 if not args.input:
@@ -252,14 +271,15 @@ def runIPRpython(Input):
 eggnog_out = os.path.join(outputdir, 'annotate_misc', 'annotations.eggnog.txt')
 lib.log.info("Annotating proteins with EggNog 4.5 database")
 if not lib.checkannotations(eggnog_out):
-    lib.runEggNog(Proteins, args.cpus, 1e-10, os.path.join(outputdir, 'annotate_misc'), eggnog_out)
+    lib.runEggNog(Proteins, os.path.join(parentdir, 'DB', args.eggnog_db+'_4.5.hmm'), os.path.join(parentdir, 'DB', args.eggnog_db+'.annotations.tsv'), args.cpus, 1e-10, os.path.join(outputdir, 'annotate_misc'), eggnog_out)
 num_annotations = lib.line_count(eggnog_out)
 lib.log.info('{0:,}'.format(num_annotations) + ' annotations added')
 #run BUSCO OGS search
 busco_out = os.path.join(outputdir, 'annotate_misc', 'annotations.busco.txt')
-lib.log.info("Annotating proteins with BUSCO models")
+lib.log.info("Annotating proteins with BUSCO %s models" % args.busco_db)
+buscoDB = os.path.join(parentdir, 'DB', args.busco_db)
 if not lib.checkannotations(busco_out):
-    lib.runBUSCO(Proteins, args.cpus, os.path.join(outputdir, 'annotate_misc'), busco_out)
+    lib.runBUSCO(Proteins, buscoDB, args.cpus, os.path.join(outputdir, 'annotate_misc'), busco_out)
 num_annotations = lib.line_count(busco_out)
 lib.log.info('{0:,}'.format(num_annotations) + ' annotations added')
 #run signalP if installed, have to manually install, so test if exists first, then run it if it does