updates to v0.2.2

Author: Jon Palmer

bin/funannotate-compare.py

@@ -107,6 +107,7 @@ def __init__(self,prog):
 busco = []
 gbkfilenames = []
 scinames = []
+signalp = []
 num_input = len(args.input)
 if num_input == 0:
     lib.log.error("Error, you did not specify an input, -i")
@@ -140,6 +141,7 @@ def __init__(self,prog):
         pfam.append(lib.getStatsfromDbxref(GBK, 'PFAM'))
         cazy.append(lib.getStatsfromNote(GBK, 'CAZy'))
         busco.append(lib.getStatsfromNote(GBK, 'BUSCO'))
+        signalp.append(lib.getStatsfromNote(GBK, 'SECRETED'))
         lib.parseGOterms(GBK, go_folder, stats[i][0].replace(' ', '_'))
         lib.gb2proteinortho(GBK, protortho, stats[i][0].replace(' ', '_'))
         eggnog.append(lib.getEggNogfromNote(GBK))
@@ -401,6 +403,28 @@ def __init__(self,prog):
     output.write(lib.FOOTER)
 ########################################################
 
+####SignalP############################
+#flip the dict and just count number for each
+signalpDict = lib.busco_dictFlip(signalp)
+
+if not os.path.isdir(os.path.join(args.out, 'signalp')):
+    os.makedirs(os.path.join(args.out, 'signalp'))
+sig = {}
+for i in range(0,len(scinames)):
+    if scinames[i] not in sig:
+        sig[scinames[i]] = len(signalpDict[i])     
+sigdf = pd.DataFrame([sig])
+sigdf.transpose().to_csv(os.path.join(args.out, 'signalp', 'signalp.csv'))
+lib.drawbarplot(sigdf, os.path.join(args.out, 'signalp', 'signalp.pdf'))
+
+#create html output
+with open(os.path.join(args.out, 'signalp.html'), 'w') as output:
+    output.write(lib.HEADER)
+    output.write(lib.SIGNALP)
+    output.write(lib.FOOTER)
+
+########################################################
+
 ####GO Terms, GO enrichment############################
 if not os.path.isdir(os.path.join(args.out, 'go_enrichment')):
     os.makedirs(os.path.join(args.out, 'go_enrichment'))
@@ -636,8 +660,9 @@ def __init__(self,prog):
 meropsDict = lib.dictFlip(merops)  
 cazyDict = lib.dictFlip(cazy)
 
+
 table = []
-header = ['GeneID','length','description', 'Ortho Group', 'EggNog', 'BUSCO','Protease family', 'CAZyme family', 'InterPro Domains', 'PFAM Domains', 'GO terms', 'SecMet Cluster', 'SMCOG']
+header = ['GeneID','length','description', 'Ortho Group', 'EggNog', 'BUSCO', 'Secreted', 'Protease family', 'CAZyme family', 'InterPro Domains', 'PFAM Domains', 'GO terms', 'SecMet Cluster', 'SMCOG']
 for y in range(0,num_input):
     outputname = os.path.join(args.out, 'annotations', scinames[y]+'.all.annotations.tsv')
     with open(outputname, 'w') as output:
@@ -657,6 +682,10 @@ def __init__(self,prog):
                             IPRdomains = "; ".join(iprDict.get(ID))
                         else:
                             IPRdomains = ''
+                        if ID in signalpDict[y]:
+                            signalphit = signalpDict[y].get(ID)[0]
+                        else:
+                            signalphit = ''
                         if ID in pfamDict:
                             pfamdomains = "; ".join(pfamDict.get(ID))
                         else:
@@ -693,7 +722,7 @@ def __init__(self,prog):
                                     if i.startswith('SMCOG:'):
                                         smcog = i
 
-                        final_result = [ID, str(length), description, orthogroup, egg, buscogroup, meropsdomains, cazydomains, IPRdomains, pfamdomains, goTerms, cluster, smcog]
+                        final_result = [ID, str(length), description, orthogroup, egg, buscogroup, signalphit, meropsdomains, cazydomains, IPRdomains, pfamdomains, goTerms, cluster, smcog]
                         output.write("%s\n" % ('\t'.join(final_result)))        
 ############################################
 

bin/funannotate-functional.py

@@ -252,7 +252,16 @@ def runIPRpython(Input):
 if not lib.checkannotations(busco_out):
     lib.runBUSCO(Proteins, 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') 
+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
+signalp_out = os.path.join(outputdir, 'annotate_misc', 'annotations.signalp.txt')
+if lib.which('signalp'):
+    if not lib.checkannotations(signalp_out):
+        lib.signalP(Proteins, os.path.join(outputdir, 'annotate_misc'), signalp_out)
+    num_annotations = lib.line_count(signalp_out)
+    lib.log.info('{0:,}'.format(num_annotations) + ' annotations added')
+else:
+    lib.log.info("SignalP not installed, skipping")
 
 if not args.skip_iprscan:
     if not args.iprscan:
@@ -345,8 +354,10 @@ def runIPRpython(Input):
     for file in os.listdir(os.path.join(outputdir, 'annotate_misc')):
         if file.startswith('annotations'):
             file = os.path.join(outputdir, 'annotate_misc', file)
-            with open(file, 'rU') as input:
+            with open(file) as input:
                 for line in input:
+                    if 'go.obo' in line:  #new goatools adds this damn line in my output, remove it here
+                        continue
                     total_annotations += 1
                     if line not in lines_seen:
                         output.write(line)

funannotate.py

@@ -31,7 +31,7 @@ def fmtcols(mylist, cols):
              for i in range(0,num_lines))
     return "\n".join(lines)
 
-version = '0.2.1'
+version = '0.2.2'
 
 default_help = """
 Usage:       funannotate <command> <arguments>

lib/library.py

@@ -91,11 +91,14 @@ def multipleReplace(text, wordDict):
 def which(name):
     try:
         with open(os.devnull) as devnull:
-            diff = ['tbl2asn', 'dustmasker', 'mafft']
+            diff = ['tbl2asn', 'dustmasker', 'mafft', 'signalp']
             if not any(name in x for x in diff):
                 subprocess.Popen([name], stdout=devnull, stderr=devnull).communicate()
             else:
-                subprocess.Popen([name, '--version'], stdout=devnull, stderr=devnull).communicate()
+                if name == 'signalp':
+                    subprocess.Popen([name, '-V'], stdout=devnull, stderr=devnull).communicate()
+                else:
+                    subprocess.Popen([name, '--version'], stdout=devnull, stderr=devnull).communicate()
     except OSError as e:
         if e.errno == os.errno.ENOENT:
             return False
@@ -573,6 +576,29 @@ def dbCANsearch(input, cpus, evalue, tmpdir, output):
                                 query = query + '-T1'
                             out.write("%s\tnote\tCAZy:%s\n" % (query, hit))
 
+def signalP(input, tmpdir, output):
+    log.info("Predicting secreted proteins with SignalP")
+    FNULL = open(os.devnull, 'w')
+    signalp_result = os.path.join(tmpdir, 'signalp.txt')
+    with open(signalp_result, 'w') as out:
+        subprocess.call(['signalp', '-t', 'euk', '-f', 'short', input], stdout = out, stderr=FNULL)
+    #parse output and turn into annotation file
+    with open(output, 'w') as signalp:
+        with open(signalp_result, 'rU') as results:
+            for line in results:
+                line = line.replace('\n', '')
+                if line.startswith('#'):
+                    continue
+                col = line.split(' ') #not tab delimited
+                col = filter(None, col) #clean up empty spaces
+                if col[9] == 'Y': #then there is signal peptide
+                    ID = col[0]
+                    if not ID.endswith('-T1'):
+                        ID = ID + '-T1'
+                    start = 1
+                    end = int(col[2]) - 1
+                    signalp.write("%s\tnote\tSECRETED:SignalP(%i-%i)\n" % (ID, start, end))
+                
 def RepeatModelMask(input, cpus, tmpdir, output, debug):
     log.info("Loading sequences and soft-masking genome")
     FNULL = open(os.devnull, 'w')
@@ -1263,6 +1289,21 @@ def drawHeatmap(df, color, output, annotate):
         item.set_fontsize(4)
     fig.savefig(output, format='pdf', dpi=1000, bbox_inches='tight')
     plt.close(fig)
+
+def drawbarplot(df, output):
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        import matplotlib.pyplot as plt
+        import seaborn as sns
+    #num = len(df.columns) + 1
+    sns.set(style="darkgrid")
+    fig = plt.figure()
+    ax = sns.barplot(data=df, palette=pref_colors)
+    plt.xlabel('Genomes')
+    plt.ylabel('Secreted Proteins')
+    fig.savefig(output, format='pdf', dpi=1000, bbox_inches='tight')
+    plt.close(fig) 
+
 
 def distance2mds(df, distance, type, output):
     import numpy as np
@@ -1515,18 +1556,20 @@ def ortho2phylogeny(folder, df, num, dict, cpus, bootstrap, tmpdir, outgroup, sp
             <span class="icon-bar"></span>
             <span class="icon-bar"></span>
             <span class="icon-bar"></span>
+            <span class="icon-bar"></span>
           </button>
           <a class="navbar-brand" href="index.html">Funannotate</a>
         </div>
         <div id="navbar" class="collapse navbar-collapse">
           <ul class="nav navbar-nav">
-            <li class="active"><a href="stats.html">Stats</a></li>
+            <li><a href="stats.html">Stats</a></li>
             <li><a href="phylogeny.html">Phylogeny</a></li>
             <li><a href="orthologs.html">Orthologs</a></li>
             <li><a href="interpro.html">InterProScan</a></li>
             <li><a href="pfam.html">PFAM</a></li>
             <li><a href="merops.html">Merops</a></li>
             <li><a href="cazy.html">CAZymes</a></li>
+            <li><a href="signalp.html">SignalP</a></li>
             <li><a href="go.html">GO ontology</a></li>
             <li><a href="citation.html">Citation</a></li>
           </ul>
@@ -1549,6 +1592,7 @@ def ortho2phylogeny(folder, df, num, dict, cpus, bootstrap, tmpdir, outgroup, sp
          <p><a href='phylogeny.html'>Maximum likelihood Phylogeny (RAxML)</a></p>
          <p><a href='merops.html'>MEROPS Protease Stats</a></p>
          <p><a href='cazy.html'>CAZyme carbohydrate activating enzyme Stats</a></p>
+         <p><a href='signal.html'>Secreted proteins (SignalP)</a></p>
          <p><a href='interpro.html'>InterProScan Domain Stats</a></p>
          <p><a href='pfam.html'>PFAM Domain Stats</a></p>
          <p><a href='go.html'>Gene Ontology Enrichment Analysis</a></p>
@@ -1593,6 +1637,14 @@ def ortho2phylogeny(folder, df, num, dict, cpus, bootstrap, tmpdir, outgroup, sp
         <a href='pfam/PFAM.nmds.pdf'><img src="pfam/PFAM.nmds.pdf" height="500" /></a></div>
         <div class="table-responsive">
 '''
+SIGNALP = '''
+    <div class="container">
+      <div class="starter-template">
+        <h2 class="sub-header">Secreted Proteins per Genome Results</h2>
+        <div class='row'>
+        <a href='signalp/signalp.pdf'><img src="signalp/signalp.pdf" height="500" /></a></div>
+        <div class="table-responsive">
+'''
 CAZY = '''
     <div class="container">
       <div class="starter-template">
@@ -1666,6 +1718,7 @@ def ortho2phylogeny(folder, df, num, dict, cpus, bootstrap, tmpdir, outgroup, sp
             <li><a href="pfam.html">PFAM</a></li>
             <li><a href="merops.html">Merops</a></li>
             <li><a href="cazy.html">CAZymes</a></li>
+            <li><a href="signalp.html">SignalP</a></li>
             <li><a href="go.html">GO ontology</a></li>
             <li><a href="citation.html">Citation</a></li>
             <li class="dropdown">