0
+ − 1 #!/usr/bin/env python
+ − 2 # -*- coding: UTF-8 -*-
+ − 3
+ − 4 import os, sys, subprocess, commands
+ − 5 import random, shutil
+ − 6 import zipfile
+ − 7 from Bio import SeqIO
+ − 8 import tempfile
+ − 9
+ − 10 blastdbpath = '/home/galaxy/bin/antismash-1.1.0/db'
+ − 11 pfamdbpath = '/home/galaxy/bin/antismash-1.1.0/db'
+ − 12 antismash_path = '/home/galaxy/bin/antismash-1.1.0/antismash.py'
+ − 13
+ − 14
+ − 15 def anitSMASH(args, sequence_path, glimmer_path):
+ − 16 #./antismash.py Tue6071_genome.fasta --geneclustertypes 1 --fullhmm y
+ − 17 print sequence_path
+ − 18 rint = random.randint(1,10000000)
+ − 19 tmp_dir = '/tmp/galaxy_%s' % rint
+ − 20 os.mkdir(tmp_dir)
+ − 21 os.mkdir(os.path.join( tmp_dir, 'geneprediction' ))
+ − 22 os.chdir(tmp_dir)
+ − 23 print 'IFORMAT:',args.iformat
+ − 24 if args.iformat == 'fasta':
+ − 25 new_input_path = os.path.join(tmp_dir, os.path.basename(sequence_path) + '.fasta')
+ − 26 elif args.iformat == 'genbank':
+ − 27 new_input_path = os.path.join(tmp_dir, os.path.basename(sequence_path) + '.gbk')
+ − 28 else:
+ − 29 new_input_path = os.path.join(tmp_dir, os.path.basename(sequence_path) + '.embl')
+ − 30
+ − 31
+ − 32 # try to generate the same name as in antismash.py
+ − 33 genomename = ".".join( (os.path.basename(sequence_path) + '.'+ args.iformat).split(".")[:-1] )
+ − 34 for i in """!"#$%&()*+,./:;=>?@[]^`{|}'""":
+ − 35 genomename = genomename.replace(i,"")
+ − 36 result_path = os.path.join( tmp_dir, genomename )
+ − 37
+ − 38 shutil.copy(sequence_path, new_input_path )
+ − 39
+ − 40 if args.eukaryotic:
+ − 41 taxon = '--taxon e'
+ − 42 else:
+ − 43 taxon = '--taxon p'
+ − 44
+ − 45 if args.clusterblast:
+ − 46 clusterblast = '--clusterblast y'
+ − 47 else:
+ − 48 clusterblast = '--clusterblast n'
+ − 49
+ − 50 if args.smcogs:
+ − 51 smcogs = '--smcogs y'
+ − 52 else:
+ − 53 smcogs = '--smcogs n'
+ − 54
+ − 55 if args.fullhmm:
+ − 56 fullhmm = '--fullhmm y'
+ − 57 else:
+ − 58 fullhmm = '--fullhmm n'
+ − 59
+ − 60 if args.fullblast:
+ − 61 fullblast = '--fullblast y'
+ − 62 else:
+ − 63 fullblast = '--fullblast n'
+ − 64
+ − 65 h = [antismash_path, new_input_path,
+ − 66 '--geneclustertypes %s' % args.geneclustertypes,
+ − 67 taxon,
+ − 68 clusterblast,
+ − 69 smcogs,
+ − 70 fullhmm,
+ − 71 fullblast,
+ − 72 '--blastdbpath %s' % blastdbpath,
+ − 73 '--pfamdbpath %s' % pfamdbpath,
+ − 74 '--cores 10',
+ − 75 ]
+ − 76 if args.iformat == 'fasta':
+ − 77 h.append('--glimmer_prediction %s' % glimmer_path)
+ − 78
+ − 79 a = ' '.join(h)
+ − 80
+ − 81 subprocess.call(a, shell=True)
+ − 82 print a
+ − 83 embl_path = os.path.join(result_path, '%s.final.embl' % genomename)
+ − 84 if os.path.exists(embl_path) and args.embl_path:
+ − 85 args.embl_path.write( open(embl_path).read() )
+ − 86
+ − 87 gff_line = "%s\tantismash\t%s\t%s\t%s\t.\t%s\t.\t%s"
+ − 88
+ − 89 clustername_mapping = {}
+ − 90 genecluster_path = os.path.join(result_path, 'clusterblast/geneclusters.txt')
+ − 91 if os.path.exists(genecluster_path):
+ − 92 for line in open( genecluster_path):
+ − 93 token = line.split('\t')
+ − 94 clustername_mapping[token[2]] = token[3]
+ − 95 old_cluster_name = False
+ − 96 parent_cluster_name = ""
+ − 97 child_lines = ""
+ − 98 starts = []
+ − 99 ends = []
+ − 100 for line in open( os.path.join(result_path, 'clusterblast/geneclusterprots.fasta') ):
+ − 101 if line.startswith('>'):
+ − 102 columns = line[1:].strip().split('|')
+ − 103 args.geneclusterprots.write( line.replace('|%s|' % columns[1], '|%s|%s|' % (columns[1],clustername_mapping[columns[1]])) )
+ − 104
+ − 105 """
+ − 106 "Handlampe" kann man zu einer Person sagen wenn es nicht mehr zum Armleuchter reicht. Es ist also ein abgeschwächte Form für "Depp".
+ − 107 """
+ − 108 if old_cluster_name == False:
+ − 109 # inital clustername setting
+ − 110 # try to catch the line in which the name changes
+ − 111 old_cluster_name = columns[1]
+ − 112 parent_cluster_name = columns[1]
+ − 113
+ − 114 seqname = columns[0]
+ − 115 strand = columns[3]
+ − 116 (start, end) = columns[2].split('-')
+ − 117 starts.append(int(start))
+ − 118 ends.append(int(end))
+ − 119
+ − 120 if old_cluster_name != columns[1]:
+ − 121 # begin of a new cluster, caclulate the parent line for the previous clusterproteins
+ − 122 # TODO: maybe the assumption is not correct, to take the smallest start and the largest end
+ − 123 starts.sort()
+ − 124 ends.sort()
+ − 125 genecluster_start = starts[0]
+ − 126 genecluster_end = ends[-1]
+ − 127 group = 'ID=%s;Name=%s;Clustertype=%s;Note=%s\n' % (parent_cluster_name, columns[4], clustername_mapping[columns[1]], columns[5])
+ − 128 oline = gff_line % (
+ − 129 seqname, 'genecluster', genecluster_start, genecluster_end, strand, group
+ − 130 )
+ − 131 args.geneclusterprots_gff.write(oline + child_lines)
+ − 132 starts = []
+ − 133 ends = []
+ − 134 old_cluster_name = columns[1]
+ − 135 parent_cluster_name = columns[1]
+ − 136 child_lines = ''
+ − 137
+ − 138 group = 'ID=%s;Parent=%s;Name=%s;Clustertype=%s;Note=%s\n' % (columns[4], columns[1], columns[4], clustername_mapping[columns[1]], columns[5])
+ − 139 child_lines += gff_line % (
+ − 140 seqname, 'gene', start, end, strand, group
+ − 141 )
+ − 142
+ − 143 else:
+ − 144 args.geneclusterprots.write( line )
+ − 145
+ − 146 starts.sort()
+ − 147 ends.sort()
+ − 148 genecluster_start = starts[0]
+ − 149 genecluster_end = ends[-1]
+ − 150 group = 'ID=%s;Name=%s;Clustertype=%s;Note=%s\n' % (parent_cluster_name, columns[4], clustername_mapping[columns[1]], columns[5])
+ − 151 oline = gff_line % (
+ − 152 seqname, 'genecluster', genecluster_start, genecluster_end, strand, group
+ − 153 )
+ − 154 args.geneclusterprots_gff.write(oline + child_lines)
+ − 155
+ − 156 # remove tmp directory
+ − 157 shutil.rmtree(tmp_dir)
+ − 158
+ − 159
+ − 160 def arg_parse():
+ − 161 import argparse
+ − 162 parser = argparse.ArgumentParser(prog = 'antiSMASH-Wrapper')
+ − 163 parser.add_argument('--version', action='version', version='%(prog)s 0.01')
+ − 164 parser.add_argument('--geneclustertypes')
+ − 165 parser.add_argument('--clusterblast', action='store_true')
+ − 166 parser.add_argument('--eukaryotic', action='store_true')
+ − 167 parser.add_argument('--fullhmm', action='store_true')
+ − 168 parser.add_argument('--smcogs', action='store_true')
+ − 169 parser.add_argument('--fullblast', action='store_true')
+ − 170
+ − 171 parser.add_argument('--input', '-i', help='FASTA Sequence File')
+ − 172 parser.add_argument('--glimmer_prediction', help='Glimmer Prediction File')
+ − 173
+ − 174 parser.add_argument('--input_format', '-f', dest="iformat", help='input format, if the input is a fasta file you need a glimmer file as additional input')
+ − 175
+ − 176 parser.add_argument('--zip', help='output: all files as zip file')
+ − 177 parser.add_argument('--html_file', help='output: the path to the index html file')
+ − 178 parser.add_argument('--html_path', help='output: the path to the output html dir')
+ − 179 parser.add_argument('--embl_path', help='output: the path to the embl output file', type=argparse.FileType('w'))
+ − 180 parser.add_argument('--geneclusterprots', help='output: Genecluster Fasta File', type=argparse.FileType('w'))
+ − 181 parser.add_argument('--geneclusterprots_gff', help='output: Genecluster GFF', type=argparse.FileType('w'))
+ − 182
+ − 183
+ − 184 args = parser.parse_args()
+ − 185 return args
+ − 186
+ − 187
+ − 188 def extract_glimmerHMM_results(sequence_id, glimmerHMM):
+ − 189 """
+ − 190 Returns the given annoations from a sequence_id as complete GFF3 formated string.
+ − 191 """
+ − 192 result = "##gff-version 3\n"
+ − 193 found_annoations = False
+ − 194 for line in open(glimmerHMM):
+ − 195 if line.startswith(sequence_id + '\t'):
+ − 196 result += line
+ − 197 found_annoations = True
+ − 198 if found_annoations:
+ − 199 return result
+ − 200 else:
+ − 201 return False
+ − 202
+ − 203 def extract_glimmer_results(sequence_id, glimmer):
+ − 204 """
+ − 205 Returns the given annoations from a sequence_id as glimmer result file.
+ − 206 """
+ − 207 result = ''
+ − 208 found_annoations = False
+ − 209 for line in open(glimmer):
+ − 210 # that construct matches the pseudo fasta header and the corresponding orfs
+ − 211 if line.split()[0].startswith(sequence_id):
+ − 212 result += line
+ − 213 found_annoations = True
+ − 214 if found_annoations:
+ − 215 return result
+ − 216 else:
+ − 217 return False
+ − 218
+ − 219
+ − 220 if __name__ == '__main__':
+ − 221 args = arg_parse()
+ − 222 geneclusterprots = ""
+ − 223 temp_dir = tempfile.mkdtemp()
+ − 224 args.geneclusterprots_gff.write("##gff-version 3\n")
+ − 225
+ − 226
+ − 227
+ − 228 for record in SeqIO.parse(open(args.input, "rU"), args.iformat) :
+ − 229 # create two tem files and fill it with the information from one sequence
+ − 230 record.id = record.id or record.name
+ − 231 tmp_sequence = os.path.join(temp_dir, record.id) #fasta -file
+ − 232 SeqIO.write(record, open(tmp_sequence, 'w+'), args.iformat)
+ − 233
+ − 234 tmp_glimmer = ''
+ − 235 if args.iformat == 'fasta':
+ − 236 # We need to annotate the sequences manually
+ − 237 if args.eukaryotic:
+ − 238 #print os.path.join(temp_dir, record.id + '.gff')
+ − 239 tmp_glimmer = os.path.join(temp_dir, record.id + '.gff')
+ − 240 annotations = extract_glimmerHMM_results(record.id, args.glimmer_prediction)
+ − 241 #print annotations, ':', record.id, args.glimmer_prediction
+ − 242 if not annotations:
+ − 243 continue
+ − 244 open(tmp_glimmer, 'w+').write( annotations )
+ − 245 else:
+ − 246 tmp_glimmer = os.path.join(temp_dir, record.id + '.tsv')
+ − 247 annotations = extract_glimmer_results(record.id, args.glimmer_prediction)
+ − 248 if not annotations:
+ − 249 continue
+ − 250 open(tmp_glimmer, 'w+').write( annotations )
+ − 251
+ − 252 anitSMASH(args, tmp_sequence, tmp_glimmer)
+ − 253 shutil.rmtree(temp_dir)
+ − 254
+ − 255
+ − 256
+ − 257
+ − 258
+ − 259
+ − 260
+ − 261
+ − 262
+ − 263
+ − 264
+ − 265
