Mercurial > repos > vimalkumarvelayudhan > viga
comparison VIGA.py @ 0:231e4c669675 draft
Initial commit - v0.10.3 git commit deeded0
| author | vimalkumarvelayudhan |
|---|---|
| date | Tue, 27 Feb 2018 14:16:54 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:231e4c669675 |
|---|---|
| 1 #!/usr/bin/env python | |
| 2 | |
| 3 # -*- coding: utf-8 -*- | |
| 4 | |
| 5 # VIGA - De-novo VIral Genome Annotator | |
| 6 # | |
| 7 # Copyright (C) 2017 - Enrique Gonzalez-Tortuero | |
| 8 # Vimalkumar Velayudhan | |
| 9 # | |
| 10 # This program is free software: you can redistribute it and/or modify | |
| 11 # it under the terms of the GNU General Public License as published by | |
| 12 # the Free Software Foundation, either version 3 of the License, or | |
| 13 # (at your option) any later version. | |
| 14 # | |
| 15 # This program is distributed in the hope that it will be useful, | |
| 16 # but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 18 # GNU General Public License for more details. | |
| 19 # | |
| 20 # You should have received a copy of the GNU General Public License | |
| 21 # along with this program. If not, see <http://www.gnu.org/licenses/>. | |
| 22 | |
| 23 # Importing python libraries | |
| 24 from __future__ import print_function | |
| 25 import argparse | |
| 26 import csv | |
| 27 import fileinput | |
| 28 import fractions | |
| 29 import glob | |
| 30 import numpy | |
| 31 import os | |
| 32 import os.path | |
| 33 import re | |
| 34 import sys | |
| 35 import shutil | |
| 36 import subprocess | |
| 37 from BCBio import GFF | |
| 38 from Bio import SeqIO | |
| 39 from Bio import SeqFeature | |
| 40 from Bio.Alphabet import IUPAC | |
| 41 from Bio.Seq import Seq | |
| 42 from Bio.SeqFeature import FeatureLocation | |
| 43 from Bio.SeqRecord import SeqRecord | |
| 44 from Bio.SeqUtils.ProtParam import ProteinAnalysis | |
| 45 from collections import OrderedDict, defaultdict | |
| 46 from itertools import product | |
| 47 from scipy import signal | |
| 48 from time import strftime | |
| 49 | |
| 50 # Preparing functions | |
| 51 def batch_iterator(iterator, batch_size): | |
| 52 entry = True | |
| 53 while entry: | |
| 54 batch = [] | |
| 55 while len(batch) < batch_size: | |
| 56 try: | |
| 57 entry = iterator.next() | |
| 58 except StopIteration: | |
| 59 entry = None | |
| 60 if entry is None: | |
| 61 break | |
| 62 batch.append(entry) | |
| 63 if batch: | |
| 64 yield batch | |
| 65 | |
| 66 def check_peaks(peaks, length): | |
| 67 # Checking if origin/terminus peaks are too close or too far apart. In that case, they are probably wrong | |
| 68 closest, farthest = int(length * float(0.45)), int(length * float(0.55)) | |
| 69 pairs = [] | |
| 70 for pair in list(product(*peaks)): | |
| 71 ### added this to make sure gets origin and ter right | |
| 72 tr, pk = sorted(list(pair), key = lambda x: x[1], reverse = False) # trough and peak | |
| 73 a = (tr[0] - pk[0]) % length | |
| 74 b = (pk[0] - tr[0]) % length | |
| 75 pt = abs(tr[1] - pk[1]) # distance between values | |
| 76 if (a <= farthest and a >= closest) or (b <=farthest and b >= closest): | |
| 77 pairs.append([pt, tr, pk]) | |
| 78 if len(pairs) == 0: | |
| 79 return [False, False] | |
| 80 pt, tr, pk = sorted(pairs, reverse = True)[0] | |
| 81 return [tr[0], pk[0]] | |
| 82 | |
| 83 def cmd_exists(cmd): | |
| 84 return subprocess.call("type " + cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0 | |
| 85 | |
| 86 def GCskew(name, length, seq, window, slide): | |
| 87 replacements = {'G':1, 'C':-1, 'A':0, 'T':0, 'N':0} | |
| 88 gmc = [] # G - C | |
| 89 for base in seq: | |
| 90 try: | |
| 91 gmc.append(replacements[base]) | |
| 92 except: | |
| 93 gmc.append(0) | |
| 94 # convert to G + C | |
| 95 gpc = [abs(i) for i in gmc] # G + C | |
| 96 # calculate sliding windows for (G - C) and (G + C) | |
| 97 weights = numpy.ones(window)/window | |
| 98 gmc = [[i, c] for i, c in enumerate(signal.fftconvolve(gmc, weights, 'same').tolist())] | |
| 99 gpc = [[i, c] for i, c in enumerate(signal.fftconvolve(gpc, weights, 'same').tolist())] | |
| 100 # calculate gc skew and cummulative gc skew sum | |
| 101 skew = [[], []] # x and y for gc skew | |
| 102 c_skew = [[], []] # x and y for gc skew cummulative sums | |
| 103 cs = 0 # cummulative sum | |
| 104 # select windows to use based on slide | |
| 105 for i, m in gmc[0::slide]: | |
| 106 p = gpc[i][1] | |
| 107 if p == 0: | |
| 108 gcs = 0 | |
| 109 else: | |
| 110 gcs = m/p | |
| 111 cs += gcs | |
| 112 skew[0].append(i) | |
| 113 c_skew[0].append(i) | |
| 114 skew[1].append(gcs) | |
| 115 c_skew[1].append(cs) | |
| 116 ori, ter = find_ori_ter(c_skew, length) | |
| 117 return ori, ter, skew, c_skew | |
| 118 | |
| 119 def eprint(*args, **kwargs): | |
| 120 print(*args, file=sys.stderr, **kwargs) | |
| 121 | |
| 122 def find_ori_ter(c_skew, length): | |
| 123 # Find origin and terminus of replication based on cumulative GC skew min and max peaks | |
| 124 c_skew_min = signal.argrelextrema(numpy.asarray(c_skew[1]), numpy.less, order = 1)[0].tolist() | |
| 125 c_skew_max = signal.argrelextrema(numpy.asarray(c_skew[1]), numpy.greater, order = 1)[0].tolist() | |
| 126 # return False if no peaks were detected | |
| 127 if len(c_skew_min) == 0 or len(c_skew_min) == 0: | |
| 128 return [False, False] | |
| 129 else: | |
| 130 c_skew_min = [[c_skew[0][i], c_skew[1][i]] for i in c_skew_min] | |
| 131 c_skew_max = [[c_skew[0][i], c_skew[1][i]] for i in c_skew_max] | |
| 132 ori, ter = check_peaks([c_skew_min, c_skew_max], length) | |
| 133 return ori, ter | |
| 134 | |
| 135 def stringSplitByNumbers(x): | |
| 136 r = re.compile('(\d+)') | |
| 137 l = r.split(x) | |
| 138 return [int(y) if y.isdigit() else y for y in l] | |
| 139 | |
| 140 # Defining the program version | |
| 141 version = "0.10.3" | |
| 142 | |
| 143 # Processing the parameters | |
| 144 parser = argparse.ArgumentParser(description='VIGA is an automatic de novo VIral Genome Annotator.') | |
| 145 basic_group = parser.add_argument_group('Basic options for VIGA [REQUIRED]') | |
| 146 | |
| 147 basic_group.add_argument("--input", dest="inputfile", type=str, required=True, help='Input file as a FASTA file', metavar="FASTAFILE") | |
| 148 basic_group.add_argument("--rfamdb", dest="rfamdatabase", type=str, required=True, help='RFAM Database that will be used for the ribosomal RNA prediction. RFAMDB should be in the format "/full/path/to/rfamdb/Rfam.cm" and must be compressed accordingly (see INFERNAL manual) before running the script.', metavar="RFAMDB") | |
| 149 basic_group.add_argument("--modifiers", dest="modifiers", type=str, required=True, help='Input file as a plain text file with the modifiers per every FASTA header according to SeqIn (https://www.ncbi.nlm.nih.gov/Sequin/modifiers.html). All modifiers must be written in a single line and are separated by a single space character. No space should be placed besides the = sign. For example: [organism=Serratia marcescens subsp. marcescens] [sub-species=marcescens] [strain=AH0650_Sm1] [topology=linear] [moltype=DNA] [tech=wgs] [gcode=11] [country=Australia] [isolation-source=sputum]. This line will be copied and printed along with the record name as the definition line of every contig sequence.', metavar="TEXTFILE") | |
| 150 | |
| 151 advanced_group = parser.add_argument_group('Advanced options for VIGA [OPTIONAL]') | |
| 152 advanced_group.add_argument("--readlength", dest="read_length", type=int, default=101, help='Read length for the circularity prediction (default: 101 bp)', metavar="INT") | |
| 153 advanced_group.add_argument("--windowsize", dest="window", type=int, default=100, help='Window size used to determine the origin of replication in circular contigs according to the cumulative GC skew (default: 100 bp)', metavar="INT") | |
| 154 advanced_group.add_argument("--slidingsize", dest="slide", type=int, default=10, help='Window size used to determine the origin of replication in circular contigs according to the cumulative GC skew (default: 10 bp)', metavar="INT") | |
| 155 advanced_group.add_argument("--out", dest="rootoutput", type=str, help='Name of the outputs files (without extension)', metavar="OUTPUTNAME") | |
| 156 advanced_group.add_argument("--locus", dest="locus", type=str, default='LOC', help='Name of the sequences. If the input is a multiFASTA file, please put a general name as the program will add the number of the contig at the end of the name (Default: %(default)s)', metavar="STRING") | |
| 157 advanced_group.add_argument("--gff", dest="gffprint", action='store_true', default=False, help='Printing the output as GFF3 file (Default: False)') | |
| 158 advanced_group.add_argument("--threads", dest="ncpus", default=1, help='Number of threads/cpus (Default: %(default)s cpu)', metavar="INT") | |
| 159 advanced_group.add_argument("--nohmmer", dest="nohmmer", action='store_true', default=False, help='Running only BLAST to predict protein function. (Default: False)') | |
| 160 advanced_group.add_argument("--noblast", dest="noblast", action='store_true', default=False, help='Running DIAMOND instead of BLAST to predict protein function according to homology. This will be less sensitive but faster than BLAST. (Default: False)') | |
| 161 advanced_group.add_argument("--blastdb", dest="blastdatabase", type=str, help='BLAST Database that will be used for the protein function prediction. The database must be an amino acid one, not nucleotidic. This argument is mandatory if --noblast option is disabled', metavar="BLASTDB") | |
| 162 advanced_group.add_argument("--diamonddb", dest="diamonddatabase", type=str, help='DIAMOND Database that will be used for the protein function prediction. The database must be created from a amino acid FASTA file as indicated in https://github.com/bbuchfink/diamond. This argument is mandatory when --noblast option is enabled', metavar="DIAMONDDB") | |
| 163 advanced_group.add_argument("--blastevalue", dest="blastevalue", default=0.00001, help='BLAST/DIAMOND e-value threshold (Default: 0.00001)', metavar="FLOAT") | |
| 164 advanced_group.add_argument("--hmmdb", dest="hmmdatabase", type=str, help='PHMMER Database that will be used for the protein function prediction according to Hidden Markov Models. In this case, HMMDB must be in FASTA format (e.g. UniProt: "', metavar="HMMDB") | |
| 165 advanced_group.add_argument("--hmmerevalue", dest="hmmerevalue", default=0.001, help='PHMMER e-value threshold (Default: 0.001)', metavar="FLOAT") | |
| 166 | |
| 167 type_choices = {'BCT': 'Prokaryotic chromosome', 'CON': 'Contig', 'PHG': 'Phages', 'VRL': 'Eukaryotic/Archaea virus'} | |
| 168 type_help = ('GenBank Division: One of the following codes - {0}. (Default: %(default)s)'.format(', '.join('{0} ({1})'.format(k, v) for k, v in type_choices.items()))) | |
| 169 advanced_group.add_argument("--typedata", dest="typedata", type=str, default='CON', help=type_help, metavar="BCT|CON|VRL|PHG") | |
| 170 | |
| 171 gcode_choices = {'1': 'Standard genetic code [Eukaryotic]', '2': 'Vertebrate mitochondrial code', '3': 'Yeast mitochondrial code', '4': 'Mycoplasma/Spiroplasma and Protozoan/mold/coelenterate mitochondrial code', '5': 'Invertebrate mitochondrial code', '6': 'Ciliate, dasycladacean and hexamita nuclear code', '9': 'Echinoderm/flatworm mitochondrial code', '10': 'Euplotid nuclear code', '11': 'Bacteria/Archaea/Phages/Plant plastid', '12': 'Alternative yeast nuclear code', '13': 'Ascidian mitochondrial code', '14': 'Alternative flatworm mitochondrial code', '16': 'Chlorophycean mitochondrial code', '21': 'Trematode mitochondrial code', '22': 'Scedenesmus obliquus mitochondrial code', '23': 'Thraustochytrium mitochondrial code', '24': 'Pterobranquia mitochondrial code', '25': 'Gracilibacteria & Candidate division SR1', '26': 'Pachysolen tannophilus nuclear code', '27': 'Karyorelict nuclear code', '28': 'Condylostoma nuclear code', '29': 'Mesodinium nuclear code', '30': 'Peritrich nuclear code', '31': 'Blastocrithidia nuclear code'} | |
| 172 gcode_help = ('Number of GenBank translation table. At this moment, the available options are {0}. (Default: %(default)s)'.format('{}'.format(', '.join('{0} ({1})'.format(k, v) for k, v in gcode_choices.items())))) | |
| 173 advanced_group.add_argument("--gcode", dest="gcode", type=str, default='11', help=gcode_help, metavar="NUMBER") | |
| 174 | |
| 175 advanced_group.add_argument('--mincontigsize', dest="mincontigsize", type=int, default = 200, help = 'Minimum contig length to be considered in the final files (Default: 200 bp)', metavar="INT") | |
| 176 advanced_group.add_argument("--idthr", dest="idthreshold", default=50.00, help='ID threshold (Default: 50.0)', metavar="FLOAT") | |
| 177 advanced_group.add_argument("--coverthr", dest="covthreshold", default=50.00, help='Coverage threshold (Default: 50.0)', metavar="FLOAT") | |
| 178 advanced_group.add_argument("--diffid", dest="diffid", default=5.00, help='Max allowed difference between the ID percentages of BLAST and HMMER. (Default = 5.00; Not recommended to change such value)', metavar="FLOAT (>0.01)") | |
| 179 advanced_group.add_argument("--minrepeat", dest="minrepeat", type=int, default=16, help="Minimum repeat length for CRISPR detection (Default: 16)", metavar="INT") | |
| 180 advanced_group.add_argument("--maxrepeat", dest="maxrepeat", type=int, default=64, help="Maximum repeat length for CRISPR detection (Default: 64)") | |
| 181 advanced_group.add_argument("--minspacer", dest="minspacer", type=int, default=8, help="Minimum spacer length for CRISPR detection (Default: 8)") | |
| 182 advanced_group.add_argument("--maxspacer", dest="maxspacer", type=int, default=64, help="Maximum spacer length for CRISPR detection (Default: 64)") | |
| 183 advanced_group.add_argument("--blastexh", dest="blastexh", action='store_true', default=False, help='Use of exhaustive BLAST to predict the proteins by homology according to Fozo et al. (2010) Nucleic Acids Res (Default=FALSE)') | |
| 184 | |
| 185 args = parser.parse_args() | |
| 186 | |
| 187 root_output = args.rootoutput | |
| 188 if not root_output: | |
| 189 root_output = '{}_annotated'.format(os.path.splitext(args.inputfile)[0]) | |
| 190 | |
| 191 if args.noblast == False and args.blastdatabase == None: | |
| 192 sys.exit('You MUST specify BLAST database using the parameter --blastdb if you are not using --noblast option') | |
| 193 | |
| 194 if args.noblast == True and args.diamonddatabase == None: | |
| 195 sys.exit('You MUST specify DIAMOND database using the parameter --diamonddb if you are using --noblast option') | |
| 196 | |
| 197 if args.nohmmer == False and args.hmmdatabase == None: | |
| 198 sys.exit('You MUST specify HMMER database using the parameter --hmmdb if you are not using --nohmmer option') | |
| 199 | |
| 200 # Printing the header of the program | |
| 201 eprint("This is VIGA %s" % str(version)) | |
| 202 eprint("Written by Enrique Gonzalez Tortuero & Vimalkumar Velayudhan") | |
| 203 eprint("Homepage is https://github.com/EGTortuero/virannot") | |
| 204 eprint("Local time: ", strftime("%a, %d %b %Y %H:%M:%S")) | |
| 205 eprint("\n\n") | |
| 206 | |
| 207 # checking the presence of the programs in the system | |
| 208 | |
| 209 if not cmd_exists("lastz")==True: | |
| 210 sys.exit("You must install LASTZ before using this script") | |
| 211 elif not cmd_exists("cmscan")==True: | |
| 212 sys.exit("You must install INFERNAL before using this script") | |
| 213 elif not cmd_exists("prodigal")==True: | |
| 214 sys.exit("You must install prodigal before using this script") | |
| 215 elif not cmd_exists("parallel")==True: | |
| 216 sys.exit("You must install GNU Parallel before using this script") | |
| 217 elif not cmd_exists("blastp")==True: | |
| 218 sys.exit("You must install BLAST before using this script") | |
| 219 elif not cmd_exists("diamond")==True: | |
| 220 sys.exit("You must install DIAMOND before using this script") | |
| 221 elif not cmd_exists("phmmer")==True: | |
| 222 sys.exit("You must install HMMER 3 before using this script") | |
| 223 elif not cmd_exists("aragorn")==True: | |
| 224 sys.exit("You must install ARAGORN before using this script") | |
| 225 elif not cmd_exists("pilercr")==True: | |
| 226 sys.exit("You must install PILERCR before using this script") | |
| 227 elif not cmd_exists("trf")==True: | |
| 228 sys.exit("You must install Tandem Repeats Finder before using this script") | |
| 229 elif not cmd_exists("irf")==True: | |
| 230 sys.exit("You must install Inverted Repeats Finder before using this script") | |
| 231 | |
| 232 eprint("Data type is {0} and GenBank translation table no is {1}\n".format(args.typedata, args.gcode)) | |
| 233 | |
| 234 # Correcting the original file (long headers problem + multiple FASTA files) | |
| 235 record_iter = SeqIO.parse(open(args.inputfile, "rU"),"fasta") | |
| 236 counter = 1 | |
| 237 newnamessequences = {} | |
| 238 for i, batch in enumerate(batch_iterator(record_iter, 1)): | |
| 239 seq_index = "CONTIG_%i" % (i + 1) | |
| 240 filename = "%s.temp.fna" % seq_index | |
| 241 newfilename = "%s.fna" % seq_index | |
| 242 with open(filename, "w") as handle: | |
| 243 count = SeqIO.write(batch, filename, "fasta") | |
| 244 | |
| 245 with open(filename, "rU") as original, open(newfilename, "w") as corrected: | |
| 246 sequences = SeqIO.parse(original, "fasta", IUPAC.ambiguous_dna) | |
| 247 for record in sequences: | |
| 248 original_name = record.id | |
| 249 record.id = "%s_%i" % (args.locus, counter) | |
| 250 record.description = record.description | |
| 251 counter += 1 | |
| 252 newnamessequences[record.id] = original_name | |
| 253 eprint("WARNING: The name of the sequence %s was corrected as %s" % (original_name, record.id)) | |
| 254 SeqIO.write(record, corrected, "fasta") | |
| 255 | |
| 256 with open("logfile.txt", "w") as logfile: | |
| 257 logfile.write("#Original\tNew\n") | |
| 258 for oldname in sorted(newnamessequences, key = stringSplitByNumbers): | |
| 259 logfile.write("%s\t%s\n" % (oldname, newnamessequences[oldname])) | |
| 260 os.remove(filename) | |
| 261 | |
| 262 for newfile in sorted(glob.glob("CONTIG_*.fna")): | |
| 263 | |
| 264 # Predicting the shape of the contig (code based on Alex Crits-Christoph's find_circular.py script [https://github.com/alexcritschristoph/VICA/blob/master/find_circular.py]) | |
| 265 eprint("Predicting the shape of the contig using LASTZ") | |
| 266 genomeshape = {} | |
| 267 with open(newfile, "rU") as targetfasta: | |
| 268 Sequence = SeqIO.parse(newfile, "fasta") | |
| 269 for record in Sequence: | |
| 270 seq_beginning = str(record.seq[0:args.read_length]) | |
| 271 seq_ending = str(record.seq[len(record.seq)-args.read_length:len(record.seq)]) | |
| 272 combined_seqs = SeqRecord(Seq(seq_beginning + seq_ending, IUPAC.ambiguous_dna), id = record.description) | |
| 273 SeqIO.write(combined_seqs, "temporal_circular.fasta", "fasta") | |
| 274 outputlastz = subprocess.check_output(["lastz", "temporal_circular.fasta", "--self", "--notrivial", "--nomirror", "--format=general-:start1,end1,start2,end2,score,strand1,strand2,identity,length1"]) | |
| 275 resultslastz = outputlastz.split("\n") | |
| 276 for resultlastz in resultslastz: | |
| 277 if resultlastz != '': | |
| 278 start1 = resultlastz.split()[0] | |
| 279 end1 = resultlastz.split()[1] | |
| 280 start2 = resultlastz.split()[2] | |
| 281 end2 = resultlastz.split()[3] | |
| 282 strand1 = resultlastz.split()[5] | |
| 283 strand2 = resultlastz.split()[6] | |
| 284 identity = resultlastz.split()[7] | |
| 285 length = int(resultlastz.split()[9]) | |
| 286 if strand1 == strand2 and length > 0.4 * args.read_length and float(fractions.Fraction(identity)) > 0.95 and int(start1) < 5 and int(start2) > args.read_length and int(end1) < args.read_length and int(end2) > args.read_length * 2 * 0.9: | |
| 287 genomeshape['genomeshape'] = "circular" | |
| 288 try: | |
| 289 if genomeshape['identity'] >= float(fractions.Fraction(identity)): | |
| 290 genomeshape['identity'] = float(fractions.Fraction(identity)) | |
| 291 genomeshape['length'] = length | |
| 292 except KeyError: | |
| 293 genomeshape['identity'] = float(fractions.Fraction(identity)) | |
| 294 genomeshape['length'] = length | |
| 295 else: | |
| 296 continue | |
| 297 if strand1 == strand2 and length > 0.4 * args.read_length and float(fractions.Fraction(identity)) > 0.95 and int(start1) < 5 and int(start2) > args.read_length and int(end1) < args.read_length and int(end2) > args.read_length * 2 * 0.9: | |
| 298 genomeshape['genomeshape'] = "circular" | |
| 299 try: | |
| 300 if genomeshape['identity'] >= float(fractions.Fraction(identity)): | |
| 301 genomeshape['identity'] = float(fractions.Fraction(identity)) | |
| 302 genomeshape['length'] = length | |
| 303 except KeyError: | |
| 304 genomeshape['identity'] = float(fractions.Fraction(identity)) | |
| 305 genomeshape['length'] = length | |
| 306 try: | |
| 307 if genomeshape['genomeshape'] == "": | |
| 308 genomeshape['genomeshape'] = "linear" | |
| 309 except KeyError: | |
| 310 genomeshape['genomeshape'] = "linear" | |
| 311 else: | |
| 312 genomeshape['genomeshape'] = "circular" | |
| 313 with open("temp.fasta", "w") as correctedcircular: | |
| 314 Corrseq = str(record.seq[int(genomeshape['length'])//2:-int(genomeshape['length'])//2]) | |
| 315 Newseq = SeqRecord(Seq(Corrseq, IUPAC.ambiguous_dna), id = record.description) | |
| 316 SeqIO.write(Newseq, correctedcircular, "fasta") | |
| 317 os.rename("temp.fasta", "temp2.fasta") | |
| 318 eprint("LASTZ predicted that %s is %s\n" % (newfile, genomeshape['genomeshape'])) | |
| 319 os.remove("temporal_circular.fasta") | |
| 320 | |
| 321 # Calculate the cumulative GCskew in circular contigs to determine the origin of replication (Based on iRep -- Brown CT, Olm MR, Thomas BC, Banfield JF (2016) Measurement of bacterial replication rates in microbial communities. Nature Biotechnology 34: 1256-63.) | |
| 322 if genomeshape['genomeshape'] == "circular": | |
| 323 for record in SeqIO.parse("temp2.fasta", "fasta"): | |
| 324 length_contig = len(record.seq) | |
| 325 #if length < min_len: | |
| 326 # print('%s: Too Short' % (name), file=sys.stderr) | |
| 327 # continue | |
| 328 oric, term, gcskew, cgcskew = GCskew(record.id, length_contig, record.seq, args.window, args.slide) | |
| 329 valoric = oric | |
| 330 if oric == False: | |
| 331 oric, term = 'n/a', 'n/a' | |
| 332 else: | |
| 333 oric, term = '{:,}'.format(oric), '{:,}'.format(term) | |
| 334 eprint('%s -> Origin: %s Terminus: %s' % (record.id, oric, term)) | |
| 335 #print('\t'.join(['# Name', 'Position', 'GC Skew', 'Cumulative GC Skew'])) | |
| 336 #for i, pos in enumerate(gcskew[0]): | |
| 337 # out = [record.id, pos, gcskew[1][i], cgcskew[1][i]] | |
| 338 # print('\t'.join([str(i) for i in out])) | |
| 339 if valoric != False: | |
| 340 firstpartseq = str(record.seq[valoric:-1]) | |
| 341 secondpartseq = str(record.seq[0:(valoric-1)]) | |
| 342 combinedcorrectedseq = SeqRecord(Seq(firstpartseq + secondpartseq, IUPAC.ambiguous_dna), id = record.description) | |
| 343 SeqIO.write(combinedcorrectedseq, newfile, "fasta") | |
| 344 else: | |
| 345 eprint("As the program was unable to predict the origin of replication, %s was considered as is without correcting!" % record.id) | |
| 346 os.rename("temp2.fasta", newfile) | |
| 347 if os.path.isfile("temp2.fasta"): | |
| 348 os.remove("temp2.fasta") | |
| 349 | |
| 350 # Predicting genes using PRODIGAL | |
| 351 eprint("\nRunning Prodigal to predict the genes in %s" % newfile) | |
| 352 for record in SeqIO.parse(newfile, "fasta"): | |
| 353 length_contig = len(record.seq) | |
| 354 if (length_contig >= 100000): | |
| 355 if genomeshape['genomeshape'] == 'linear': | |
| 356 subprocess.call(["prodigal", "-a", "pretemp.faa", "-i", newfile, "-o", "/dev/null", "-g", args.gcode, "-c", "-q"]) | |
| 357 else: | |
| 358 subprocess.call(["prodigal", "-a", "pretemp.faa", "-i", newfile, "-o", "/dev/null", "-g", args.gcode, "-q"]) | |
| 359 else: | |
| 360 if genomeshape['genomeshape'] == 'linear': | |
| 361 subprocess.call(["prodigal", "-a", "pretemp.faa", "-i", newfile, "-o", "/dev/null", "-p", "meta", "-g", args.gcode, "-c", "-q"]) | |
| 362 else: | |
| 363 subprocess.call(["prodigal", "-a", "pretemp.faa", "-i", newfile, "-o", "/dev/null", "-p", "meta", "-g", args.gcode, "-q"]) | |
| 364 num_seqs = len(list(SeqIO.parse("pretemp.faa", "fasta"))) | |
| 365 eprint("PRODIGAL was able to predict %i genes in %s\n" % (num_seqs, newfile)) | |
| 366 | |
| 367 with open("pretemp.faa", "rU") as originalfaa, open("temp.faa", "w") as correctedfaa: | |
| 368 sequences = SeqIO.parse(originalfaa, "fasta") | |
| 369 for record in sequences: | |
| 370 record.seq = record.seq.rstrip("*") | |
| 371 SeqIO.write(record, correctedfaa, "fasta") | |
| 372 | |
| 373 faa_file = "%s.faa" % newfile # TO DEBUG | |
| 374 shutil.copyfile("temp.faa", faa_file) # TO DEBUG | |
| 375 os.remove("pretemp.faa") | |
| 376 | |
| 377 # Predicting the function of the proteins based on homology using BLAST | |
| 378 equivalence = {} | |
| 379 record_iter = SeqIO.parse(open("temp.faa"),"fasta") | |
| 380 for i, batch in enumerate(batch_iterator(record_iter, 1)): | |
| 381 seq_index = "SEQ_%i" % (i + 1) | |
| 382 filename = "%s.faa" % seq_index | |
| 383 with open(filename, "w") as handle: | |
| 384 count = SeqIO.write(batch, handle, "fasta") | |
| 385 equivalence[seq_index] = batch[0].id | |
| 386 | |
| 387 if args.blastexh==True: | |
| 388 eprint("Running BLAST to predict the genes according to homology inference in %s using exhaustive mode (see Fozo et al. (2010) Nucleic Acids Res for details)" % newfile) | |
| 389 subprocess.call(['blastp', '-query', "temp.faa", '-db', args.blastdatabase, '-evalue', str(args.blastevalue), '-outfmt', '6 qseqid sseqid pident length qlen slen qstart qend evalue bitscore stitle', '-out', 'temp_blast.csv', '-max_target_seqs', '10', '-word_size', '2', '-gapopen', '8', '-gapextend', '2', '-matrix', '"PAM70"', '-comp_based_stats', '"0"', "-num_threads", str(args.ncpus)]) | |
| 390 eprint("Done. BLAST was done to predict the genes by homology\n") | |
| 391 blast_log = "%s.blast.log" % newfile # TO DEBUG | |
| 392 shutil.copyfile("temp_blast.csv", blast_log) # TO DEBUG | |
| 393 elif args.noblast==True: | |
| 394 eprint("Running DIAMOND to predict the genes according to homology inference in %s using default parameters" % newfile) | |
| 395 with open("temp.faa", "r") as tempfile: | |
| 396 first_line = tempfile.readline() | |
| 397 if first_line.startswith(">"): | |
| 398 subprocess.call(['diamond', 'blastp', '-q', 'temp.faa', '-d', args.diamonddatabase, '-e', str(args.blastevalue), '-f', '6', 'qseqid', 'sseqid', 'pident', 'length', 'qlen', 'slen', 'qstart', 'qend', 'evalue', 'bitscore', 'stitle', '-o', 'temp_blast.csv', '-k', '10', "-p", str(args.ncpus), '--quiet']) | |
| 399 else: | |
| 400 open("temp_blast.csv", 'a').close() | |
| 401 blast_log = "%s.blast.log" % newfile # TO DEBUG | |
| 402 shutil.copyfile("temp_blast.csv", blast_log) # TO DEBUG | |
| 403 eprint("Done. DIAMOND was done to predict the genes by homology\n") | |
| 404 else: | |
| 405 eprint("Running BLAST to predict the genes according to homology inference in %s using default parameters" % newfile) | |
| 406 subprocess.call(['blastp', '-query', "temp.faa", '-db', args.blastdatabase, '-evalue', str(args.blastevalue), '-outfmt', '6 qseqid sseqid pident length qlen slen qstart qend evalue bitscore stitle', '-out', 'temp_blast.csv', '-max_target_seqs', '10', "-num_threads", str(args.ncpus)]) | |
| 407 blast_log = "%s.blast.log" % newfile | |
| 408 shutil.copyfile("temp_blast.csv", blast_log) # TO DEBUG | |
| 409 eprint("Done. BLAST was done to predict the genes by homology\n") # TO DEBUG | |
| 410 | |
| 411 # Parsing the results from BLAST | |
| 412 with open("temp_blast.csv", "rU") as blastresults: | |
| 413 hypotheticalpat = re.compile(r'(((((?i)hypothetical)|(?i)uncharacteri(z|s)ed|(?i)predicted)) protein)|((?i)ORF|((?i)unnamed protein product|(?i)gp\d+))') | |
| 414 reader = csv.DictReader(blastresults, delimiter='\t', fieldnames=['qseqid','sseqid','pident','length','qlen','slen','qstart','qend','evalue','bitscore','stitle']) | |
| 415 information_proteins_blast = {} | |
| 416 for row in reader: | |
| 417 perc_cover = round(100.00*(float(row['length'])/float(row['qlen'])),2) | |
| 418 perc_id = float(row['pident']) | |
| 419 infoprot_blast = {} | |
| 420 infoprot_blast['sseqid'] = row['sseqid'] | |
| 421 infoprot_blast['pident'] = perc_id | |
| 422 infoprot_blast['pcover'] = perc_cover | |
| 423 infoprot_blast['evalue'] = row['evalue'] | |
| 424 infoprot_blast['descr'] = row['stitle'] | |
| 425 try: | |
| 426 data = information_proteins_blast[row['qseqid']] | |
| 427 except KeyError: | |
| 428 if not re.search(hypotheticalpat, infoprot_blast['descr']) and float(perc_id) >= float(args.idthreshold) and float(perc_cover) >= float(args.covthreshold) and float(row['evalue']) <= float(args.blastevalue): | |
| 429 information_proteins_blast[row['qseqid']] = infoprot_blast | |
| 430 else: | |
| 431 continue | |
| 432 else: | |
| 433 if not re.search(hypotheticalpat, infoprot_blast['descr']) and float(perc_id) >= float(args.idthreshold) and float(perc_id) >= float(infoprot_blast['pident']) and float(perc_cover) >= float(args.covthreshold) and float(perc_cover) >= float(infoprot_blast['pcover']) and float(row['evalue']) <= float(args.blastevalue): | |
| 434 information_proteins_blast[row['qseqid']] = infoprot_blast | |
| 435 | |
| 436 ## Predicting the function of the proteins based on HMM predictions using phmmer | |
| 437 if args.nohmmer == False: | |
| 438 with open("commands.sh", "w") as commands: | |
| 439 for singleprot in sorted(glob.glob("SEQ_*.faa")): | |
| 440 hhmtable = "%s.tbl" % singleprot | |
| 441 eprint("Creating file to run parallel PHMMER") | |
| 442 eprint("Adding %s to run PHMMER." % singleprot) | |
| 443 line2write = ' '.join(["phmmer", "--cpu", "1", "--domtblout", hhmtable, "-E", str(args.hmmerevalue), "-o", "/dev/null", singleprot, args.hmmdatabase, '\n']) | |
| 444 commands.write(line2write) | |
| 445 | |
| 446 eprint("Running parallel PHMMER") | |
| 447 subprocess.call(['parallel', '-j', str(args.ncpus)], stdin=open('commands.sh')) | |
| 448 eprint("Done. PHMMER was done to predict the function of the genes according to Hidden Markov Models\n") | |
| 449 os.remove("commands.sh") | |
| 450 | |
| 451 # Parsing the results from HMMER | |
| 452 information_proteins_hmmer = {} | |
| 453 hypotheticalpat = re.compile(r'(((((?i)hypothetical)|(?i)uncharacteri(z|s)ed|(?i)predicted)) protein)|((?i)ORF|((?i)unnamed protein product|(?i)gp\d+))') | |
| 454 for singletbl in sorted(glob.glob("*.faa.tbl")): | |
| 455 rootname = singletbl.replace(".faa.tbl", "") | |
| 456 with open(singletbl) as tblfile: | |
| 457 infoprot_hmmer = {} | |
| 458 for line in tblfile: | |
| 459 if line.startswith("#"): | |
| 460 continue | |
| 461 else: | |
| 462 try: | |
| 463 pat = re.compile("^(\S+)\s+\S\s+\d+\s+(\S+)\s+\S\s+(\d+)\s+((?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)\s+\S+\s+\S+\s+\S+\s+\S+\s+(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?\s+\S+\s+\S+\s+\S+\s+(\d+)\s+(\d+)\s+\d+\s+\d+\s+\S+\s+\S+\s+(\S+)\s+(.*)") | |
| 464 matchname, lociname, length, evaluehh, start, end, pident, description = re.match(pat, line).groups() | |
| 465 except AttributeError: | |
| 466 continue | |
| 467 else: | |
| 468 length = float(length) | |
| 469 pident = 100.00*float(pident) | |
| 470 protarea = 100.00*(((float(end)-1.00) - (float(start)-1.00))/length) | |
| 471 try: | |
| 472 data2 = infoprot_hmmer['lociname'] | |
| 473 except KeyError: | |
| 474 if not re.search(hypotheticalpat, description) and float(protarea) >= float(args.covthreshold) and float(evaluehh) <= float(args.hmmerevalue) and float(pident) >= 50.00: | |
| 475 infoprot_hmmer['lociname'] = lociname | |
| 476 infoprot_hmmer['name'] = matchname | |
| 477 infoprot_hmmer['evalue'] = float(evaluehh) | |
| 478 infoprot_hmmer['pcover'] = float(protarea) | |
| 479 infoprot_hmmer['pident'] = float(pident) | |
| 480 infoprot_hmmer['descr'] = description | |
| 481 else: | |
| 482 try: | |
| 483 if not re.search(hypotheticalpat, description) and float(protarea) >= float(args.covthreshold) and float(evaluehh) <= float(args.hmmerevalue) and float(pident) >= 50.00 and float(pident) >= infoprot_hmmer['pident']: | |
| 484 infoprot_hmmer['lociname'] = lociname | |
| 485 infoprot_hmmer['name'] = matchname | |
| 486 infoprot_hmmer['evalue'] = float(evaluehh) | |
| 487 infoprot_hmmer['pcover'] = float(protarea) | |
| 488 infoprot_hmmer['pident'] = float(pident) | |
| 489 infoprot_hmmer['descr'] = description | |
| 490 except KeyError: | |
| 491 continue | |
| 492 else: | |
| 493 if not re.search(hypotheticalpat, description) and float(protarea) >= float(args.covthreshold) and float(evaluehh) <= float(args.hmmerevalue) and float(pident) >= 50.00 and float(pident) >= infoprot_hmmer['pident']: | |
| 494 infoprot_hmmer['lociname'] = lociname | |
| 495 infoprot_hmmer['name'] = matchname | |
| 496 infoprot_hmmer['evalue'] = float(evaluehh) | |
| 497 infoprot_hmmer['pcover'] = float(protarea) | |
| 498 infoprot_hmmer['pident'] = float(pident) | |
| 499 infoprot_hmmer['descr'] = description | |
| 500 information_proteins_hmmer[rootname] = infoprot_hmmer | |
| 501 | |
| 502 #Storing protein information in memory | |
| 503 with open("temp.faa", "rU") as protsfile: | |
| 504 tempprotsdict = {} | |
| 505 for protseq in SeqIO.parse(protsfile, "fasta"): | |
| 506 tempindprot = {} | |
| 507 dataprot = protseq.description.split(' # ') | |
| 508 modseq = str(protseq.seq).replace("X","") # Removing all ambiguous amino acids to avoid problems with ProteinAnalysis module | |
| 509 analysed_seq = ProteinAnalysis(modseq) | |
| 510 tempindprot['length'] = len(protseq.seq) | |
| 511 tempindprot['isoelectricpoint'] = analysed_seq.isoelectric_point() | |
| 512 tempindprot['molweightkda'] = analysed_seq.molecular_weight()/1000.00 | |
| 513 tempindprot['instability'] = analysed_seq.instability_index() | |
| 514 tempindprot['protein_id'] = dataprot[0] | |
| 515 tempindprot['strand'] = int(dataprot[3]) | |
| 516 tempindprot['begin'] = int(dataprot[1])-1 | |
| 517 tempindprot['end'] = int(dataprot[2]) | |
| 518 tempprotsdict[dataprot[0]] = tempindprot | |
| 519 | |
| 520 # Creation of table | |
| 521 debugtable = "%s.csv" % newfile | |
| 522 with open(debugtable, "w") as tablefile: | |
| 523 if args.nohmmer == False: | |
| 524 print("\t".join(["Identifier", "Start", "Stop", "Strand", "size_aa", "pI", "Mol_weight_kDa", "Instability_index", "ID_BLAST", "Descr_BLAST", "evalue_BLAST", "%ID_BLAST", "%Cover_BLAST", "ID_HMMER", "Descr_HMMER", "evalue_HMMER", "%ID_HMMER", "%Cover_HMMER"]), file=tablefile) | |
| 525 keylist = information_proteins_hmmer.keys() | |
| 526 keylist.sort() | |
| 527 for keyB in keylist: | |
| 528 keyB = keyB.replace(".faa.tbl", "") | |
| 529 try: | |
| 530 print("\t".join([equivalence[keyB], str(tempprotsdict[equivalence[keyB]]['begin']), str(tempprotsdict[equivalence[keyB]]['end']), str(tempprotsdict[equivalence[keyB]]['strand']), str(tempprotsdict[equivalence[keyB]]['length']), str(tempprotsdict[equivalence[keyB]]['isoelectricpoint']), str(tempprotsdict[equivalence[keyB]]['molweightkda']), str(tempprotsdict[equivalence[keyB]]['instability']), information_proteins_blast[equivalence[keyB]]['sseqid'], information_proteins_blast[equivalence[keyB]]['descr'], str(information_proteins_blast[equivalence[keyB]]['evalue']), str(information_proteins_blast[equivalence[keyB]]['pident']), str(information_proteins_blast[equivalence[keyB]]['pcover']), information_proteins_hmmer[keyB]['name'], information_proteins_hmmer[keyB]['descr'], str(information_proteins_hmmer[keyB]['evalue']), str(information_proteins_hmmer[keyB]['pident']), str(information_proteins_hmmer[keyB]['pcover'])]), file=tablefile) | |
| 531 except KeyError: | |
| 532 try: | |
| 533 print("\t".join([equivalence[keyB], str(tempprotsdict[equivalence[keyB]]['begin']), str(tempprotsdict[equivalence[keyB]]['end']), str(tempprotsdict[equivalence[keyB]]['strand']), str(tempprotsdict[equivalence[keyB]]['length']), str(tempprotsdict[equivalence[keyB]]['isoelectricpoint']), str(tempprotsdict[equivalence[keyB]]['molweightkda']), str(tempprotsdict[equivalence[keyB]]['instability']), information_proteins_blast[equivalence[keyB]]['sseqid'], information_proteins_blast[equivalence[keyB]]['descr'], str(information_proteins_blast[equivalence[keyB]]['evalue']), str(information_proteins_blast[equivalence[keyB]]['pident']), str(information_proteins_blast[equivalence[keyB]]['pcover']), "None", "None", "NA", "NA", "NA"]), file=tablefile) | |
| 534 except KeyError: | |
| 535 try: | |
| 536 print("\t".join([equivalence[keyB], str(tempprotsdict[equivalence[keyB]]['begin']), str(tempprotsdict[equivalence[keyB]]['end']), str(tempprotsdict[equivalence[keyB]]['strand']), str(tempprotsdict[equivalence[keyB]]['length']), str(tempprotsdict[equivalence[keyB]]['isoelectricpoint']), str(tempprotsdict[equivalence[keyB]]['molweightkda']), str(tempprotsdict[equivalence[keyB]]['instability']), "None", "None", "NA", "NA", "NA", information_proteins_hmmer[keyB]['name'], information_proteins_hmmer[keyB]['descr'], str(information_proteins_hmmer[keyB]['evalue']), str(information_proteins_hmmer[keyB]['pident']), str(information_proteins_hmmer[keyB]['pcover'])]), file=tablefile) | |
| 537 except KeyError: | |
| 538 print("\t".join([equivalence[keyB], str(tempprotsdict[equivalence[keyB]]['begin']), str(tempprotsdict[equivalence[keyB]]['end']), str(tempprotsdict[equivalence[keyB]]['strand']), str(tempprotsdict[equivalence[keyB]]['length']), str(tempprotsdict[equivalence[keyB]]['isoelectricpoint']), str(tempprotsdict[equivalence[keyB]]['molweightkda']), str(tempprotsdict[equivalence[keyB]]['instability']), "None", "None", "NA", "NA", "NA", "None", "None", "NA", "NA", "NA"]), file=tablefile) | |
| 539 else: | |
| 540 print("\t".join(["Identifier", "Start", "Stop", "Strand", "size_aa", "pI", "Mol_weight_kDa", "Instability_index", "ID_BLAST", "Descr_BLAST", "evalue_BLAST", "%ID_BLAST", "%Cover_BLAST"]), file=tablefile) | |
| 541 keylist = equivalence.values() | |
| 542 keylist.sort() | |
| 543 for keyB in keylist: | |
| 544 try: | |
| 545 print("\t".join([keyB, str(tempprotsdict[keyB]['begin']), str(tempprotsdict[keyB]['end']), str(tempprotsdict[keyB]['strand']), str(tempprotsdict[keyB]['length']), str(tempprotsdict[keyB]['isoelectricpoint']), str(tempprotsdict[keyB]['molweightkda']), str(tempprotsdict[keyB]['instability']), information_proteins_blast[keyB]['sseqid'], information_proteins_blast[keyB]['descr'], str(information_proteins_blast[keyB]['evalue']), str(information_proteins_blast[keyB]['pident']), str(information_proteins_blast[keyB]['pcover'])]), file=tablefile) | |
| 546 except KeyError: | |
| 547 print("\t".join([keyB, str(tempprotsdict[keyB]['begin']), str(tempprotsdict[keyB]['end']), str(tempprotsdict[keyB]['strand']), str(tempprotsdict[keyB]['length']), str(tempprotsdict[keyB]['isoelectricpoint']), str(tempprotsdict[keyB]['molweightkda']), str(tempprotsdict[keyB]['instability']), "None", "None", "NA", "NA", "NA"]), file=tablefile) | |
| 548 | |
| 549 # Algorithm of decisions (which one: BLAST/HMMER?) | |
| 550 multipleprots = {} | |
| 551 Hypotheticalpat = re.compile(r'(((((?i)hypothetical)|(?i)uncharacteri(z|s)ed|(?i)predicted)) protein)|((?i)ORF|((?i)unnamed protein product|(?i)gp\d+))') | |
| 552 if args.nohmmer == False: | |
| 553 keylist = information_proteins_hmmer.keys() | |
| 554 keylist.sort() | |
| 555 for keyB in keylist: | |
| 556 keyB = keyB.replace(".faa.tbl", "") | |
| 557 singleprot = {} | |
| 558 singleprot['name'] = equivalence[keyB] | |
| 559 if (equivalence[keyB] in information_proteins_blast) and (keyB in information_proteins_hmmer): | |
| 560 if re.search(Hypotheticalpat, information_proteins_blast[equivalence[keyB]]['descr']): | |
| 561 try: | |
| 562 if re.search(Hypotheticalpat, information_proteins_hmmer[keyB]['descr']): | |
| 563 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 564 else: | |
| 565 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 566 except KeyError: | |
| 567 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 568 else: | |
| 569 try: | |
| 570 if (float(information_proteins_blast[equivalence[keyB]]['pident'])>float(information_proteins_hmmer[keyB]['pident'])) and (float(information_proteins_blast[equivalence[keyB]]['pcover'])>float(information_proteins_hmmer[keyB]['pcover'])): | |
| 571 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 572 elif (float(information_proteins_blast[equivalence[keyB]]['pident'])<float(information_proteins_hmmer[keyB]['pident'])) and (float(information_proteins_blast[equivalence[keyB]]['pcover'])<float(information_proteins_hmmer[keyB]['pcover'])): | |
| 573 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 574 elif (float(information_proteins_blast[equivalence[keyB]]['pident'])>float(information_proteins_hmmer[keyB]['pident'])) and (float(information_proteins_blast[equivalence[keyB]]['pcover'])<float(information_proteins_hmmer[keyB]['pcover'])): | |
| 575 if (float(information_proteins_blast[equivalence[keyB]]['pident'])-float(information_proteins_hmmer[keyB]['pident']) >= args.diffid): | |
| 576 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 577 else: | |
| 578 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 579 else: | |
| 580 if (float(information_proteins_hmmer[keyB]['pident'])-float(information_proteins_blast[equivalence[keyB]]['pident']) >= args.diffid): | |
| 581 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 582 else: | |
| 583 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 584 except KeyError: | |
| 585 try: | |
| 586 if (float(information_proteins_blast[equivalence[keyB]]['pcover'])>float(information_proteins_hmmer[keyB]['pcover'])): | |
| 587 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 588 else: | |
| 589 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 590 except KeyError: | |
| 591 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 592 elif equivalence[keyB] in information_proteins_blast: | |
| 593 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 594 elif keyB in information_proteins_hmmer: | |
| 595 try: | |
| 596 singleprot['descr'] = information_proteins_hmmer[keyB]['descr'] | |
| 597 except KeyError: | |
| 598 singleprot['descr'] = 'Hypothetical protein' | |
| 599 else: | |
| 600 singleprot['descr'] = information_proteins_blast[equivalence[keyB]]['descr'] | |
| 601 multipleprots[keyB] = singleprot | |
| 602 else: | |
| 603 keylist = equivalence.values() | |
| 604 keylist.sort() | |
| 605 for keyB in keylist: | |
| 606 singleprot = {} | |
| 607 singleprot['name'] = keyB | |
| 608 try: | |
| 609 if information_proteins_blast[keyB]['descr'] == None: | |
| 610 singleprot['descr'] = 'Hypothetical protein' | |
| 611 elif re.search(Hypotheticalpat, information_proteins_blast[keyB]['descr']): | |
| 612 singleprot['descr'] = 'Conserved hypothetical protein' | |
| 613 else: | |
| 614 singleprot['descr'] = information_proteins_blast[keyB]['descr'] | |
| 615 except KeyError: | |
| 616 singleprot['descr'] = 'Hypothetical protein' | |
| 617 multipleprots[keyB] = singleprot | |
| 618 | |
| 619 #Storing protein information in memory | |
| 620 with open("temp.faa", "rU") as protsfile: | |
| 621 protsdict = {} | |
| 622 for protseq in SeqIO.parse(protsfile, "fasta"): | |
| 623 indprot = {} | |
| 624 dataprot = protseq.description.split(' # ') | |
| 625 indprot['translation'] = protseq.seq | |
| 626 indprot['protein_id'] = dataprot[0] | |
| 627 indprot['strand'] = int(dataprot[3]) | |
| 628 indprot['begin'] = int(dataprot[1])-1 | |
| 629 indprot['end'] = int(dataprot[2]) | |
| 630 for keyOmega in sorted(multipleprots): | |
| 631 if multipleprots[keyOmega]['name'] == dataprot[0]: | |
| 632 indprot['product'] = multipleprots[keyOmega]['descr'] | |
| 633 protsdict[dataprot[0]] = indprot | |
| 634 | |
| 635 # Predicting the rRNA sequences | |
| 636 with open(newfile, "rU") as targetfasta, open("/dev/null", "w") as apocalypse: | |
| 637 eprint("Running INFERNAL+RFAM to predict rRNA-like sequences in %s" % newfile) | |
| 638 subprocess.call(["cmscan", "--rfam", "--cut_ga", "--nohmmonly", "--tblout", "rrnafile.csv", "--cpu", args.ncpus, args.rfamdatabase, newfile], stdout=apocalypse) | |
| 639 | |
| 640 #Storing rRNA information in memory | |
| 641 with open("rrnafile.csv", "rU") as rrnafile: | |
| 642 namedict = {"SSU_rRNA_archaea": "16s_rRNA", "SSU_rRNA_bacteria": "16s_rRNA", "SSU_rRNA_eukarya": "18s_rRNA", "SSU_rRNA_microsporidia": "16s_rRNA", "LSU_rRNA_archaea": "23s_rRNA", "LSU_rRNA_bacteria": "23s_rRNA", "LSU_rRNA_eukarya": "28s_rRNA", "5S_rRNA": "5s_rRNA"} | |
| 643 rRNAdict = defaultdict(list) | |
| 644 for line in rrnafile: | |
| 645 if not line.startswith("#"): | |
| 646 InfoLINE = re.sub("\s+", ",", line) | |
| 647 line_splitted = InfoLINE.split(",") | |
| 648 item_type = line_splitted[0] | |
| 649 if item_type.startswith(('LSU', 'SSU', '5S')): | |
| 650 strand = 0 | |
| 651 if line_splitted[9] == "+": | |
| 652 strand = 1 | |
| 653 else: | |
| 654 strand = -1 | |
| 655 rRNAdict[item_type].append({'score': float(line_splitted[14]), 'product': namedict[line_splitted[0]], 'begin': int(line_splitted[7]), 'end': int(line_splitted[8]), 'strand': strand}) | |
| 656 | |
| 657 subunits = {'LSU': {'max_score': 0 }, 'SSU': {'max_score': 0 }, '5S': {'max_score': 0 }} | |
| 658 for type_rRNA, rRNA_data in rRNAdict.items(): | |
| 659 max_score = max([item['score'] for item in rRNAdict[type_rRNA]]) | |
| 660 for item in ('LSU', 'SSU'): | |
| 661 if type_rRNA.startswith(item): | |
| 662 if max_score > subunits[item]['max_score']: | |
| 663 subunits[item]['listdata'] = rRNA_data | |
| 664 subunits[item]['max_score'] = max_score | |
| 665 if type_rRNA.startswith('5S'): | |
| 666 subunits['5S']['listdata'] = rRNA_data | |
| 667 subunits['5S']['max_score'] = max_score | |
| 668 | |
| 669 for rRNA in subunits: | |
| 670 i = 0 | |
| 671 try: | |
| 672 lengthlist = len(subunits[rRNA]['listdata']) | |
| 673 except KeyError: | |
| 674 continue | |
| 675 else: | |
| 676 while i < lengthlist: | |
| 677 eprint("%s harbours a %s from %i to %i" % (newfile, subunits[rRNA]['listdata'][i]['product'], int(subunits[rRNA]['listdata'][i]['begin']), int(subunits[rRNA]['listdata'][i]['end']))) | |
| 678 i += 1 | |
| 679 | |
| 680 # Predicting the tRNA sequences | |
| 681 eprint("Running ARAGORN to predict tRNA-like sequences in %s" % newfile) | |
| 682 genetictable = "-gc%s" % str(args.gcode) | |
| 683 with open("trnafile.fasta", "w") as trnafile: | |
| 684 if genomeshape['genomeshape'] == "circular": | |
| 685 subprocess.call(["aragorn", "-c", "-fon", genetictable, newfile], stdout=trnafile) | |
| 686 else: | |
| 687 subprocess.call(["aragorn", "-l", "-fon", genetictable, newfile], stdout=trnafile) | |
| 688 num_tRNA = len(list(SeqIO.parse("trnafile.fasta", "fasta"))) | |
| 689 eprint("ARAGORN was able to predict %i tRNAs in %s\n" % (num_tRNA, newfile)) | |
| 690 | |
| 691 #Storing tRNA and tmRNA information in memory | |
| 692 with open("trnafile.fasta", "rU") as trnafile: | |
| 693 tRNAdict = {} | |
| 694 tmRNAdict = {} | |
| 695 for tRNAseq in SeqIO.parse(trnafile, "fasta"): | |
| 696 indtRNA = {} | |
| 697 indtmRNA = {} | |
| 698 tRNA_information = tRNAseq.description.split(" ") | |
| 699 tRNApat = re.compile("^tRNA-") | |
| 700 if tRNA_information[1] == "tmRNA": | |
| 701 if str(tRNA_information[2]) == "(Permuted)": | |
| 702 indtmRNA['product'] = "tmRNA" | |
| 703 tmRNA_coords = str(tRNA_information[3]) | |
| 704 Beginningrevcomppat = re.compile("^c") | |
| 705 if re.match(Beginningrevcomppat, tmRNA_coords): | |
| 706 indtmRNA['strand'] = -1 | |
| 707 tmRNA_coords = tmRNA_coords.replace("c[","").replace("]","").split(",") | |
| 708 else: | |
| 709 indtmRNA['strand'] = 1 | |
| 710 tmRNA_coords = tmRNA_coords.replace("[","").replace("]","").split(",") | |
| 711 indtmRNA['begin'] = int(tmRNA_coords[0]) | |
| 712 indtmRNA['end'] = int(tmRNA_coords[1]) | |
| 713 tmRNAdict[tRNAseq.id] = indtmRNA | |
| 714 else: | |
| 715 indtmRNA['product'] = "tmRNA" | |
| 716 tmRNA_coords = str(tRNA_information[2]) | |
| 717 Beginningrevcomppat = re.compile("^c") | |
| 718 if re.match(Beginningrevcomppat, tmRNA_coords): | |
| 719 indtmRNA['strand'] = -1 | |
| 720 tmRNA_coords = tmRNA_coords.replace("c[","").replace("]","").split(",") | |
| 721 else: | |
| 722 indtmRNA['strand'] = 1 | |
| 723 tmRNA_coords = tmRNA_coords.replace("[","").replace("]","").split(",") | |
| 724 indtmRNA['begin'] = int(tmRNA_coords[0]) | |
| 725 indtmRNA['end'] = int(tmRNA_coords[1]) | |
| 726 tmRNAdict[tRNAseq.id] = indtmRNA | |
| 727 elif re.match(tRNApat, tRNA_information[1]): | |
| 728 indtRNA['product'] = re.sub("\(\w{3}\)", "", tRNA_information[1]) | |
| 729 tRNA_coords = tRNA_information[2] | |
| 730 Beginningrevcomppat = re.compile("^c") | |
| 731 if re.match(Beginningrevcomppat, tRNA_coords): | |
| 732 indtRNA['strand'] = -1 | |
| 733 tRNA_coords = tRNA_coords.replace("c[","").replace("]","").split(",") | |
| 734 else: | |
| 735 indtRNA['strand'] = 1 | |
| 736 tRNA_coords = tRNA_coords.replace("[","").replace("]","").split(",") | |
| 737 indtRNA['begin'] = int(tRNA_coords[0]) | |
| 738 indtRNA['end'] = int(tRNA_coords[1]) | |
| 739 tRNAdict[tRNAseq.id] = indtRNA | |
| 740 | |
| 741 #Predicting CRISPR repeats and others | |
| 742 eprint("Running PILERCR to predict CRISPR repeats in %s" % newfile) | |
| 743 subprocess.call(["pilercr", "-in", newfile, "-out", "crisprfile.txt", "-noinfo", "-minrepeat", str(args.minrepeat), "-maxrepeat", str(args.maxrepeat), "-minspacer", str(args.minspacer), "-maxspacer", str(args.maxspacer)]) | |
| 744 eprint("Predicting repeats in the sequences using TRF and IRF") | |
| 745 with open("/dev/null", "w") as stderr: | |
| 746 subprocess.call(["trf", newfile, "2", "7", "7", "80", "10", "50", "500", "-h"], stderr=stderr) | |
| 747 os.rename("%s.2.7.7.80.10.50.500.dat" % newfile, "trf_temp.dat") | |
| 748 with open("/dev/null", "w") as stderr: | |
| 749 subprocess.call(["irf", newfile, "2", "3", "5", "80", "10", "40", "500000", "10000", "-d", "-h"], stderr=stderr) | |
| 750 os.rename("%s.2.3.5.80.10.40.500000.10000.dat" % newfile, "irf_temp.dat") | |
| 751 | |
| 752 # Storing CRISPR repeats information | |
| 753 information_CRISPR = {} | |
| 754 with open("crisprfile.txt", "rU") as crisprfile: | |
| 755 for line in crisprfile: | |
| 756 if "SUMMARY BY POSITION" in line: | |
| 757 for line in crisprfile: | |
| 758 information_crispr_repeat = {} | |
| 759 try: | |
| 760 patC = re.compile('^\s+(\d+)\s+.{16}\s+(\d+)\s+(\d+)\s+\d+\s+\d+\s+\d+\s+\d?\s+(\w+)') | |
| 761 key, start, length, seq = re.match(patC, line).groups() | |
| 762 except AttributeError: | |
| 763 continue | |
| 764 else: | |
| 765 information_crispr_repeat['start'] = start | |
| 766 information_crispr_repeat['end'] = int(start) + int(length) | |
| 767 information_crispr_repeat['repeatseq'] = seq | |
| 768 information_crispr_repeat['repeatend'] = int(start) + len(seq) | |
| 769 information_CRISPR[key] = information_crispr_repeat | |
| 770 | |
| 771 # Storing tandem repeats information | |
| 772 information_TRF = {} | |
| 773 count = 1 | |
| 774 with open("trf_temp.dat", "rU") as trfile: | |
| 775 for line in trfile: | |
| 776 information_tandem_repeat = {} | |
| 777 try: | |
| 778 patT = re.compile('^(\d+)\s(\d+)\s\d+\s\d+\.\d+\s') | |
| 779 start, end = re.match(patT, line).groups() | |
| 780 except AttributeError: | |
| 781 continue | |
| 782 else: | |
| 783 information_tandem_repeat['start'] = start | |
| 784 information_tandem_repeat['end'] = end | |
| 785 information_TRF[count] = information_tandem_repeat | |
| 786 count += 1 | |
| 787 | |
| 788 # Storing inverted repeats information | |
| 789 information_IRF = {} | |
| 790 count = 1 | |
| 791 with open("irf_temp.dat", "rU") as irfile: | |
| 792 for line in irfile: | |
| 793 information_inverted_repeat = {} | |
| 794 try: | |
| 795 patI = re.compile('^(\d+)\s(\d+)\s\d+\s\d+\s\d+') | |
| 796 start, end = re.match(patI, line).groups() | |
| 797 except AttributeError: | |
| 798 continue | |
| 799 else: | |
| 800 information_inverted_repeat['start'] = start | |
| 801 information_inverted_repeat['end'] = end | |
| 802 information_IRF[count] = information_inverted_repeat | |
| 803 count += 1 | |
| 804 | |
| 805 # Creating a new Genbank and GFF file | |
| 806 eprint("Creating the output files") | |
| 807 newtempgbk = "%s.temp.gbk" % newfile | |
| 808 with open(newfile, "rU") as basefile, open(newtempgbk, "w"): | |
| 809 for record in SeqIO.parse(basefile, "fasta", IUPAC.ambiguous_dna): | |
| 810 whole_sequence = SeqRecord(record.seq) | |
| 811 whole_sequence.id = str(record.id) | |
| 812 whole_sequence.annotations['data_file_division'] = args.typedata.upper() | |
| 813 whole_sequence.annotations['date'] = strftime("%d-%b-%Y").upper() | |
| 814 for protein in sorted(protsdict, key = stringSplitByNumbers): | |
| 815 putative_start = int(protsdict[protein]['begin']) | |
| 816 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 817 end_pos = SeqFeature.ExactPosition(protsdict[protein]['end']) | |
| 818 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos, strand=protsdict[protein]['strand']) | |
| 819 qualifiersgene = OrderedDict([('locus_tag', protsdict[protein]['protein_id'])]) | |
| 820 new_data_gene = SeqFeature.SeqFeature(feature_location, type = "gene", strand = protsdict[protein]['strand'], qualifiers = qualifiersgene) | |
| 821 whole_sequence.features.append(new_data_gene) | |
| 822 qualifiers = [('locus_tag', protsdict[protein]['protein_id']), ('product', protsdict[protein]['product']), ('protein_id', protsdict[protein]['protein_id']), ('translation', protsdict[protein]['translation'])] | |
| 823 feature_qualifiers = OrderedDict(qualifiers) | |
| 824 new_data_cds = SeqFeature.SeqFeature(feature_location, type = "CDS", strand = protsdict[protein]['strand'], qualifiers = feature_qualifiers) | |
| 825 whole_sequence.features.append(new_data_cds) | |
| 826 for rRNA in sorted(subunits, key = stringSplitByNumbers): | |
| 827 i = 0 | |
| 828 try: | |
| 829 lengthlist = len(subunits[rRNA]['listdata']) | |
| 830 except KeyError: | |
| 831 continue | |
| 832 else: | |
| 833 while i < lengthlist: | |
| 834 putative_start = int(subunits[rRNA]['listdata'][i]['begin']) | |
| 835 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 836 end_pos = SeqFeature.ExactPosition(subunits[rRNA]['listdata'][i]['end']) | |
| 837 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos, strand=subunits[rRNA]['listdata'][i]['strand']) | |
| 838 new_data_gene = SeqFeature.SeqFeature(feature_location, type = "gene", strand = subunits[rRNA]['listdata'][i]['strand']) | |
| 839 whole_sequence.features.append(new_data_gene) | |
| 840 qualifiers = [('product', subunits[rRNA]['listdata'][i]['product'])] | |
| 841 feature_qualifiers = OrderedDict(qualifiers) | |
| 842 new_data_rRNA = SeqFeature.SeqFeature(feature_location, type = "rRNA", strand = subunits[rRNA]['listdata'][i]['strand'], qualifiers = feature_qualifiers) | |
| 843 whole_sequence.features.append(new_data_rRNA) | |
| 844 i += 1 | |
| 845 for tRNA in sorted(tRNAdict, key = stringSplitByNumbers): | |
| 846 putative_start = int(tRNAdict[tRNA]['begin']) | |
| 847 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 848 end_pos = SeqFeature.ExactPosition(tRNAdict[tRNA]['end']) | |
| 849 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos, strand=tRNAdict[tRNA]['strand']) | |
| 850 new_data_gene = SeqFeature.SeqFeature(feature_location, type = "gene", strand = tRNAdict[tRNA]['strand']) | |
| 851 whole_sequence.features.append(new_data_gene) | |
| 852 qualifiers = [('product', tRNAdict[tRNA]['product'])] | |
| 853 feature_qualifiers = OrderedDict(qualifiers) | |
| 854 new_data_tRNA = SeqFeature.SeqFeature(feature_location, type = "tRNA", strand = tRNAdict[tRNA]['strand'], qualifiers = feature_qualifiers) | |
| 855 whole_sequence.features.append(new_data_tRNA) | |
| 856 for tmRNA in sorted(tmRNAdict, key = stringSplitByNumbers): | |
| 857 putative_start = int(tmRNAdict[tmRNA]['begin']) | |
| 858 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 859 end_pos = SeqFeature.ExactPosition(tmRNAdict[tmRNA]['end']) | |
| 860 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos, strand=tmRNAdict[tmRNA]['strand']) | |
| 861 new_data_gene = SeqFeature.SeqFeature(feature_location, type = "gene", strand = tmRNAdict[tmRNA]['strand']) | |
| 862 whole_sequence.features.append(new_data_gene) | |
| 863 qualifiers = [('product', tmRNAdict[tmRNA]['product'])] | |
| 864 feature_qualifiers = OrderedDict(qualifiers) | |
| 865 new_data_tmRNA = SeqFeature.SeqFeature(feature_location, type = "tmRNA", strand = tmRNAdict[tmRNA]['strand'], qualifiers = feature_qualifiers) | |
| 866 whole_sequence.features.append(new_data_tmRNA) | |
| 867 for CRISPR in sorted(information_CRISPR, key = stringSplitByNumbers): | |
| 868 putative_start = int(information_CRISPR[CRISPR]['start']) | |
| 869 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 870 end_pos = SeqFeature.ExactPosition(information_CRISPR[CRISPR]['end']) | |
| 871 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos) | |
| 872 qualifiers = [('rpt_family', 'CRISPR'), ('rpt_type', 'direct'), ('rpt_unit_range', "%i..%i" % (int(information_CRISPR[CRISPR]['start']), int(information_CRISPR[CRISPR]['repeatend']))), ('rpt_unit_seq', information_CRISPR[CRISPR]['repeatseq'])] | |
| 873 feature_qualifiers = OrderedDict(qualifiers) | |
| 874 new_data_CRISPRrepeat = SeqFeature.SeqFeature(feature_location, type = "repeat_region", qualifiers = feature_qualifiers) | |
| 875 whole_sequence.features.append(new_data_CRISPRrepeat) | |
| 876 for tandem in sorted(information_TRF): | |
| 877 putative_start = int(information_TRF[tandem]['start']) | |
| 878 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 879 end_pos = SeqFeature.ExactPosition(information_TRF[tandem]['end']) | |
| 880 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos) | |
| 881 qualifiers = [('rpt_type', 'direct')] | |
| 882 feature_qualifiers = OrderedDict(qualifiers) | |
| 883 new_data_tandemrepeat = SeqFeature.SeqFeature(feature_location, type = "repeat_region", qualifiers = feature_qualifiers) | |
| 884 whole_sequence.features.append(new_data_tandemrepeat) | |
| 885 for inverted in sorted(information_IRF): | |
| 886 putative_start = int(information_IRF[inverted]['start']) | |
| 887 start_pos = SeqFeature.ExactPosition(putative_start) | |
| 888 end_pos = SeqFeature.ExactPosition(information_IRF[inverted]['end']) | |
| 889 feature_location = SeqFeature.FeatureLocation(start_pos, end_pos) | |
| 890 qualifiers = [('rpt_type', 'inverted')] | |
| 891 feature_qualifiers = OrderedDict(qualifiers) | |
| 892 new_data_invertedrepeat = SeqFeature.SeqFeature(feature_location, type = "repeat_region", qualifiers = feature_qualifiers) | |
| 893 whole_sequence.features.append(new_data_invertedrepeat) | |
| 894 SeqIO.write(whole_sequence, newtempgbk, "genbank") | |
| 895 | |
| 896 newgbk = "%s.gbk" % newfile | |
| 897 with open(newtempgbk, "rU") as gbktempfile, open(newgbk, "w") as gbkrealfile: | |
| 898 newpat = re.compile("D|RNA\s+(CON|PHG|VRL|BCT)") | |
| 899 for line in gbktempfile: | |
| 900 if line.startswith("LOCUS ") and re.search(newpat, line): | |
| 901 if genomeshape['genomeshape'] == "linear": | |
| 902 newline = re.sub("bp DNA\s+", "bp DNA linear ", line) | |
| 903 else: | |
| 904 newline = re.sub("bp DNA\s+", "bp DNA circular ", line) | |
| 905 gbkrealfile.write(newline) | |
| 906 else: | |
| 907 gbkrealfile.write(line) | |
| 908 | |
| 909 for f in glob.glob("*.temp.gbk"): | |
| 910 os.remove(f) | |
| 911 | |
| 912 if args.gffprint==True: | |
| 913 newgff = "%s.gff" % newfile | |
| 914 with open(newgff, "w") as outgff, open(newgbk, "rU") as ingbk: | |
| 915 GFF.write(SeqIO.parse(ingbk, "genbank"), outgff) | |
| 916 | |
| 917 # Removing intermediate files | |
| 918 os.remove(newfile) | |
| 919 os.remove("temp.faa") | |
| 920 os.remove("temp_blast.csv") | |
| 921 os.remove("crisprfile.txt") | |
| 922 os.remove("trnafile.fasta") | |
| 923 os.remove("rrnafile.csv") | |
| 924 os.remove("trf_temp.dat") | |
| 925 os.remove("irf_temp.dat") | |
| 926 for f in glob.glob("SEQ*"): | |
| 927 os.remove(f) | |
| 928 | |
| 929 # Joining all GENBANK files into one | |
| 930 listgbk = sorted(glob.glob('CONTIG_*.gbk')) | |
| 931 gbkoutputfile = "%s.gbk" % root_output | |
| 932 with open(gbkoutputfile, 'w') as finalgbk: | |
| 933 for fname in listgbk: | |
| 934 with open(fname) as infile: | |
| 935 for line in infile: | |
| 936 finalgbk.write(line) | |
| 937 | |
| 938 for tempgbk in glob.glob("CONTIG_*.gbk"): | |
| 939 os.remove(tempgbk) | |
| 940 | |
| 941 # Joining all GFF files into one | |
| 942 if args.gffprint==True: | |
| 943 listgff = sorted(glob.glob('CONTIG_*.gff')) | |
| 944 gffoutputfile = "%s.gff" % root_output | |
| 945 with open(gffoutputfile, 'w') as finalgff: | |
| 946 for fname in listgff: | |
| 947 with open(fname) as infile: | |
| 948 for line in infile: | |
| 949 finalgff.write(line) | |
| 950 for tempgff in glob.glob("CONTIG_*.gff"): | |
| 951 os.remove(tempgff) | |
| 952 | |
| 953 # Joining all TABLE files into one | |
| 954 listcsv = sorted(glob.glob('CONTIG_*.csv')) | |
| 955 tbloutputfile = "%s.csv" % root_output | |
| 956 with open(tbloutputfile, 'w') as finaltable: | |
| 957 for fname in listcsv: | |
| 958 with open(fname) as infile: | |
| 959 for line in infile: | |
| 960 finaltable.write(line) | |
| 961 | |
| 962 for temptbl in glob.glob("CONTIG_*.csv"): | |
| 963 os.remove(temptbl) | |
| 964 | |
| 965 # Preparing sequences for GenBank submission (Original code from Wan Yu's gbk2tbl.py script [https://github.com/wanyuac/BINF_toolkit/blob/master/gbk2tbl.py]) | |
| 966 allowed_qualifiers = ['locus_tag', 'gene', 'product', 'pseudo', 'protein_id', 'gene_desc', 'old_locus_tag', 'note', 'inference', 'organism', 'mol_type', 'strain', 'sub_species', 'isolation-source', 'country'] | |
| 967 newfastafile = "%s.fasta" % root_output | |
| 968 newtablefile = "%s.tbl" % root_output | |
| 969 with open(args.modifiers, "rU") as modifiers, open(gbkoutputfile, "r") as genbank_fh, open(newfastafile, "w") as fasta_fh, open(newtablefile, "w") as feature_fh: | |
| 970 info = modifiers.readline() | |
| 971 wholelist = list(SeqIO.parse(genbank_fh, 'genbank')) | |
| 972 for record in wholelist: | |
| 973 if len(record) <= args.mincontigsize: | |
| 974 eprint("WARNING: Skipping small contig %s" % record.id) | |
| 975 continue | |
| 976 record.description = "%s %s" % (record.id, info) | |
| 977 SeqIO.write([record], fasta_fh, 'fasta') | |
| 978 print('>Feature %s' % (record.name), file=feature_fh) | |
| 979 for line in record.features: | |
| 980 if line.strand == 1: | |
| 981 print('%d\t%d\t%s' % (line.location.nofuzzy_start + 1, line.location.nofuzzy_end, line.type), file=feature_fh) | |
| 982 else: | |
| 983 print('%d\t%d\t%s' % (line.location.nofuzzy_end, line.location.nofuzzy_start + 1, line.type), file=feature_fh) | |
| 984 for (key, values) in line.qualifiers.iteritems(): | |
| 985 if key not in allowed_qualifiers: | |
| 986 continue | |
| 987 for v in values: | |
| 988 print('\t\t\t%s\t%s' % (key, v), file=feature_fh) | |
| 989 | |
| 990 # Final statement | |
| 991 eprint("Genome annotation done!") | |
| 992 eprint("The GenBank file is %s" % gbkoutputfile) | |
| 993 if args.gffprint==True: | |
| 994 eprint("The GFF3 file is %s" % gffoutputfile) | |
| 995 eprint("The table file for GenBank submission is %s" % tbloutputfile) | |
| 996 eprint("The FASTA file for GenBank submission is %s" % newfastafile) | |
| 997 eprint("The table file with all protein information is %s" % newtablefile) |