[med-svn] [Git][med-team/python-cgecore][master] 8 commits: New upstream version 1.5.0
Andreas Tille
gitlab at salsa.debian.org
Thu Aug 1 21:29:53 BST 2019
Andreas Tille pushed to branch master at Debian Med / python-cgecore
Commits:
78248d1c by Andreas Tille at 2019-08-01T20:23:50Z
New upstream version 1.5.0
- - - - -
7031dbdc by Andreas Tille at 2019-08-01T20:23:53Z
Update upstream source from tag 'upstream/1.5.0'
Update to upstream version '1.5.0'
with Debian dir 43db26c29429410f0ccf58661160c3126e2074d2
- - - - -
996fd1f5 by Andreas Tille at 2019-08-01T20:23:53Z
New upstream version
- - - - -
2619a12b by Andreas Tille at 2019-08-01T20:23:53Z
debhelper-compat 12
- - - - -
b3f8e504 by Andreas Tille at 2019-08-01T20:23:56Z
Standards-Version: 4.4.0
- - - - -
f9588477 by Andreas Tille at 2019-08-01T20:23:56Z
Testsuite: autopkgtest-pkg-python
- - - - -
b8e6d542 by Andreas Tille at 2019-08-01T20:24:02Z
Set upstream metadata fields: Name.
Fixes lintian: upstream-metadata-file-is-missing
See https://lintian.debian.org/tags/upstream-metadata-file-is-missing.html for more details.
- - - - -
8cc4a6e9 by Andreas Tille at 2019-08-01T20:26:30Z
Upload to unstable
- - - - -
12 changed files:
- PKG-INFO
- README.md
- cgecore.egg-info/PKG-INFO
- cgecore/__init__.py
- cgecore/blaster/blaster.py
- cgecore/blaster/run_blaster.py
- cgecore/cgefinder.py
- cgecore/cmdline.py
- debian/changelog
- − debian/compat
- debian/control
- + debian/upstream/metadata
Changes:
=====================================
PKG-INFO
=====================================
@@ -1,6 +1,6 @@
Metadata-Version: 1.0
Name: cgecore
-Version: 1.4.4
+Version: 1.5.0
Summary: Center for Genomic Epidemiology Core Module
Home-page: https://bitbucket.org/genomicepidemiology/cge_core_module
Author: Center for Genomic Epidemiology
@@ -24,10 +24,12 @@ Description: # cge_core_module
# Install package locally
python2 setup.py install
+
python3 setup.py install
# Distribute to PyPi
python3 setup.py sdist bdist_wheel
+
twine upload dist/*
*deprecated:*
=====================================
README.md
=====================================
@@ -1,7 +1,7 @@
# cge_core_module
-Core module for the Center for Genomic Epidemiology
-
+Core module for the Center for Genomic Epidemiology
+
This module contains classes and functions needed to run the service wrappers and pipeline scripts
The pypi project can be found here:
@@ -16,10 +16,12 @@ https://pypi.org/project/cgecore/
# Install package locally
python2 setup.py install
+
python3 setup.py install
# Distribute to PyPi
python3 setup.py sdist bdist_wheel
+
twine upload dist/*
*deprecated:*
=====================================
cgecore.egg-info/PKG-INFO
=====================================
@@ -1,6 +1,6 @@
Metadata-Version: 1.0
Name: cgecore
-Version: 1.4.4
+Version: 1.5.0
Summary: Center for Genomic Epidemiology Core Module
Home-page: https://bitbucket.org/genomicepidemiology/cge_core_module
Author: Center for Genomic Epidemiology
@@ -24,10 +24,12 @@ Description: # cge_core_module
# Install package locally
python2 setup.py install
+
python3 setup.py install
# Distribute to PyPi
python3 setup.py sdist bdist_wheel
+
twine upload dist/*
*deprecated:*
=====================================
cgecore/__init__.py
=====================================
@@ -13,7 +13,7 @@ from .argumentparsing import (check_file_type, get_arguments, get_string,
)
#####################
-__version__ = "1.4.4"
+__version__ = "1.5.0"
__all__ = [
"argumentparsing",
"cmdline",
=====================================
cgecore/blaster/blaster.py
=====================================
@@ -14,588 +14,617 @@ import collections
class Blaster():
- def __init__(self, inputfile, databases, db_path, out_path='', min_cov=0.6,
- threshold=0.9, blast='blastn', cut_off=True,
- max_target_seqs=50000, reuse_results=False,
- allowed_overlap=0):
-
- min_cov = 100 * float(min_cov)
- threshold = 100 * float(threshold)
-
- # For alignment data storage
- self.gene_align_query = dict() # Sequence alignment lines
- self.gene_align_homo = dict() # Sequence alignment homolog string
- self.gene_align_sbjct = dict() # Sequence alignment allele string
- self.results = dict() # Results
-
- # TODO: Add excluded results to this dictionay
- self.results["excluded"] = dict()
-
- for db in databases:
- # Adding the path to the database and output
- db_file = "%s/%s.fsa" % (db_path, db)
- tmp_out_path = "%s/tmp" % (out_path)
- out_file = "%s/out_%s.xml" % (tmp_out_path, db)
- os.makedirs(tmp_out_path, exist_ok=True)
- os.chmod(tmp_out_path, 0o775)
-
- # Running blast
- if (os.path.isfile(out_file) and os.access(out_file, os.R_OK)
- and reuse_results):
- print("Found " + out_file + " skipping DB.")
- out, err = (b'', b'')
- else:
- cmd = ("%s -subject %s -query %s -out %s -outfmt '5'"
- " -perc_identity %s -max_target_seqs %s"
- " -dust 'no'" % (blast, db_file, inputfile,
- out_file, threshold, max_target_seqs))
-
- process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE,
- stderr=subprocess.PIPE)
- out, err = process.communicate()
-
- # Get results file
- try:
- # Test if output file exist
- result_handle = open(out_file, "r")
- except IOError:
- sys.exit(("Error: BLAST did not run as expected. " +
- "The expected output file, {}, did not exist.\n" +
- "BLAST finished with the following response:" +
- "\n{}\n{}").format(os.path.abspath(out_file),
- out.decode("utf-8"),
- err.decode("utf-8")))
-
- # Test if blast output is empty
- if os.stat(out_file).st_size == 0:
- sys.exit(("Error: BLAST did not run as expected. " +
- "The output file {} was empty.\n" +
- "BLAST finished with the following response:" +
- "\n{}\n{}").format(os.path.abspath(out_file),
- out.decode("utf-8"),
- err.decode("utf-8")))
-
- # Get blast output
- blast_records = NCBIXML.parse(result_handle)
-
- # Declaring variables for saving the results
-
- # results for each gene
- gene_results = dict()
-
- # For finding the best hits
- best_hsp = dict()
-
- # Keeping track of gene split
- gene_split = collections.defaultdict(dict)
-
- # Making the dicts for sequence outputs
- self.gene_align_query[db] = dict()
- self.gene_align_homo[db] = dict()
- self.gene_align_sbjct[db] = dict()
-
- # Parsing over the hits and only keeping the best
- for blast_record in blast_records:
- query = blast_record.query
- blast_record.alignments.sort(key=lambda align: (
- -max((len(hsp.query)
- * (int(hsp.identities) / float(len(hsp.query)))
- for hsp in align.hsps)))
- )
-
- for alignment in blast_record.alignments:
- # Setting the e-value as 1 and bit as 0 to get the best
- # HSP fragment
- best_e_value = 1
- best_bit = 0
-
- for hsp in alignment.hsps:
-
- if hsp.expect < best_e_value or hsp.bits > best_bit:
- best_e_value = hsp.expect
- best_bit = hsp.bits
-
- tmp = alignment.title.split(" ")
- sbjct_header = tmp[1]
- bit = hsp.bits
- sbjct_length = alignment.length
- sbjct_start = hsp.sbjct_start
- sbjct_end = hsp.sbjct_end
- gaps = hsp.gaps
- query_string = str(hsp.query)
- homo_string = str(hsp.match)
- sbjct_string = str(hsp.sbjct)
- contig_name = query.replace(">", "")
- query_start = hsp.query_start
- query_end = hsp.query_end
- HSP_length = len(query_string)
- perc_ident = int(hsp.identities) / float(HSP_length) * 100
- strand = 0
- coverage = ((int(HSP_length) - int(gaps))
- / float(sbjct_length))
- perc_coverage = (((int(HSP_length) - int(gaps))
- / float(sbjct_length)) * 100)
-
- # cal_score is later used to select the best hit
- cal_score = perc_ident * coverage
-
- hit_id = "%s:%s..%s:%s:%f" % (contig_name, query_start,
- query_end, sbjct_header,
- cal_score)
-
- # If the hit is on the other strand
- if sbjct_start > sbjct_end:
- tmp = sbjct_start
- sbjct_start = sbjct_end
- sbjct_end = tmp
-
- query_string = self.reversecomplement(query_string)
- homo_string = homo_string[::-1]
- sbjct_string = self.reversecomplement(sbjct_string)
- strand = 1
-
- # Save hit
- if (cut_off and perc_coverage > 20) or cut_off == False:
- best_hsp = {'evalue': hsp.expect,
- 'sbjct_header': sbjct_header,
- 'bit': bit,
- 'perc_ident': perc_ident,
- 'sbjct_length': sbjct_length,
- 'sbjct_start': sbjct_start,
- 'sbjct_end': sbjct_end,
- 'gaps': gaps,
- 'query_string': query_string,
- 'homo_string': homo_string,
- 'sbjct_string': sbjct_string,
- 'contig_name': contig_name,
- 'query_start': query_start,
- 'query_end': query_end,
- 'HSP_length': HSP_length,
- 'coverage': coverage,
- 'cal_score': cal_score,
- 'hit_id': hit_id,
- 'strand': strand,
- 'perc_coverage': perc_coverage
- }
-
- # Saving the result if any
- if best_hsp:
- save = 1
-
- # If there are other gene alignments they are compared
- if gene_results:
- tmp_gene_split = gene_split
- tmp_results = gene_results
-
- # Compare the hit results
- save, gene_split, gene_results = self.compare_results(
- save, best_hsp, tmp_results, tmp_gene_split,
- allowed_overlap)
-
- # If the hit is not overlapping with other hit
- # seqeunces it is kept
- if save == 1:
- gene_results[hit_id] = best_hsp
-
- result_handle.close()
-
- # If the hit does not cover the entire database reference the
- # missing seqence data are extracted
- keys = list(gene_results.keys())
-
- for hit_id in keys:
- hit = gene_results[hit_id]
-
- # Calculate possible split gene coverage
- perc_coverage = hit['perc_coverage']
-
- if(hit['sbjct_header'] in gene_split
- and len(gene_split[hit['sbjct_header']]) > 1):
-
- # Calculate new length
- new_length = self.calculate_new_length(gene_split, gene_results,
- hit)
- hit['split_length'] = new_length
-
- # Calculate new coverage
- perc_coverage = new_length / float(hit['sbjct_length']) * 100
-
- # If the hit is above the minimum length threshold it is kept
- if perc_coverage >= min_cov:
- if hit['coverage'] == 1:
- self.gene_align_query[db][hit_id] = hit['query_string']
- self.gene_align_homo[db][hit_id] = hit['homo_string']
- self.gene_align_sbjct[db][hit_id] = hit['sbjct_string']
- elif hit['coverage'] != 1:
- # Getting the whole database sequence
- for seq_record in SeqIO.parse(db_file, "fasta"):
- if seq_record.description.replace(" ", "") == hit['sbjct_header'].replace(" ", ""):
- start_seq = str(seq_record.seq)[:int(hit["sbjct_start"])-1]
- end_seq = str(seq_record.seq)[int(hit["sbjct_end"]):]
- self.gene_align_sbjct[db][hit_id] = start_seq + hit['sbjct_string'] + end_seq
- #self.gene_align_sbjct[db][hit_id] = str(seq_record.seq)
+ def __init__(self, inputfile, databases, db_path, out_path='', min_cov=0.6,
+ threshold=0.9, blast='blastn', cut_off=True,
+ max_target_seqs=50000, reuse_results=False,
+ allowed_overlap=0):
+
+ min_cov = 100 * float(min_cov)
+ threshold = 100 * float(threshold)
+
+ # For alignment data storage
+ self.gene_align_query = dict() # Sequence alignment lines
+ self.gene_align_homo = dict() # Sequence alignment homolog string
+ self.gene_align_sbjct = dict() # Sequence alignment allele string
+ self.results = dict() # Results
+
+ # TODO: Add excluded results to this dictionay
+ self.results["excluded"] = dict()
+
+ for db in databases:
+ # Adding the path to the database and output
+ db_file = "%s/%s.fsa" % (db_path, db)
+ tmp_out_path = "%s/tmp" % (out_path)
+ out_file = "%s/out_%s.xml" % (tmp_out_path, db)
+ os.makedirs(tmp_out_path, exist_ok=True)
+ os.chmod(tmp_out_path, 0o775)
+
+ # Running blast
+ if (os.path.isfile(out_file)
+ and os.access(out_file, os.R_OK)
+ and reuse_results):
+ print("Found " + out_file + " skipping DB.")
+ out, err = (b'', b'')
+ else:
+ cmd = ("%s -subject %s -query %s -out %s -outfmt '5'"
+ " -perc_identity %s -max_target_seqs %s"
+ " -dust 'no'" % (blast, db_file, inputfile,
+ out_file, threshold, max_target_seqs))
+
+ process = subprocess.Popen(cmd, shell=True,
+ stdout=subprocess.PIPE,
+ stderr=subprocess.PIPE)
+ out, err = process.communicate()
+
+ # Get results file
+ try:
+ # Test if output file exist
+ result_handle = open(out_file, "r")
+ except IOError:
+ sys.exit(("Error: BLAST did not run as expected. "
+ "The expected output file, {}, did not exist.\n"
+ "BLAST finished with the following response:"
+ "\n{}\n{}").format(os.path.abspath(out_file),
+ out.decode("utf-8"),
+ err.decode("utf-8")))
+
+ # Test if blast output is empty
+ if os.stat(out_file).st_size == 0:
+ sys.exit(("Error: BLAST did not run as expected. "
+ "The output file {} was empty.\n"
+ "BLAST finished with the following response:"
+ "\n{}\n{}").format(os.path.abspath(out_file),
+ out.decode("utf-8"),
+ err.decode("utf-8")))
+
+ # Get blast output
+ blast_records = NCBIXML.parse(result_handle)
+
+ # Declaring variables for saving the results
+
+ # results for each gene
+ gene_results = dict()
+
+ # For finding the best hits
+ best_hsp = dict()
+
+ # Keeping track of gene split
+ gene_split = collections.defaultdict(dict)
+
+ # Making the dicts for sequence outputs
+ self.gene_align_query[db] = dict()
+ self.gene_align_homo[db] = dict()
+ self.gene_align_sbjct[db] = dict()
+
+ # Parsing over the hits and only keeping the best
+ for blast_record in blast_records:
+
+ # OLD CODE TO BE REMOVED
+ # query = blast_record.query
+ # blast_record.alignments.sort(key=lambda align: (
+ # -max((len(hsp.query)
+ # * (int(hsp.identities) / float(len(hsp.query)))
+ # for hsp in align.hsps)))
+ # )
+
+ # Sort BLAST alignments by the hsp in each alignment with the
+ # highest number of identical nucleotides (hsp.identities)
+ blast_record.alignments.sort(
+ key=lambda align: (max((int(hsp.identities)
+ for hsp in align.hsps))),
+ reverse=True)
+
+ query = blast_record.query
+
+ for alignment in blast_record.alignments:
+ # Setting the e-value as 1 and bit as 0 to get the best
+ # HSP fragment
+ best_e_value = 1
+ best_bit = 0
+
+ for hsp in alignment.hsps:
+
+ if hsp.expect < best_e_value or hsp.bits > best_bit:
+ best_e_value = hsp.expect
+ best_bit = hsp.bits
+
+ tmp = alignment.title.split(" ")
+ sbjct_header = tmp[1]
+ # DEBUG
+ print("Found: {}".format(sbjct_header))
+ bit = hsp.bits
+ sbjct_length = alignment.length
+ sbjct_start = hsp.sbjct_start
+ sbjct_end = hsp.sbjct_end
+ gaps = hsp.gaps
+ query_string = str(hsp.query)
+ homo_string = str(hsp.match)
+ sbjct_string = str(hsp.sbjct)
+ contig_name = query.replace(">", "")
+ query_start = hsp.query_start
+ query_end = hsp.query_end
+ HSP_length = len(query_string)
+ perc_ident = (int(hsp.identities)
+ / float(HSP_length) * 100)
+ strand = 0
+ coverage = ((int(HSP_length) - int(gaps))
+ / float(sbjct_length))
+ perc_coverage = (((int(HSP_length) - int(gaps))
+ / float(sbjct_length)) * 100)
+
+ # cal_score is later used to select the best hit
+ cal_score = perc_ident * coverage
+
+ hit_id = "%s:%s..%s:%s:%f" % (
+ contig_name, query_start, query_end,
+ sbjct_header, cal_score)
+
+ # If the hit is on the other strand
+ if sbjct_start > sbjct_end:
+ tmp = sbjct_start
+ sbjct_start = sbjct_end
+ sbjct_end = tmp
+
+ query_string = self.reversecomplement(
+ query_string)
+ homo_string = homo_string[::-1]
+ sbjct_string = self.reversecomplement(
+ sbjct_string)
+ strand = 1
+
+ # Save hit
+ if((cut_off and perc_coverage > 20)
+ or cut_off is False):
+
+ best_hsp = {'evalue': hsp.expect,
+ 'sbjct_header': sbjct_header,
+ 'bit': bit,
+ 'perc_ident': perc_ident,
+ 'sbjct_length': sbjct_length,
+ 'sbjct_start': sbjct_start,
+ 'sbjct_end': sbjct_end,
+ 'gaps': gaps,
+ 'query_string': query_string,
+ 'homo_string': homo_string,
+ 'sbjct_string': sbjct_string,
+ 'contig_name': contig_name,
+ 'query_start': query_start,
+ 'query_end': query_end,
+ 'HSP_length': HSP_length,
+ 'coverage': coverage,
+ 'cal_score': cal_score,
+ 'hit_id': hit_id,
+ 'strand': strand,
+ 'perc_coverage': perc_coverage
+ }
+
+ # Saving the result if any
+ if best_hsp:
+ save = 1
+
+ # If there are other gene alignments they are compared
+ if gene_results:
+ tmp_gene_split = gene_split
+ tmp_results = gene_results
+
+ # Compare the hit results
+ save, gene_split, gene_results = (
+ self.compare_results(save, best_hsp,
+ tmp_results,
+ tmp_gene_split,
+ allowed_overlap)
+ )
+
+ # If the hit is not overlapping with other hit
+ # seqeunces it is kept
+ if save == 1:
+ # DEBUG
+ print("Saving: {}".format(hit_id))
+ gene_results[hit_id] = best_hsp
+
+ result_handle.close()
+
+ # If the hit does not cover the entire database reference the
+ # missing seqence data are extracted
+ keys = list(gene_results.keys())
+
+ for hit_id in keys:
+ hit = gene_results[hit_id]
+
+ # Calculate possible split gene coverage
+ perc_coverage = hit['perc_coverage']
+
+ if(hit['sbjct_header'] in gene_split
+ and len(gene_split[hit['sbjct_header']]) > 1):
+
+ # Calculate new length
+ new_length = self.calculate_new_length(gene_split, gene_results,
+ hit)
+ hit['split_length'] = new_length
+
+ # Calculate new coverage
+ perc_coverage = new_length / float(hit['sbjct_length']) * 100
+
+ # If the hit is above the minimum length threshold it is kept
+ if perc_coverage >= min_cov:
+ if hit['coverage'] == 1:
+ self.gene_align_query[db][hit_id] = hit['query_string']
+ self.gene_align_homo[db][hit_id] = hit['homo_string']
+ self.gene_align_sbjct[db][hit_id] = hit['sbjct_string']
+ elif hit['coverage'] != 1:
+ # Getting the whole database sequence
+ for seq_record in SeqIO.parse(db_file, "fasta"):
+ if seq_record.description.replace(" ", "") == hit['sbjct_header'].replace(" ", ""):
+ start_seq = str(seq_record.seq)[:int(hit["sbjct_start"])-1]
+ end_seq = str(seq_record.seq)[int(hit["sbjct_end"]):]
+ self.gene_align_sbjct[db][hit_id] = start_seq + hit['sbjct_string'] + end_seq
+ #self.gene_align_sbjct[db][hit_id] = str(seq_record.seq)
+ break
+
+ # Getting the whole contig to extract extra query seqeunce
+ contig = ''
+ for seq_record in SeqIO.parse(inputfile, "fasta"):
+ if seq_record.description.replace(" ", "") == hit['contig_name'].replace(" ", ""):
+ contig = str(seq_record.seq)
+ break
+
+ # Extract extra sequence from query
+ query_seq, homo_seq = self.get_query_align(hit, contig)
+
+ # Saving the new alignment sequences
+ self.gene_align_query[db][hit_id] = query_seq
+ self.gene_align_homo[db][hit_id] = homo_seq
+
+ else:
+ del gene_results[hit_id]
+ if hit['sbjct_header'] in gene_split:
+ del gene_split[hit['sbjct_header']]
+
+ # Save the database result
+ if gene_results:
+ self.results[db] = gene_results
+ else:
+ self.results[db] = "No hit found"
+
+ @staticmethod
+ def reversecomplement(seq):
+ # Make reverse complement strand
+ trans = str.maketrans("ATGC", "TACG")
+ return seq.translate(trans)[::-1]
+
+ @staticmethod
+ def compare_results(save, best_hsp, tmp_results, tmp_gene_split,
+ allowed_overlap):
+ """
+ Function for comparing hits and saving only the best hit
+ """
+
+ # Get data for comparison
+ hit_id = best_hsp['hit_id']
+ new_start_query = best_hsp['query_start']
+ new_end_query = best_hsp['query_end']
+ new_start_sbjct = int(best_hsp['sbjct_start'])
+ new_end_sbjct = int(best_hsp['sbjct_end'])
+ new_score = best_hsp['cal_score']
+ new_db_hit = best_hsp['sbjct_header']
+ new_contig = best_hsp['contig_name']
+ new_HSP = best_hsp['HSP_length']
+
+ # See if the best HSP fragment overlap with another allignment
+ # and keep the allignment with the highest score - if the new
+ # fragment is not providing new sequence
+ keys = list(tmp_results.keys())
+ for hit in keys:
+ hit_data = tmp_results[hit]
+ old_start_query = hit_data['query_start']
+ old_end_query = hit_data['query_end']
+ old_start_sbjct = int(hit_data['sbjct_start'])
+ old_end_sbjct = int(hit_data['sbjct_end'])
+ old_score = hit_data['cal_score']
+ old_db_hit = hit_data['sbjct_header']
+ old_contig = hit_data['contig_name']
+ old_HSP = hit_data['HSP_length']
+
+ remove_old = 0
+
+ # If they align to the same gene in the database they are
+ # compared
+ if new_db_hit == old_db_hit:
+
+ # If the hit provids additional sequence it is kept and
+ # the new coverage is saved otherwise the one with the
+ # highest score is kept
+ if(new_start_sbjct < old_start_sbjct
+ or new_end_sbjct > old_end_sbjct):
+
+ # Save the hits as split
+ tmp_gene_split[old_db_hit][hit_id] = 1
+ if hit not in tmp_gene_split[old_db_hit]:
+ tmp_gene_split[old_db_hit][hit] = 1
+
+ else:
+ if new_score > old_score:
+ # Set to remove old hit
+ remove_old = 1
+
+ # Save a split if the new hit still creats one
+ if(new_db_hit in tmp_gene_split
+ and hit_id not in tmp_gene_split[new_db_hit]):
+
+ tmp_gene_split[new_db_hit][hit_id] = 1
+ else:
+ save = 0
+
+ # If the old and new hit is not identical the
+ # possible saved gene split for the new hit is
+ # removed
+ if hit_id != hit:
+ if(new_db_hit in tmp_gene_split
+ and hit_id in tmp_gene_split[new_db_hit]):
+
+ del tmp_gene_split[new_db_hit][hit_id]
break
- # Getting the whole contig to extract extra query seqeunce
- contig = ''
- for seq_record in SeqIO.parse(inputfile, "fasta"):
- if seq_record.description.replace(" ", "") == hit['contig_name'].replace(" ", ""):
- contig = str(seq_record.seq)
- break
+ # If the hits comes form the same part of the contig
+ # sequnce but match different genes only the best hit is
+ # kept
+ if new_contig == old_contig:
+ print("Same contig: {} == {}".format(new_contig, old_contig))
+ print("\t{} vs {}".format(new_db_hit, old_db_hit))
+ # Check if saved hits overlaps with current hit
+ hit_union_length = (max(old_end_query, new_end_query)
+ - min(old_start_query, new_start_query))
+ hit_lengths_sum = ((old_end_query - old_start_query)
+ + (new_end_query - new_start_query))
+ overlap_len = (hit_lengths_sum - hit_union_length)
- # Extract extra sequence from query
- query_seq, homo_seq = self.get_query_align(hit, contig)
+ if overlap_len < allowed_overlap:
+ print("\tignore overlap ({}): {}".format(overlap_len, new_db_hit))
+ continue
+ print("\toverlap found ({}): {}".format(overlap_len, new_db_hit))
- # Saving the new alignment sequences
- self.gene_align_query[db][hit_id] = query_seq
- self.gene_align_homo[db][hit_id] = homo_seq
+ # If the two hits cover the exact same place on the
+ # contig only the percentage of identity is compared
+ if(old_start_query == new_start_query
+ and old_end_query == new_end_query):
- else:
- del gene_results[hit_id]
- if hit['sbjct_header'] in gene_split:
- del gene_split[hit['sbjct_header']]
-
- # Save the database result
- if gene_results:
- self.results[db] = gene_results
- else:
- self.results[db] = "No hit found"
-
- @staticmethod
- def reversecomplement(seq):
- # Make reverse complement strand
- trans = str.maketrans("ATGC", "TACG")
- return seq.translate(trans)[::-1]
-
- @staticmethod
- def compare_results(save, best_hsp, tmp_results, tmp_gene_split,
- allowed_overlap):
- """
- Function for comparing hits and saving only the best hit
- """
-
- # Get data for comparison
- hit_id = best_hsp['hit_id']
- new_start_query = best_hsp['query_start']
- new_end_query = best_hsp['query_end']
- new_start_sbjct = int(best_hsp['sbjct_start'])
- new_end_sbjct = int(best_hsp['sbjct_end'])
- new_score = best_hsp['cal_score']
- new_db_hit = best_hsp['sbjct_header']
- new_contig = best_hsp['contig_name']
- new_HSP = best_hsp['HSP_length']
-
- # See if the best HSP fragment overlap with another allignment
- # and keep the allignment with the highest score - if the new
- # fragment is not providing new sequence
- keys = list(tmp_results.keys())
- for hit in keys:
- hit_data = tmp_results[hit]
- old_start_query = hit_data['query_start']
- old_end_query = hit_data['query_end']
- old_start_sbjct = int(hit_data['sbjct_start'])
- old_end_sbjct = int(hit_data['sbjct_end'])
- old_score = hit_data['cal_score']
- old_db_hit = hit_data['sbjct_header']
- old_contig = hit_data['contig_name']
- old_HSP = hit_data['HSP_length']
-
- remove_old = 0
-
- # If they align to the same gene in the database they are
- # compared
- if new_db_hit == old_db_hit:
-
- # If the hit provids additional sequence it is kept and
- # the new coverage is saved otherwise the one with the
- # highest score is kept
- if(new_start_sbjct < old_start_sbjct
- or new_end_sbjct > old_end_sbjct):
-
- # Save the hits as split
- tmp_gene_split[old_db_hit][hit_id] = 1
- if hit not in tmp_gene_split[old_db_hit]:
- tmp_gene_split[old_db_hit][hit] = 1
+ if best_hsp['perc_ident'] > hit_data['perc_ident']:
- else:
- if new_score > old_score:
- # Set to remove old hit
- remove_old = 1
-
- # Save a split if the new hit still creats one
- if(new_db_hit in tmp_gene_split
- and hit_id not in tmp_gene_split[new_db_hit]):
-
- tmp_gene_split[new_db_hit][hit_id] = 1
- else:
- save = 0
-
- # If the old and new hit is not identical the
- # possible saved gene split for the new hit is
- # removed
- if hit_id != hit:
- if(new_db_hit in tmp_gene_split
- and hit_id in tmp_gene_split[new_db_hit]):
-
- del tmp_gene_split[new_db_hit][hit_id]
- break
-
- # If the hits comes form the same part of the contig
- # sequnce but match different genes only the best hit is
- # kept
- if new_contig == old_contig:
- # Check if saved hits overlaps with current hit
- overlap_len = (((old_end_query - old_start_query)
- + (new_end_query - new_start_query))
- - (max(old_end_query, new_end_query)
- - min(old_start_query, new_start_query)))
-
- if overlap_len < allowed_overlap:
- continue
+ # Set to remove old hit
+ remove_old = 1
- # If the two hits cover the exact same place on the
- # contig only the percentage of identity is compared
- if(old_start_query == new_start_query
- and old_end_query == new_end_query):
+ # Save a split if the new hit still creats one
+ if(new_db_hit in tmp_gene_split
+ and hit_id not in tmp_gene_split[new_db_hit]):
- if best_hsp['perc_ident'] > hit_data['perc_ident']:
-
- # Set to remove old hit
- remove_old = 1
-
- # Save a split if the new hit still creats one
- if(new_db_hit in tmp_gene_split
- and hit_id not in tmp_gene_split[new_db_hit]):
-
- tmp_gene_split[new_db_hit][hit_id] = 1
-
- elif best_hsp['perc_ident'] == hit_data['perc_ident']:
- # Save both
-
- # Save a split if the new hit still creats one
- if(new_db_hit in tmp_gene_split
- and hit_id not in tmp_gene_split[new_db_hit]):
-
- tmp_gene_split[new_db_hit][hit_id] = 1
- else:
- save = 0
-
- # Remove new gene from gene split if present
- if(new_db_hit in tmp_gene_split
- and hit_id in tmp_gene_split[new_db_hit]):
-
- del tmp_gene_split[new_db_hit][hit_id]
- break
-
- # If new hit overlaps with the saved hit
- elif((max(old_end_query, new_end_query)
- - min(old_start_query, new_start_query))
- <= ((old_end_query - old_start_query)
- + (new_end_query - new_start_query))):
-
- if new_score > old_score:
- # Set to remove old gene
- remove_old = 1
- # Save a split if the new hit still creats one
- if(new_db_hit in tmp_gene_split
- and hit_id not in tmp_gene_split[new_db_hit]):
-
- tmp_gene_split[new_db_hit][hit_id] = 1
-
- elif new_score == old_score:
- # If both genes are of same coverage
- # and identity is the same implyed by new_score == old_score
- # hit is chosen based on length
- if((int(best_hsp['perc_coverage']) ==
- int(hit_data['perc_coverage']))
- and new_HSP > old_HSP):
-
- # Set to remove old gene
- remove_old = 1
- elif((int(best_hsp['perc_coverage']) ==
- int(hit_data['perc_coverage']))
- and old_HSP > new_HSP):
-
- # Remove current hit
- save = 0
- elif((int(best_hsp['perc_coverage']) ==
- int(hit_data['perc_coverage']))
- and old_HSP==new_HSP):
- #Both hits has same coverage, and same identity
- # and same length, how to choose only one hit?
- pass
-
- # TODO
- # If new_score == old_score but identity and coverages are not the same.
- # which gene should be chosen?? Now they are both keept.
-
- # Save a split if the new hit creats one - both
- # hits are saved
- if(new_db_hit in tmp_gene_split
- and hit_id not in tmp_gene_split[new_db_hit]):
-
- tmp_gene_split[new_db_hit][hit_id] = 1
- else:
- # Remove new gene from gene split if present
- if(new_db_hit in tmp_gene_split
- and hit_id in tmp_gene_split[new_db_hit]):
-
- del tmp_gene_split[new_db_hit][hit_id]
-
- save = 0
- break
-
- # Remove old hit if new hit is better
- if remove_old == 1:
- del tmp_results[hit]
-
- # Remove gene from gene split if present
- if(old_db_hit in tmp_gene_split
- and hit in tmp_gene_split[old_db_hit]):
-
- del tmp_gene_split[old_db_hit][hit]
-
- return save, tmp_gene_split, tmp_results
-
- @staticmethod
- def calculate_new_length(gene_split, gene_results, hit):
- """
- Function for calcualting new length if the gene is split on
- several contigs
- """
- # Looping over splitted hits and calculate new length
- first = 1
- for split in gene_split[hit['sbjct_header']]:
-
- new_start = int(gene_results[split]['sbjct_start'])
- new_end = int(gene_results[split]['sbjct_end'])
-
- # Get the first HSP
- if first == 1:
- new_length = int(gene_results[split]['HSP_length'])
- old_start = new_start
- old_end = new_end
- first = 0
- continue
- if new_start < old_start:
- new_length = new_length + (old_start - new_start)
- old_start = new_start
-
- if new_end > old_end:
- new_length = new_length + (new_end - old_end)
- old_end = new_end
-
- return(new_length)
-
- @staticmethod
- def get_query_align(hit, contig):
- """
- Function for extracting extra seqeunce data to the query
- alignment if the full reference length are not covered
- """
-
- # Getting data needed to extract sequences
- query_seq = hit['query_string']
- homo_seq = hit['homo_string']
- sbjct_start = int(hit['sbjct_start'])
- sbjct_end = int(hit['sbjct_end'])
- query_start = int(hit['query_start'])
- query_end = int(hit['query_end'])
- length = int(hit['sbjct_length'])
-
- # If the alignment doesn't start at the first position data is
- # added to the begnning
- if sbjct_start != 1:
- missing = sbjct_start - 1
-
- if(query_start >= missing and hit['strand'] != 1
- or hit['strand'] == 1 and missing <= (len(contig) - query_end)):
-
- # Getting the query sequence.
- # If the the hit is on the other strand the characters
- # are reversed.
- if hit['strand'] == 1:
- start_pos = query_end
- end_pos = query_end + missing
- chars = contig[start_pos:end_pos]
- chars = Blaster.reversecomplement(chars)
- else:
- start_pos = query_start - missing - 1
- end_pos = query_start - 1
- chars = contig[start_pos:end_pos]
-
- query_seq = chars + str(query_seq)
- else:
- # Getting the query sequence.
- # If the the hit is on the other strand the characters
- # are reversed.
- if hit['strand'] == 1:
- if query_end == len(contig):
- query_seq = "-" * missing + str(query_seq)
- else:
- start_pos = query_end
- chars = contig[start_pos:]
- chars = Blaster.reversecomplement(chars)
-
- query_seq = ("-" * (missing - len(chars))
- + chars + str(query_seq))
- elif query_start < 3:
- query_seq = "-" * missing + str(query_seq)
- else:
- end_pos = query_start - 2
- chars = contig[0:end_pos]
-
- query_seq = ("-" * (missing - len(chars))
- + chars + str(query_seq))
-
- # Adding to the homo sequence
- spaces = " " * missing
- homo_seq = str(spaces) + str(homo_seq)
-
- # If the alignment dosen't end and the last position data is
- # added to the end
- if sbjct_end < length:
- missing = length - sbjct_end
-
- if(missing <= (len(contig) - query_end) and hit['strand'] != 1
- or hit['strand'] == 1 and query_start >= missing):
-
- # Getting the query sequence.
- # If the the hit is on the other strand the characters
- # are reversed.
- if hit['strand'] == 1:
- start_pos = query_start - missing - 1
- end_pos = query_start - 1
- chars = contig[start_pos:end_pos]
- chars = Blaster.reversecomplement(chars)
- else:
- start_pos = query_end
- end_pos = query_end + missing
- chars = contig[start_pos:end_pos]
-
- query_seq = query_seq + chars
- else:
- # If the hit is on the other strand the characters are reversed
- if hit['strand'] == 1:
- if query_start < 3:
- query_seq = query_seq + "-" * missing
- else:
- end_pos = query_start - 2
- chars = contig[0:end_pos]
- chars = Blaster.reversecomplement(chars)
-
- query_seq = (query_seq
- + chars + "-" * (missing - len(chars)))
- elif query_end == len(contig):
- query_seq = query_seq + "-" * missing
- else:
- start_pos = query_end
- chars = contig[start_pos:]
+ tmp_gene_split[new_db_hit][hit_id] = 1
- query_seq = query_seq + chars + "-" * (missing - len(chars))
+ elif best_hsp['perc_ident'] == hit_data['perc_ident']:
+ # Save both
- # Adding to the homo sequence
- spaces = " " * int(missing)
- homo_seq = str(homo_seq) + str(spaces)
+ # Save a split if the new hit still creats one
+ if(new_db_hit in tmp_gene_split
+ and hit_id not in tmp_gene_split[new_db_hit]):
- return query_seq, homo_seq
+ tmp_gene_split[new_db_hit][hit_id] = 1
+ else:
+ save = 0
+
+ # Remove new gene from gene split if present
+ if(new_db_hit in tmp_gene_split
+ and hit_id in tmp_gene_split[new_db_hit]):
+
+ del tmp_gene_split[new_db_hit][hit_id]
+ break
+
+ # If new hit overlaps with the saved hit
+ elif(hit_union_length <= hit_lengths_sum):
+ print("\t{} <= {}".format(hit_union_length, hit_lengths_sum))
+ print("\t\tScores: {} cmp {}".format(new_score, old_score))
+
+ if new_score > old_score:
+ # Set to remove old gene
+ remove_old = 1
+ # Save a split if the new hit still creats one
+ if(new_db_hit in tmp_gene_split
+ and hit_id not in tmp_gene_split[new_db_hit]):
+
+ tmp_gene_split[new_db_hit][hit_id] = 1
+
+ elif new_score == old_score:
+ # If both genes are of same coverage
+ # and identity is the same implied by new_score == old_score
+ # hit is chosen based on length
+ if((int(best_hsp['perc_coverage']) ==
+ int(hit_data['perc_coverage']))
+ and new_HSP > old_HSP):
+
+ # Set to remove old gene
+ remove_old = 1
+ elif((int(best_hsp['perc_coverage']) ==
+ int(hit_data['perc_coverage']))
+ and old_HSP > new_HSP):
+
+ # Remove current hit
+ save = 0
+ elif((int(best_hsp['perc_coverage']) ==
+ int(hit_data['perc_coverage']))
+ and old_HSP==new_HSP):
+ # Both hits has same coverage, and same identity
+ # and same length, how to choose only one hit?
+ pass
+
+ # TODO
+ # If new_score == old_score but identity and coverages are not the same.
+ # which gene should be chosen?? Now they are both keept.
+
+ # Save a split if the new hit creats one - both
+ # hits are saved
+ if(new_db_hit in tmp_gene_split
+ and hit_id not in tmp_gene_split[new_db_hit]):
+
+ tmp_gene_split[new_db_hit][hit_id] = 1
+ else:
+ # Remove new gene from gene split if present
+ if(new_db_hit in tmp_gene_split
+ and hit_id in tmp_gene_split[new_db_hit]):
+
+ del tmp_gene_split[new_db_hit][hit_id]
+
+ save = 0
+ break
+
+ # Remove old hit if new hit is better
+ if remove_old == 1:
+ del tmp_results[hit]
+
+ # Remove gene from gene split if present
+ if(old_db_hit in tmp_gene_split
+ and hit in tmp_gene_split[old_db_hit]):
+
+ del tmp_gene_split[old_db_hit][hit]
+
+ return save, tmp_gene_split, tmp_results
+
+ @staticmethod
+ def calculate_new_length(gene_split, gene_results, hit):
+ """
+ Function for calcualting new length if the gene is split on
+ several contigs
+ """
+ # Looping over splitted hits and calculate new length
+ first = 1
+ for split in gene_split[hit['sbjct_header']]:
+
+ new_start = int(gene_results[split]['sbjct_start'])
+ new_end = int(gene_results[split]['sbjct_end'])
+
+ # Get the first HSP
+ if first == 1:
+ new_length = int(gene_results[split]['HSP_length'])
+ old_start = new_start
+ old_end = new_end
+ first = 0
+ continue
+ if new_start < old_start:
+ new_length = new_length + (old_start - new_start)
+ old_start = new_start
+
+ if new_end > old_end:
+ new_length = new_length + (new_end - old_end)
+ old_end = new_end
+
+ return(new_length)
+
+ @staticmethod
+ def get_query_align(hit, contig):
+ """
+ Function for extracting extra seqeunce data to the query
+ alignment if the full reference length are not covered
+ """
+
+ # Getting data needed to extract sequences
+ query_seq = hit['query_string']
+ homo_seq = hit['homo_string']
+ sbjct_start = int(hit['sbjct_start'])
+ sbjct_end = int(hit['sbjct_end'])
+ query_start = int(hit['query_start'])
+ query_end = int(hit['query_end'])
+ length = int(hit['sbjct_length'])
+
+ # If the alignment doesn't start at the first position data is
+ # added to the begnning
+ if sbjct_start != 1:
+ missing = sbjct_start - 1
+
+ if(query_start >= missing and hit['strand'] != 1
+ or hit['strand'] == 1 and missing <= (len(contig) - query_end)):
+
+ # Getting the query sequence.
+ # If the the hit is on the other strand the characters
+ # are reversed.
+ if hit['strand'] == 1:
+ start_pos = query_end
+ end_pos = query_end + missing
+ chars = contig[start_pos:end_pos]
+ chars = Blaster.reversecomplement(chars)
+ else:
+ start_pos = query_start - missing - 1
+ end_pos = query_start - 1
+ chars = contig[start_pos:end_pos]
+
+ query_seq = chars + str(query_seq)
+ else:
+ # Getting the query sequence.
+ # If the the hit is on the other strand the characters
+ # are reversed.
+ if hit['strand'] == 1:
+ if query_end == len(contig):
+ query_seq = "-" * missing + str(query_seq)
+ else:
+ start_pos = query_end
+ chars = contig[start_pos:]
+ chars = Blaster.reversecomplement(chars)
+
+ query_seq = ("-" * (missing - len(chars))
+ + chars + str(query_seq))
+ elif query_start < 3:
+ query_seq = "-" * missing + str(query_seq)
+ else:
+ end_pos = query_start - 2
+ chars = contig[0:end_pos]
+
+ query_seq = ("-" * (missing - len(chars))
+ + chars + str(query_seq))
+
+ # Adding to the homo sequence
+ spaces = " " * missing
+ homo_seq = str(spaces) + str(homo_seq)
+
+ # If the alignment dosen't end and the last position data is
+ # added to the end
+ if sbjct_end < length:
+ missing = length - sbjct_end
+
+ if(missing <= (len(contig) - query_end) and hit['strand'] != 1
+ or hit['strand'] == 1 and query_start >= missing):
+
+ # Getting the query sequence.
+ # If the the hit is on the other strand the characters
+ # are reversed.
+ if hit['strand'] == 1:
+ start_pos = query_start - missing - 1
+ end_pos = query_start - 1
+ chars = contig[start_pos:end_pos]
+ chars = Blaster.reversecomplement(chars)
+ else:
+ start_pos = query_end
+ end_pos = query_end + missing
+ chars = contig[start_pos:end_pos]
+
+ query_seq = query_seq + chars
+ else:
+ # If the hit is on the other strand the characters are reversed
+ if hit['strand'] == 1:
+ if query_start < 3:
+ query_seq = query_seq + "-" * missing
+ else:
+ end_pos = query_start - 2
+ chars = contig[0:end_pos]
+ chars = Blaster.reversecomplement(chars)
+
+ query_seq = (query_seq
+ + chars + "-" * (missing - len(chars)))
+ elif query_end == len(contig):
+ query_seq = query_seq + "-" * missing
+ else:
+ start_pos = query_end
+ chars = contig[start_pos:]
+
+ query_seq = query_seq + chars + "-" * (missing - len(chars))
+
+ # Adding to the homo sequence
+ spaces = " " * int(missing)
+ homo_seq = str(homo_seq) + str(spaces)
+
+ return query_seq, homo_seq
=====================================
cgecore/blaster/run_blaster.py
=====================================
@@ -36,6 +36,11 @@ parser.add_argument("-t", "--threshold", dest="threshold",
help="Blast threshold for identity",
type=float,
default=0.90)
+parser.add_argument("--overlap",
+ help=("Allow hits/genes to overlap with this number of "
+ "nucleotides. Default: 30."),
+ type=int,
+ default=30)
args = parser.parse_args()
@@ -94,7 +99,8 @@ for inp_file in input_list:
# Calling blast and parsing output
blast_run = Blaster(inputfile=inp_file, databases=databases,
db_path=db_path, out_path=out_path, min_cov=min_cov,
- threshold=threshold, blast=blast)
+ threshold=threshold, blast=blast,
+ allowed_overlap=args.overlap)
results = blast_run.results
query_align = blast_run.gene_align_query
=====================================
cgecore/cgefinder.py
=====================================
@@ -32,9 +32,9 @@ class CGEFinder():
def kma(inputfile_1, out_path, databases, db_path_kma, min_cov=0.6,
threshold=0.9, kma_path="cge/kma/kma", sample_name="",
inputfile_2=None, kma_mrs=None, kma_gapopen=None,
- kma_gapextend=None, kma_penalty=None, kma_reward=None, kma_pm=None,
- kma_fpm=None, kma_memmode=False, kma_nanopore=False, debug=False,
- kma_add_args=None):
+ kma_gapextend=None, kma_penalty=None, kma_reward=None,
+ kma_apm=None, kma_memmode=False, kma_nanopore=False, debug=False,
+ kma_add_args=None, kma_cge=False, kma_1t1=False):
"""
I expect that there will only be one hit pr gene, but if there are
more, I assume that the sequence of the hits are the same in the res
@@ -73,10 +73,12 @@ class CGEFinder():
kma_cmd += " -penalty " + str(kma_penalty)
if(kma_reward is not None):
kma_cmd += " -reward " + str(kma_reward)
- if(kma_pm is not None):
- kma_cmd += " -pm " + kma_pm
- if(kma_fpm is not None):
- kma_cmd += " -fpm " + kma_fpm
+ if(kma_apm is not None):
+ kma_cmd += " -apm " + kma_apm
+ if(kma_cge):
+ kma_cmd += " -cge "
+ if(kma_1t1):
+ kma_cmd += " -1t1 "
if (kma_memmode):
kma_cmd += " -mem_mode "
if (kma_nanopore):
=====================================
cgecore/cmdline.py
=====================================
@@ -258,7 +258,7 @@ class Program:
'else echo "Done" >> {stderr}; fi\n').format(
stderr=self.stderr))
f.write('exit $ec\n')
- os.chmod(script, 0o744)
+ os.chmod(script, 0o755)
if self.queue is not None:
# Setup execution of shell script through TORQUE
=====================================
debian/changelog
=====================================
@@ -1,3 +1,13 @@
+python-cgecore (1.5.0-1) unstable; urgency=medium
+
+ * New upstream version
+ * debhelper-compat 12
+ * Standards-Version: 4.4.0
+ * Testsuite: autopkgtest-pkg-python
+ * Set upstream metadata fields: Name.
+
+ -- Andreas Tille <tille at debian.org> Thu, 01 Aug 2019 22:24:02 +0200
+
python-cgecore (1.4.4-1) unstable; urgency=medium
* Initial release (Closes: #929496)
=====================================
debian/compat deleted
=====================================
@@ -1 +0,0 @@
-12
=====================================
debian/control
=====================================
@@ -2,13 +2,14 @@ Source: python-cgecore
Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
Uploaders: Andreas Tille <tille at debian.org>
Section: python
+Testsuite: autopkgtest-pkg-python
Priority: optional
-Build-Depends: debhelper (>= 12~),
+Build-Depends: debhelper-compat (= 12),
dh-python,
python3,
python3-setuptools,
python3-biopython
-Standards-Version: 4.3.0
+Standards-Version: 4.4.0
Vcs-Browser: https://salsa.debian.org/med-team/python-cgecore
Vcs-Git: https://salsa.debian.org/med-team/python-cgecore.git
Homepage: https://bitbucket.org/genomicepidemiology/cge_core_module/
=====================================
debian/upstream/metadata
=====================================
@@ -0,0 +1 @@
+Name: cgecore
View it on GitLab: https://salsa.debian.org/med-team/python-cgecore/compare/de1163aef85089737bf65b311062e29ce35441dc...8cc4a6e93212e45fc7c927da394880e3a4fdcfe5
--
View it on GitLab: https://salsa.debian.org/med-team/python-cgecore/compare/de1163aef85089737bf65b311062e29ce35441dc...8cc4a6e93212e45fc7c927da394880e3a4fdcfe5
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20190801/0893bcf8/attachment-0001.html>
More information about the debian-med-commit
mailing list