Compare commits
28 Commits
Author | SHA1 | Date | |
---|---|---|---|
8ec13a294c | |||
3e45c34491 | |||
c2f3d90dc1 | |||
6b732d11d3 | |||
9eb833a0af | |||
6b7b0e3bd1 | |||
47691a8f58 | |||
b908b581c8 | |||
03c174fd7a | |||
2156588ff6 | |||
6ff29c6a6a | |||
51a3c68fb5 | |||
da91ffc2c7 | |||
c884615522 | |||
cb53381863 | |||
72b3e5d872 | |||
238e9f70f3 | |||
e099a16624 | |||
847c9c816d | |||
6026129ca8 | |||
169b6514b4 | |||
89b0c48141 | |||
7c02782e3c | |||
ecc4c2c78b | |||
f5413381fd | |||
3e93cfff7b | |||
6d445fe3ad | |||
824deb7e21 |
@ -1,3 +1,9 @@
|
||||
import codecs
|
||||
|
||||
def unescaped_str(arg_str):
|
||||
return arg_str.encode('latin-1', 'backslashreplace').decode('unicode-escape')
|
||||
|
||||
|
||||
def __addInputOption(optionManager):
|
||||
|
||||
optionManager.add_argument(
|
||||
@ -43,7 +49,13 @@ def __addImportInputOption(optionManager):
|
||||
action="store_const", dest="obi:inputformat",
|
||||
default=None,
|
||||
const=b'silva',
|
||||
help="Input file is in SILVA fasta format")
|
||||
help="Input file is in SILVA fasta format. If NCBI taxonomy provided with --taxonomy, taxid and scientific name will be added for each sequence.")
|
||||
|
||||
group.add_argument('--rdp-input',
|
||||
action="store_const", dest="obi:inputformat",
|
||||
default=None,
|
||||
const=b'rdp',
|
||||
help="Input file is in RDP training set fasta format. If NCBI taxonomy provided with --taxonomy, taxid and scientific name will be added for each sequence.")
|
||||
|
||||
group.add_argument('--embl-input',
|
||||
action="store_const", dest="obi:inputformat",
|
||||
@ -133,7 +145,7 @@ def __addTabularOption(optionManager):
|
||||
group.add_argument('--sep',
|
||||
action="store", dest="obi:sep",
|
||||
default="\t",
|
||||
type=str,
|
||||
type=unescaped_str,
|
||||
help="Column separator")
|
||||
|
||||
|
||||
|
@ -26,7 +26,7 @@ import sys
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Annotate sequences with their corresponding NCBI taxid found from the taxon scientific name."
|
||||
__title__="Annotate sequences with their corresponding NCBI taxid found from the taxon scientific name"
|
||||
|
||||
|
||||
|
||||
@ -45,7 +45,7 @@ def addOptions(parser):
|
||||
metavar="<TAXID_TAG>",
|
||||
default=b"TAXID",
|
||||
help="Name of the tag to store the found taxid "
|
||||
"(default: 'TAXID'.")
|
||||
"(default: 'TAXID').")
|
||||
|
||||
group.add_argument('-n', '--taxon-name-tag',
|
||||
action="store",
|
||||
@ -53,7 +53,7 @@ def addOptions(parser):
|
||||
metavar="<SCIENTIFIC_NAME_TAG>",
|
||||
default=b"SCIENTIFIC_NAME",
|
||||
help="Name of the tag giving the scientific name of the taxon "
|
||||
"(default: 'SCIENTIFIC_NAME'.")
|
||||
"(default: 'SCIENTIFIC_NAME').")
|
||||
|
||||
group.add_argument('-g', '--try-genus-match',
|
||||
action="store_true", dest="addtaxids:try_genus_match",
|
||||
@ -174,6 +174,7 @@ def run(config):
|
||||
taxid_column[i] = taxon.taxid
|
||||
found_count+=1
|
||||
elif try_genus: # try finding genus or other parent taxon from the first word
|
||||
#print(i, o_view[i].id)
|
||||
taxon_name_sp = taxon_name.split(b" ")
|
||||
taxon = taxo.get_taxon_by_name(taxon_name_sp[0], res_anc)
|
||||
if taxon is not None:
|
||||
|
@ -19,7 +19,7 @@ import time
|
||||
import sys
|
||||
|
||||
|
||||
__title__="Aligns one sequence column with itself or two sequence columns"
|
||||
__title__="Align one sequence column with itself or two sequence columns"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
@ -158,7 +158,7 @@ def run(config):
|
||||
i_view_name = i_uri.split(b"/")[0]
|
||||
i_column_name = b""
|
||||
i_element_name = b""
|
||||
if len(i_uri.split(b"/")) == 2:
|
||||
if len(i_uri.split(b"/")) >= 2:
|
||||
i_column_name = i_uri.split(b"/")[1]
|
||||
if len(i_uri.split(b"/")) == 3:
|
||||
i_element_name = i_uri.split(b"/")[2]
|
||||
@ -181,7 +181,7 @@ def run(config):
|
||||
i_dms_name_2 = i_dms_2.name
|
||||
i_uri_2 = input_2[1]
|
||||
original_i_view_name_2 = i_uri_2.split(b"/")[0]
|
||||
if len(i_uri_2.split(b"/")) == 2:
|
||||
if len(i_uri_2.split(b"/")) >= 2:
|
||||
i_column_name_2 = i_uri_2.split(b"/")[1]
|
||||
if len(i_uri_2.split(b"/")) == 3:
|
||||
i_element_name_2 = i_uri_2.split(b"/")[2]
|
||||
|
@ -23,7 +23,7 @@ import os
|
||||
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
__title__="Aligns paired-ended reads"
|
||||
__title__="Align paired-ended reads"
|
||||
|
||||
|
||||
|
||||
|
@ -16,7 +16,7 @@ import sys
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Tag a set of sequences for PCR and sequencing errors identification"
|
||||
__title__="Build a reference database for ecotag"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
@ -31,10 +31,9 @@ def addOptions(parser):
|
||||
group.add_argument('--threshold','-t',
|
||||
action="store", dest="build_ref_db:threshold",
|
||||
metavar='<THRESHOLD>',
|
||||
default=0.0,
|
||||
default=0.99,
|
||||
type=float,
|
||||
help="Score threshold as a normalized identity, e.g. 0.95 for an identity of 95%%. Default: 0.00"
|
||||
" (no threshold).")
|
||||
help="Score threshold as a normalized identity, e.g. 0.95 for an identity of 95%%. Default: 0.99.")
|
||||
|
||||
|
||||
def run(config):
|
||||
|
@ -22,7 +22,7 @@ import sys
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Concatenate views."
|
||||
__title__="Concatenate views"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
@ -97,7 +97,7 @@ def run(config):
|
||||
Column.new_column(o_view, REVERSE_QUALITY_COLUMN, OBI_QUAL, associated_column_name=REVERSE_SEQUENCE_COLUMN, associated_column_version=o_view[REVERSE_SEQUENCE_COLUMN].version)
|
||||
|
||||
# Initialize multiple elements columns
|
||||
if type(output_0)==BufferedWriter:
|
||||
if type(output_0)!=BufferedWriter:
|
||||
dict_cols = {}
|
||||
for v_uri in config["cat"]["views_to_cat"]:
|
||||
v = open_uri(v_uri)[1]
|
||||
|
@ -54,11 +54,11 @@ def addOptions(parser):
|
||||
default=False,
|
||||
help="Only sequences labeled as heads are kept in the output. Default: False")
|
||||
|
||||
group.add_argument('--cluster-tags', '-C',
|
||||
action="store_true",
|
||||
dest="clean:cluster-tags",
|
||||
default=False,
|
||||
help="Adds tags for each sequence giving its cluster's head and weight for each sample.")
|
||||
# group.add_argument('--cluster-tags', '-C',
|
||||
# action="store_true",
|
||||
# dest="clean:cluster-tags",
|
||||
# default=False,
|
||||
# help="Adds tags for each sequence giving its cluster's head and weight for each sample.")
|
||||
|
||||
group.add_argument('--thread-count','-p', # TODO should probably be in a specific option group
|
||||
action="store", dest="clean:thread-count",
|
||||
@ -142,4 +142,5 @@ def run(config):
|
||||
|
||||
i_dms.close(force=True)
|
||||
|
||||
logger("info", "Done.")
|
||||
logger("info", "Done.")
|
||||
|
@ -10,7 +10,7 @@ from obitools3.dms.capi.obiview cimport COUNT_COLUMN
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Counts sequence records"
|
||||
__title__="Count sequence records"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
@ -29,6 +29,12 @@ def addOptions(parser):
|
||||
default=False,
|
||||
help="Prints only the total count of sequence records (if a sequence has no `count` attribute, its default count is 1) (default: False).")
|
||||
|
||||
group.add_argument('-c','--count-tag',
|
||||
action="store", dest="count:countcol",
|
||||
default='COUNT',
|
||||
type=str,
|
||||
help="Name of the tag/column associated with the count information (default: COUNT).")
|
||||
|
||||
|
||||
def run(config):
|
||||
|
||||
@ -41,18 +47,20 @@ def run(config):
|
||||
if input is None:
|
||||
raise Exception("Could not read input")
|
||||
entries = input[1]
|
||||
|
||||
|
||||
countcol = config['count']['countcol']
|
||||
|
||||
count1 = len(entries)
|
||||
count2 = 0
|
||||
|
||||
if COUNT_COLUMN in entries and ((config['count']['sequence'] == config['count']['all']) or (config['count']['all'])) :
|
||||
if countcol in entries and ((config['count']['sequence'] == config['count']['all']) or (config['count']['all'])) :
|
||||
for e in entries:
|
||||
PyErr_CheckSignals()
|
||||
count2+=e[COUNT_COLUMN]
|
||||
count2+=e[countcol]
|
||||
|
||||
if COUNT_COLUMN in entries and (config['count']['sequence'] == config['count']['all']):
|
||||
if countcol in entries and (config['count']['sequence'] == config['count']['all']):
|
||||
print(count1,count2)
|
||||
elif COUNT_COLUMN in entries and config['count']['all']:
|
||||
elif countcol in entries and config['count']['all']:
|
||||
print(count2)
|
||||
else:
|
||||
print(count1)
|
||||
|
@ -258,6 +258,13 @@ def Filter_generator(options, tax_filter, i_view):
|
||||
|
||||
|
||||
def Taxonomy_filter_generator(taxo, options):
|
||||
|
||||
if (("required_ranks" in options and options["required_ranks"]) or \
|
||||
("required_taxids" in options and options["required_taxids"]) or \
|
||||
("ignored_taxids" in options and options["ignored_taxids"])) and \
|
||||
(taxo is None):
|
||||
raise RollbackException("obi grep error: can't use taxonomy options without providing a taxonomy. Rollbacking view")
|
||||
|
||||
if taxo is not None:
|
||||
def tax_filter(seq):
|
||||
good = True
|
||||
|
@ -16,7 +16,7 @@ from io import BufferedWriter
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Keep the N first lines of a view."
|
||||
__title__="Keep the N first lines of a view"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
|
@ -34,14 +34,17 @@ from obitools3.dms.capi.obidms cimport obi_import_view
|
||||
from obitools3.dms.capi.obitypes cimport obitype_t, \
|
||||
OBI_VOID, \
|
||||
OBI_QUAL, \
|
||||
OBI_STR
|
||||
OBI_STR, \
|
||||
OBI_INT
|
||||
|
||||
from obitools3.dms.capi.obierrno cimport obi_errno
|
||||
|
||||
from obitools3.apps.optiongroups import addImportInputOption, \
|
||||
addTabularInputOption, \
|
||||
addTaxdumpInputOption, \
|
||||
addMinimalOutputOption
|
||||
addMinimalOutputOption, \
|
||||
addNoProgressBarOption, \
|
||||
addTaxonomyOption
|
||||
|
||||
from obitools3.uri.decode import open_uri
|
||||
|
||||
@ -50,9 +53,10 @@ from obitools3.apps.config import logger
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
from io import BufferedWriter
|
||||
import ast
|
||||
|
||||
|
||||
__title__="Imports sequences from different formats into a DMS"
|
||||
__title__="Import sequences from different formats into a DMS"
|
||||
|
||||
|
||||
default_config = { 'destview' : None,
|
||||
@ -69,7 +73,9 @@ def addOptions(parser):
|
||||
addImportInputOption(parser)
|
||||
addTabularInputOption(parser)
|
||||
addTaxdumpInputOption(parser)
|
||||
addTaxonomyOption(parser)
|
||||
addMinimalOutputOption(parser)
|
||||
addNoProgressBarOption(parser)
|
||||
|
||||
group = parser.add_argument_group('obi import specific options')
|
||||
|
||||
@ -85,6 +91,10 @@ def addOptions(parser):
|
||||
help="If importing a view into another DMS, do it by importing each line, saving disk space if the original view "
|
||||
"has a line selection associated.")
|
||||
|
||||
# group.add_argument('--only-id',
|
||||
# action="store", dest="import:onlyid",
|
||||
# help="only id")
|
||||
|
||||
def run(config):
|
||||
|
||||
cdef tuple input
|
||||
@ -97,6 +107,7 @@ def run(config):
|
||||
cdef bint get_quality
|
||||
cdef bint NUC_SEQS_view
|
||||
cdef bint silva
|
||||
cdef bint rdp
|
||||
cdef int nb_elts
|
||||
cdef object d
|
||||
cdef View view
|
||||
@ -180,6 +191,16 @@ def run(config):
|
||||
logger("info", "Done.")
|
||||
return
|
||||
|
||||
# Open taxonomy if there is one
|
||||
if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
|
||||
taxo_uri = open_uri(config['obi']['taxoURI'])
|
||||
if taxo_uri is None or taxo_uri[2] == bytes:
|
||||
raise Exception("Couldn't open taxonomy")
|
||||
taxo = taxo_uri[1]
|
||||
|
||||
else :
|
||||
taxo = None
|
||||
|
||||
# If importing a view between two DMS and not wanting to save space if line selection in original view, use C API
|
||||
if isinstance(input[1], View) and not config['import']['space_priority']:
|
||||
if obi_import_view(input[0].name_with_full_path, o_dms.name_with_full_path, input[1].name, tobytes((config['obi']['outputURI'].split('/'))[-1])) < 0 :
|
||||
@ -192,8 +213,11 @@ def run(config):
|
||||
logger("info", "Done.")
|
||||
return
|
||||
|
||||
if entry_count >= 0:
|
||||
# Reinitialize the progress bar
|
||||
if entry_count >= 0 and config['obi']['noprogressbar'] == False:
|
||||
pb = ProgressBar(entry_count, config)
|
||||
else:
|
||||
pb = None
|
||||
|
||||
NUC_SEQS_view = False
|
||||
if isinstance(output[1], View) :
|
||||
@ -202,20 +226,25 @@ def run(config):
|
||||
NUC_SEQS_view = True
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
|
||||
# Save basic columns in variables for optimization
|
||||
if NUC_SEQS_view :
|
||||
id_col = view[ID_COLUMN]
|
||||
def_col = view[DEFINITION_COLUMN]
|
||||
seq_col = view[NUC_SEQUENCE_COLUMN]
|
||||
|
||||
# Prepare taxon scientific name if SILVA file
|
||||
if 'inputformat' in config['obi'] and config['obi']['inputformat'] == b"silva":
|
||||
# Prepare taxon scientific name and taxid refs if RDP or SILVA file
|
||||
silva = False
|
||||
rdp = False
|
||||
if 'inputformat' in config['obi'] and (config['obi']['inputformat'] == b"silva" or config['obi']['inputformat'] == b"rdp"):
|
||||
#if taxo is None:
|
||||
# raise Exception("A taxonomy (as built by 'obi import --taxdump') must be provided for SILVA and RDP files")
|
||||
silva = True
|
||||
sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
|
||||
else:
|
||||
silva = False
|
||||
|
||||
rdp = True
|
||||
if taxo is not None:
|
||||
sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
|
||||
taxid_col = Column.new_column(view, TAXID_COLUMN, OBI_INT)
|
||||
|
||||
dcols = {}
|
||||
|
||||
# First read through the entries to prepare columns with dictionaries as they are very time-expensive to rewrite
|
||||
@ -265,8 +294,13 @@ def run(config):
|
||||
# Reinitialize the input
|
||||
if isinstance(input[0], CompressedFile):
|
||||
input_is_file = True
|
||||
if entry_count >= 0:
|
||||
|
||||
# Reinitialize the progress bar
|
||||
if entry_count >= 0 and config['obi']['noprogressbar'] == False:
|
||||
pb = ProgressBar(entry_count, config)
|
||||
else:
|
||||
pb = None
|
||||
|
||||
try:
|
||||
input[0].close()
|
||||
except AttributeError:
|
||||
@ -275,6 +309,11 @@ def run(config):
|
||||
if input is None:
|
||||
raise Exception("Could not open input URI")
|
||||
|
||||
# if 'onlyid' in config['import']:
|
||||
# onlyid = tobytes(config['import']['onlyid'])
|
||||
# else:
|
||||
# onlyid = None
|
||||
|
||||
entries = input[1]
|
||||
i = 0
|
||||
for entry in entries :
|
||||
@ -292,6 +331,9 @@ def run(config):
|
||||
elif not i%50000:
|
||||
logger("info", "Imported %d entries", i)
|
||||
|
||||
# if onlyid is not None and entry.id != onlyid:
|
||||
# continue
|
||||
|
||||
try:
|
||||
|
||||
if NUC_SEQS_view:
|
||||
@ -307,11 +349,18 @@ def run(config):
|
||||
if get_quality:
|
||||
qual_col[i] = entry.quality
|
||||
|
||||
# Parse taxon scientific name if SILVA file
|
||||
if silva:
|
||||
sci_name = entry.definition.split(b";")[-1]
|
||||
sci_name_col[i] = sci_name
|
||||
|
||||
# Parse taxon scientific name if RDP file
|
||||
if (rdp or silva) and taxo is not None:
|
||||
sci_names = entry.definition.split(b";")
|
||||
for sci_name in reversed(sci_names):
|
||||
if sci_name.split()[0] != b'unidentified' and sci_name.split()[0] != b'uncultured' and sci_name.split()[0] != b'metagenome' :
|
||||
taxon = taxo.get_taxon_by_name(sci_name)
|
||||
if taxon is not None:
|
||||
sci_name_col[i] = taxon.name
|
||||
taxid_col[i] = taxon.taxid
|
||||
#print(taxid_col[i], sci_name_col[i])
|
||||
break
|
||||
|
||||
for tag in entry :
|
||||
|
||||
if tag != ID_COLUMN and tag != DEFINITION_COLUMN and tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN : # TODO dirty
|
||||
@ -323,27 +372,39 @@ def run(config):
|
||||
tag = COUNT_COLUMN
|
||||
if tag[:7] == b"merged_":
|
||||
tag = MERGED_PREFIX+tag[7:]
|
||||
|
||||
|
||||
if type(value) == bytes and value[:1]==b"[" :
|
||||
try:
|
||||
if type(eval(value)) == list:
|
||||
value = eval(value)
|
||||
#print(value)
|
||||
except:
|
||||
pass
|
||||
|
||||
if tag not in dcols :
|
||||
|
||||
|
||||
value_type = type(value)
|
||||
nb_elts = 1
|
||||
value_obitype = OBI_VOID
|
||||
dict_col = False
|
||||
|
||||
if value_type == dict or value_type == list :
|
||||
|
||||
if value_type == dict :
|
||||
nb_elts = len(value)
|
||||
elt_names = list(value)
|
||||
if value_type == dict :
|
||||
dict_col = True
|
||||
dict_col = True
|
||||
else :
|
||||
nb_elts = 1
|
||||
elt_names = None
|
||||
|
||||
|
||||
if value_type == list :
|
||||
tuples = True
|
||||
else:
|
||||
tuples = False
|
||||
|
||||
value_obitype = get_obitype(value)
|
||||
|
||||
|
||||
if value_obitype != OBI_VOID :
|
||||
dcols[tag] = (Column.new_column(view, tag, value_obitype, nb_elements_per_line=nb_elts, elements_names=elt_names, dict_column=dict_col), value_obitype)
|
||||
dcols[tag] = (Column.new_column(view, tag, value_obitype, nb_elements_per_line=nb_elts, elements_names=elt_names, dict_column=dict_col, tuples=tuples), value_obitype)
|
||||
|
||||
# Fill value
|
||||
dcols[tag][0][i] = value
|
||||
@ -371,8 +432,8 @@ def run(config):
|
||||
# Fill value
|
||||
dcols[tag][0][i] = value
|
||||
|
||||
except IndexError :
|
||||
|
||||
except (IndexError, OverflowError):
|
||||
|
||||
value_type = type(value)
|
||||
old_column = dcols[tag][0]
|
||||
old_nb_elements_per_line = old_column.nb_elements_per_line
|
||||
@ -419,18 +480,19 @@ def run(config):
|
||||
dcols[tag][0][i] = value
|
||||
|
||||
except Exception as e:
|
||||
print("\nCould not import sequence id:", entry.id, "(error raised:", e, ")")
|
||||
print("\nCould not import sequence:", repr(entry), "(error raised:", e, ")")
|
||||
if 'skiperror' in config['obi'] and not config['obi']['skiperror']:
|
||||
raise e
|
||||
else:
|
||||
pass
|
||||
i-=1 # overwrite problematic entry
|
||||
|
||||
i+=1
|
||||
|
||||
if pb is not None:
|
||||
pb(i, force=True)
|
||||
print("", file=sys.stderr)
|
||||
logger("info", "Imported %d entries", i)
|
||||
logger("info", "Imported %d entries", len(view))
|
||||
|
||||
# Save command config in View and DMS comments
|
||||
command_line = " ".join(sys.argv[1:])
|
||||
|
@ -24,7 +24,7 @@ from cpython.exc cimport PyErr_CheckSignals
|
||||
from io import BufferedWriter
|
||||
|
||||
|
||||
__title__="Assigns sequence records to the corresponding experiment/sample based on DNA tags and primers"
|
||||
__title__="Assign sequence records to the corresponding experiment/sample based on DNA tags and primers"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
@ -450,7 +450,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
|
||||
sequences[1] = sequences[1][reversematch[1][2]:]
|
||||
if not directmatch[0].forward:
|
||||
sequences[1] = sequences[1].reverse_complement
|
||||
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
|
||||
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
|
||||
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
|
||||
else:
|
||||
sequences[0] = sequences[0][reversematch[1][2]:]
|
||||
|
44
python/obitools3/commands/rm.pyx
Normal file
44
python/obitools3/commands/rm.pyx
Normal file
@ -0,0 +1,44 @@
|
||||
#cython: language_level=3
|
||||
|
||||
from obitools3.uri.decode import open_uri
|
||||
from obitools3.apps.config import logger
|
||||
from obitools3.dms import DMS
|
||||
from obitools3.apps.optiongroups import addMinimalInputOption
|
||||
from obitools3.dms.view.view cimport View
|
||||
import os
|
||||
|
||||
__title__="Delete a view"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
addMinimalInputOption(parser)
|
||||
|
||||
def run(config):
|
||||
|
||||
DMS.obi_atexit()
|
||||
|
||||
logger("info", "obi rm")
|
||||
|
||||
# Open the input
|
||||
input = open_uri(config['obi']['inputURI'])
|
||||
if input is None:
|
||||
raise Exception("Could not read input")
|
||||
|
||||
# Check that it's a view
|
||||
if isinstance(input[1], View) :
|
||||
view = input[1]
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# Get the path to the view file to remove
|
||||
path = input[0].full_path # dms path
|
||||
path+=b"/VIEWS/"
|
||||
path+=view.name
|
||||
path+=b".obiview"
|
||||
|
||||
# Close the view and the DMS
|
||||
view.close()
|
||||
input[0].close(force=True)
|
||||
|
||||
# Rm
|
||||
os.remove(path)
|
@ -36,7 +36,7 @@ NULL_VALUE = {OBI_BOOL: OBIBool_NA,
|
||||
OBI_STR: b""}
|
||||
|
||||
|
||||
__title__="Sort view lines according to the value of a given attribute."
|
||||
__title__="Sort view lines according to the value of a given attribute"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
|
@ -16,7 +16,7 @@ import sys
|
||||
from cpython.exc cimport PyErr_CheckSignals
|
||||
|
||||
|
||||
__title__="Compute basic statistics for attribute values."
|
||||
__title__="Compute basic statistics for attribute values"
|
||||
|
||||
'''
|
||||
`obi stats` computes basic statistics for attribute values of sequence records.
|
||||
@ -154,7 +154,7 @@ def run(config):
|
||||
else :
|
||||
taxo = None
|
||||
|
||||
statistics = set(config['stats']['minimum']) | set(config['stats']['maximum']) | set(config['stats']['mean'])
|
||||
statistics = set(config['stats']['minimum']) | set(config['stats']['maximum']) | set(config['stats']['mean']) | set(config['stats']['var']) | set(config['stats']['sd'])
|
||||
total = 0
|
||||
catcount={}
|
||||
totcount={}
|
||||
@ -195,7 +195,7 @@ def run(config):
|
||||
except KeyError:
|
||||
totcount[category]=totcount.get(category,0)+1
|
||||
for var in statistics:
|
||||
if var in line:
|
||||
if var in line and line[var] is not None:
|
||||
v = line[var]
|
||||
if var not in values:
|
||||
values[var]={}
|
||||
@ -238,14 +238,34 @@ def run(config):
|
||||
else:
|
||||
sdvar= "%s"
|
||||
|
||||
hcat = "\t".join([pcat % x for x in config['stats']['categories']]) + "\t" +\
|
||||
"\t".join([minvar % x for x in config['stats']['minimum']]) + "\t" +\
|
||||
"\t".join([maxvar % x for x in config['stats']['maximum']]) + "\t" +\
|
||||
"\t".join([meanvar % x for x in config['stats']['mean']]) + "\t" +\
|
||||
"\t".join([varvar % x for x in config['stats']['var']]) + "\t" +\
|
||||
"\t".join([sdvar % x for x in config['stats']['sd']]) + \
|
||||
"\t count" + \
|
||||
"\t total"
|
||||
hcat = ""
|
||||
|
||||
for x in config['stats']['categories']:
|
||||
hcat += pcat % x
|
||||
hcat += "\t"
|
||||
|
||||
for x in config['stats']['minimum']:
|
||||
hcat += minvar % x
|
||||
hcat += "\t"
|
||||
|
||||
for x in config['stats']['maximum']:
|
||||
hcat += maxvar % x
|
||||
hcat += "\t"
|
||||
|
||||
for x in config['stats']['mean']:
|
||||
hcat += meanvar % x
|
||||
hcat += "\t"
|
||||
|
||||
for x in config['stats']['var']:
|
||||
hcat += varvar % x
|
||||
hcat += "\t"
|
||||
|
||||
for x in config['stats']['sd']:
|
||||
hcat += sdvar % x
|
||||
hcat += "\t"
|
||||
|
||||
hcat += "count\ttotal"
|
||||
|
||||
print(hcat)
|
||||
sorted_stats = sorted(catcount.items(), key = lambda kv:(totcount[kv[0]]), reverse=True)
|
||||
for i in range(len(sorted_stats)):
|
||||
|
@ -15,7 +15,7 @@ from cpython.exc cimport PyErr_CheckSignals
|
||||
from io import BufferedWriter
|
||||
|
||||
|
||||
__title__="Keep the N last lines of a view."
|
||||
__title__="Keep the N last lines of a view"
|
||||
|
||||
|
||||
def addOptions(parser):
|
||||
|
@ -39,7 +39,7 @@ COL_COMMENTS_MAX_LEN = 2048
|
||||
MAX_INT = 2147483647 # used to generate random float values
|
||||
|
||||
|
||||
__title__="Tests if the obitools are working properly"
|
||||
__title__="Test if the obitools are working properly"
|
||||
|
||||
|
||||
default_config = {
|
||||
|
@ -53,6 +53,8 @@ cdef extern from "obitypes.h" nogil:
|
||||
extern const_char_p OBIQual_char_NA
|
||||
extern uint8_t* OBIQual_int_NA
|
||||
extern void* OBITuple_NA
|
||||
|
||||
extern obiint_t OBI_INT_MAX
|
||||
|
||||
const_char_p name_data_type(int data_type)
|
||||
|
||||
|
@ -77,7 +77,7 @@ cdef class View(OBIWrapper) :
|
||||
|
||||
|
||||
@staticmethod
|
||||
def import_view(object dms_1, object dms_2, object view_name_1, object view_name_2):
|
||||
def import_view(object dms_1, object dms_2, object view_name_1, object view_name_2): # TODO argument to import line by line to avoid huge AVL copy
|
||||
if obi_import_view(tobytes(dms_1), tobytes(dms_2), tobytes(view_name_1), tobytes(view_name_2)) < 0 :
|
||||
raise Exception("Error importing a view")
|
||||
|
||||
@ -600,7 +600,8 @@ cdef class View(OBIWrapper) :
|
||||
if element is not None:
|
||||
if element.comments[b"input_dms_name"] is not None :
|
||||
for i in range(len(element.comments[b"input_dms_name"])) :
|
||||
if element.comments[b"input_dms_name"][i] == element.dms.name and b"/" not in element.comments[b"input_view_name"][i]: # Same DMS and not a special element like a taxonomy
|
||||
if b"/" not in element.comments[b"input_view_name"][i] and element.comments[b"input_view_name"][i] in element.dms \
|
||||
and element.comments[b"input_dms_name"][i] == element.dms.name : # Same DMS and not a special element like a taxonomy and view was not deleted
|
||||
top_level.append(element.dms[element.comments[b"input_view_name"][i]])
|
||||
else:
|
||||
top_level.append(None)
|
||||
|
@ -192,7 +192,7 @@ def open_uri(uri,
|
||||
|
||||
config = getConfiguration()
|
||||
urip = urlparse(urib)
|
||||
|
||||
|
||||
if 'obi' not in config:
|
||||
config['obi']={}
|
||||
|
||||
@ -209,13 +209,14 @@ def open_uri(uri,
|
||||
scheme = urip.scheme
|
||||
|
||||
error = None
|
||||
|
||||
if urib != b"-" and \
|
||||
|
||||
if b'/taxonomy/' in urib or \
|
||||
(urib != b"-" and \
|
||||
(scheme==b"dms" or \
|
||||
(scheme==b"" and \
|
||||
(((not input) and "outputformat" not in config["obi"]) or \
|
||||
(input and "inputformat" not in config["obi"])))): # TODO maybe not best way
|
||||
|
||||
(input and "inputformat" not in config["obi"]))))): # TODO maybe not best way
|
||||
|
||||
if default_dms is not None and b"/" not in urip.path: # assuming view to open in default DMS (TODO discuss)
|
||||
dms=(default_dms, urip.path)
|
||||
else:
|
||||
@ -223,7 +224,7 @@ def open_uri(uri,
|
||||
|
||||
if dms is None and default_dms is not None:
|
||||
dms=(default_dms, urip.path)
|
||||
|
||||
|
||||
if dms is not None:
|
||||
if dms_only:
|
||||
return (dms[0],
|
||||
@ -248,7 +249,7 @@ def open_uri(uri,
|
||||
|
||||
if default_dms is None:
|
||||
config["obi"]["defaultdms"]=resource[0]
|
||||
|
||||
|
||||
return (resource[0],
|
||||
resource[1],
|
||||
type(resource[1]),
|
||||
@ -386,10 +387,10 @@ def open_uri(uri,
|
||||
raise MalformedURIException('Malformed header argument in URI')
|
||||
|
||||
if b"sep" in qualifiers:
|
||||
sep=tobytes(qualifiers[b"sep"][0][0])
|
||||
sep = tobytes(qualifiers[b"sep"][0][0])
|
||||
else:
|
||||
try:
|
||||
sep=tobytes(config["obi"]["sep"])
|
||||
sep = tobytes(config["obi"]["sep"])
|
||||
except KeyError:
|
||||
sep=None
|
||||
|
||||
@ -464,7 +465,7 @@ def open_uri(uri,
|
||||
if format is not None:
|
||||
if seqtype==b"nuc":
|
||||
objclass = Nuc_Seq # Nuc_Seq_Stored? TODO
|
||||
if format==b"fasta" or format==b"silva":
|
||||
if format==b"fasta" or format==b"silva" or format==b"rdp":
|
||||
if input:
|
||||
iseq = fastaNucIterator(file,
|
||||
skip=skip,
|
||||
@ -565,6 +566,6 @@ def open_uri(uri,
|
||||
|
||||
entry_count = -1
|
||||
if input:
|
||||
entry_count = count_entries(file, format)
|
||||
entry_count = count_entries(file, format, header)
|
||||
|
||||
return (file, iseq, objclass, urib, entry_count)
|
||||
|
@ -3,7 +3,7 @@
|
||||
from obitools3.dms.capi.obitypes cimport obitype_t, index_t
|
||||
|
||||
cpdef bytes format_uniq_pattern(bytes format)
|
||||
cpdef int count_entries(file, bytes format)
|
||||
cpdef int count_entries(file, bytes format, bint header)
|
||||
|
||||
cdef obi_errno_to_exception(index_t line_nb=*, object elt_id=*, str error_message=*)
|
||||
|
||||
@ -18,7 +18,7 @@ cdef object clean_empty_values_from_object(object value, exclude=*)
|
||||
|
||||
cdef obitype_t get_obitype_single_value(object value)
|
||||
cdef obitype_t update_obitype(obitype_t obitype, object new_value)
|
||||
cdef obitype_t get_obitype_iterable_value(object value)
|
||||
cdef obitype_t get_obitype_iterable_value(object value, type t)
|
||||
cdef obitype_t get_obitype(object value)
|
||||
|
||||
cdef object __etag__(bytes x, bytes nastring=*)
|
||||
|
@ -9,7 +9,8 @@ from obitools3.dms.capi.obitypes cimport is_a_DNA_seq, \
|
||||
OBI_QUAL, \
|
||||
OBI_SEQ, \
|
||||
OBI_STR, \
|
||||
index_t
|
||||
index_t, \
|
||||
OBI_INT_MAX
|
||||
|
||||
from obitools3.dms.capi.obierrno cimport OBI_LINE_IDX_ERROR, \
|
||||
OBI_ELT_IDX_ERROR, \
|
||||
@ -39,7 +40,7 @@ cpdef bytes format_uniq_pattern(bytes format):
|
||||
return None
|
||||
|
||||
|
||||
cpdef int count_entries(file, bytes format):
|
||||
cpdef int count_entries(file, bytes format, bint header):
|
||||
|
||||
try:
|
||||
sep = format_uniq_pattern(format)
|
||||
@ -74,6 +75,8 @@ cpdef int count_entries(file, bytes format):
|
||||
total_count += len(re.findall(sep, mmapped_file))
|
||||
if format != b"ngsfilter" and format != b"tabular" and format != b"embl" and format != b"genbank" and format != b"fastq":
|
||||
total_count += 1 # adding +1 for 1st entry because separators include \n (ngsfilter and tabular already count one more because of last \n)
|
||||
if format == b"tabular" and header: # not counting header as an entry
|
||||
total_count -= 1
|
||||
|
||||
except:
|
||||
if len(files) > 1:
|
||||
@ -257,38 +260,51 @@ cdef obitype_t update_obitype(obitype_t obitype, object new_value) :
|
||||
|
||||
new_type = type(new_value)
|
||||
|
||||
if obitype == OBI_INT :
|
||||
if new_type == float :
|
||||
return OBI_FLOAT
|
||||
#if new_type == NoneType: # doesn't work because Cython sucks
|
||||
if new_value == None or new_type==list or new_type==dict or new_type==tuple:
|
||||
return obitype
|
||||
|
||||
# TODO BOOL vers INT/FLOAT
|
||||
elif new_type == str or new_type == bytes :
|
||||
if new_type == str or new_type == bytes :
|
||||
if obitype == OBI_SEQ and is_a_DNA_seq(tobytes(new_value)) :
|
||||
pass
|
||||
else :
|
||||
return OBI_STR
|
||||
|
||||
elif obitype == OBI_INT :
|
||||
if new_type == float or new_value > OBI_INT_MAX :
|
||||
return OBI_FLOAT
|
||||
|
||||
return obitype
|
||||
|
||||
|
||||
cdef obitype_t get_obitype_iterable_value(object value) :
|
||||
cdef obitype_t get_obitype_iterable_value(object value, type t) :
|
||||
|
||||
cdef obitype_t value_obitype
|
||||
|
||||
value_obitype = OBI_VOID
|
||||
|
||||
for k in value :
|
||||
if value_obitype == OBI_VOID :
|
||||
value_obitype = get_obitype_single_value(value[k])
|
||||
else :
|
||||
value_obitype = update_obitype(value_obitype, value[k])
|
||||
if t == dict:
|
||||
for k in value :
|
||||
if value_obitype == OBI_VOID :
|
||||
value_obitype = get_obitype_single_value(value[k])
|
||||
else :
|
||||
value_obitype = update_obitype(value_obitype, value[k])
|
||||
|
||||
elif t == list or t == tuple:
|
||||
for v in value :
|
||||
if value_obitype == OBI_VOID :
|
||||
value_obitype = get_obitype_single_value(v)
|
||||
else :
|
||||
value_obitype = update_obitype(value_obitype, v)
|
||||
|
||||
return value_obitype
|
||||
|
||||
|
||||
cdef obitype_t get_obitype(object value) :
|
||||
|
||||
if type(value) == dict or type(value) == list or type(value) == tuple :
|
||||
return get_obitype_iterable_value(value)
|
||||
t = type(value)
|
||||
if t == dict or t == list or t == tuple :
|
||||
return get_obitype_iterable_value(value, t)
|
||||
|
||||
else :
|
||||
return get_obitype_single_value(value)
|
||||
|
@ -1,5 +1,5 @@
|
||||
major = 3
|
||||
minor = 0
|
||||
serial= '1b4'
|
||||
serial= '1b9'
|
||||
|
||||
version ="%d.%d.%s" % (major,minor,serial)
|
||||
|
@ -863,7 +863,8 @@ int build_reference_db(const char* dms_name,
|
||||
fprintf(stderr,"\rDone : 100 %% \n");
|
||||
|
||||
// Add information about the threshold used to build the DB
|
||||
snprintf(threshold_str, 5, "%f", threshold);
|
||||
#define snprintf_nowarn(...) (snprintf(__VA_ARGS__) < 0 ? abort() : (void)0)
|
||||
snprintf_nowarn(threshold_str, 5, "%f", threshold);
|
||||
|
||||
new_comments = obi_add_comment((o_view->infos)->comments, DB_THRESHOLD_KEY_IN_COMMENTS, threshold_str);
|
||||
if (new_comments == NULL)
|
||||
|
@ -77,7 +77,6 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
int32_t* shift_count_array;
|
||||
Obi_ali_p ali = NULL;
|
||||
int i, j;
|
||||
bool switched_seqs;
|
||||
bool reversed;
|
||||
int score = 0;
|
||||
Obi_blob_p blob1 = NULL;
|
||||
@ -124,6 +123,8 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
bool keep_seq2_start;
|
||||
bool keep_seq1_end;
|
||||
bool keep_seq2_end;
|
||||
bool left_ali;
|
||||
bool rev_quals = false;
|
||||
|
||||
// Check kmer size
|
||||
if ((kmer_size < 1) || (kmer_size > 4))
|
||||
@ -148,19 +149,8 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
}
|
||||
|
||||
// Choose the shortest sequence to save kmer positions in array
|
||||
switched_seqs = false;
|
||||
len1 = blob1->length_decoded_value;
|
||||
len2 = blob2->length_decoded_value;
|
||||
if (len2 < len1)
|
||||
{
|
||||
switched_seqs = true;
|
||||
temp_len = len1;
|
||||
len1 = len2;
|
||||
len2 = temp_len;
|
||||
temp_blob = blob1;
|
||||
blob1 = blob2;
|
||||
blob2 = temp_blob;
|
||||
}
|
||||
|
||||
// Force encoding on 2 bits by replacing ambiguous nucleotides by 'a's
|
||||
if (blob1->element_size == 4)
|
||||
@ -196,7 +186,47 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
else
|
||||
reversed = false;
|
||||
if (reversed)
|
||||
switched_seqs = !switched_seqs;
|
||||
// unreverse to make cases simpler. Costly but rare (direct match is reverse primer match)
|
||||
{
|
||||
if (seq2 == NULL)
|
||||
seq2 = obi_blob_to_seq(blob2);
|
||||
seq2 = reverse_complement_sequence(seq2);
|
||||
blob2 = obi_seq_to_blob(seq2);
|
||||
|
||||
if (seq1 == NULL)
|
||||
seq1 = obi_blob_to_seq(blob1);
|
||||
seq1 = reverse_complement_sequence(seq1);
|
||||
blob1 = obi_seq_to_blob(seq1);
|
||||
free_blob1 = true;
|
||||
|
||||
// Get the quality arrays
|
||||
qual1 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view1, qual_col1, idx1, 0, &qual1_len);
|
||||
if (qual1 == NULL)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 1st sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
}
|
||||
qual2 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view2, qual_col2, idx2, 0, &qual2_len);
|
||||
if (qual2 == NULL)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 2nd sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
}
|
||||
|
||||
uint8_t* newqual1 = malloc(qual1_len*sizeof(uint8_t));
|
||||
uint8_t* newqual2 = malloc(qual2_len*sizeof(uint8_t));
|
||||
|
||||
for (i=0;i<qual1_len;i++)
|
||||
newqual1[i] = qual1[qual1_len-1-i];
|
||||
|
||||
for (i=0;i<qual2_len;i++)
|
||||
newqual2[i] = qual2[qual2_len-1-i];
|
||||
|
||||
qual1 = newqual1;
|
||||
qual2 = newqual2;
|
||||
|
||||
rev_quals = true;
|
||||
}
|
||||
|
||||
// Save total length for the shift counts array
|
||||
total_len = len1 + len2 + 1; // +1 for shift 0
|
||||
@ -237,7 +267,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
return NULL;
|
||||
}
|
||||
}
|
||||
else if (len1 >= shift_array_height)
|
||||
else if (total_len >= shift_array_height)
|
||||
{
|
||||
shift_array_height = total_len;
|
||||
*shift_array_p = (int32_t*) realloc(*shift_array_p, ARRAY_LENGTH * shift_array_height * sizeof(int32_t));
|
||||
@ -291,7 +321,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
*shift_array_height_p = shift_array_height;
|
||||
*shift_count_array_length_p = shift_count_array_length;
|
||||
|
||||
// Fill array with positions of kmers in the shortest sequence
|
||||
// Fill array with positions of kmers in the first sequence
|
||||
encoding = blob1->element_size;
|
||||
kmer_count = len1 - kmer_size + 1;
|
||||
for (kmer_idx=0; kmer_idx < kmer_count; kmer_idx++)
|
||||
@ -310,7 +340,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
kmer_pos_array[(kmer*kmer_pos_array_height)+i] = kmer_idx;
|
||||
}
|
||||
|
||||
// Compare positions of kmers between both sequences and store shifts
|
||||
// Compare positions of kmers between both sequences and store shifts (a shift corresponds to pos2 - pos1)
|
||||
kmer_count = blob2->length_decoded_value - kmer_size + 1;
|
||||
for (kmer_idx=0; kmer_idx < kmer_count; kmer_idx++)
|
||||
{
|
||||
@ -374,35 +404,42 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
// The 873863 cases of hell
|
||||
if (best_shift > 0)
|
||||
{
|
||||
left_ali = false;
|
||||
overlap_len = len2 - best_shift;
|
||||
if (len1 <= overlap_len)
|
||||
{
|
||||
overlap_len = len1;
|
||||
if (! switched_seqs)
|
||||
keep_seq2_end = true;
|
||||
else
|
||||
keep_seq2_start = true;
|
||||
}
|
||||
else if (switched_seqs)
|
||||
{
|
||||
keep_seq2_start = true;
|
||||
keep_seq1_end = true;
|
||||
keep_seq2_end = true;
|
||||
}
|
||||
}
|
||||
else if (best_shift < 0)
|
||||
{
|
||||
left_ali = true;
|
||||
overlap_len = len1 + best_shift;
|
||||
if (!switched_seqs)
|
||||
{
|
||||
keep_seq1_start = true;
|
||||
keep_seq2_end = true;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
overlap_len = len1;
|
||||
if ((!switched_seqs) && (len2 > len1))
|
||||
if (len2 <= overlap_len)
|
||||
{
|
||||
overlap_len = len2;
|
||||
keep_seq1_start = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
keep_seq1_start = true;
|
||||
keep_seq2_end = true;
|
||||
}
|
||||
}
|
||||
else // if (best_shift == 0)
|
||||
{
|
||||
if (len2 >= len1)
|
||||
{
|
||||
overlap_len = len1;
|
||||
keep_seq2_end = true;
|
||||
left_ali = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
overlap_len = len2;
|
||||
left_ali = false; // discussable
|
||||
}
|
||||
}
|
||||
|
||||
ali = (Obi_ali_p) malloc(sizeof(Obi_ali_t));
|
||||
@ -433,7 +470,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
ali->direction[0] = '\0';
|
||||
else
|
||||
{
|
||||
if (((best_shift <= 0) && (!switched_seqs)) || ((best_shift > 0) && switched_seqs))
|
||||
if (left_ali)
|
||||
strcpy(ali->direction, "left");
|
||||
else
|
||||
strcpy(ali->direction, "right");
|
||||
@ -442,28 +479,28 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
// Build the consensus sequence if asked
|
||||
if (build_consensus)
|
||||
{
|
||||
// Get the quality arrays
|
||||
qual1 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view1, qual_col1, idx1, 0, &qual1_len);
|
||||
if (qual1 == NULL)
|
||||
if (! rev_quals)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 1st sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
}
|
||||
qual2 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view2, qual_col2, idx2, 0, &qual2_len);
|
||||
if (qual2 == NULL)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 2nd sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
// Get the quality arrays
|
||||
qual1 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view1, qual_col1, idx1, 0, &qual1_len);
|
||||
if (qual1 == NULL)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 1st sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
}
|
||||
qual2 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view2, qual_col2, idx2, 0, &qual2_len);
|
||||
if (qual2 == NULL)
|
||||
{
|
||||
obidebug(1, "\nError getting the quality of the 2nd sequence when computing the kmer similarity between two sequences");
|
||||
return NULL;
|
||||
}
|
||||
}
|
||||
|
||||
// Decode the first sequence if not already done
|
||||
if (seq1 == NULL)
|
||||
seq1 = obi_blob_to_seq(blob1);
|
||||
|
||||
if (! switched_seqs)
|
||||
consensus_len = len2 - best_shift;
|
||||
else
|
||||
consensus_len = len1 + best_shift;
|
||||
consensus_len = len2 - best_shift;
|
||||
|
||||
// Allocate memory for consensus sequence
|
||||
consensus_seq = (char*) malloc(consensus_len + 1 * sizeof(char)); // TODO keep malloced too maybe
|
||||
@ -557,6 +594,12 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
|
||||
free(seq2);
|
||||
free(blob2);
|
||||
|
||||
if (rev_quals)
|
||||
{
|
||||
free(qual1);
|
||||
free(qual2);
|
||||
}
|
||||
|
||||
return ali;
|
||||
}
|
||||
|
||||
|
@ -233,7 +233,7 @@ int obi_ecotag(const char* dms_name,
|
||||
// Write result (max score, threshold, LCA assigned, list of the ids of the best matches)
|
||||
|
||||
|
||||
index_t i, j, k;
|
||||
index_t i, j, k, t;
|
||||
ecotx_t* lca;
|
||||
ecotx_t* lca_in_array;
|
||||
ecotx_t* best_match;
|
||||
@ -259,16 +259,20 @@ int obi_ecotag(const char* dms_name,
|
||||
int32_t* best_match_taxids;
|
||||
int32_t* best_match_taxids_to_store;
|
||||
int best_match_count;
|
||||
int best_match_taxid_count;
|
||||
int buffer_size;
|
||||
int best_match_ids_buffer_size;
|
||||
index_t best_match_idx;
|
||||
int32_t lca_array_length;
|
||||
int32_t lca_taxid;
|
||||
int32_t taxid_best_match;
|
||||
int32_t taxid;
|
||||
int32_t taxid_to_store;
|
||||
bool assigned;
|
||||
const char* lca_name;
|
||||
const char* id;
|
||||
int id_len;
|
||||
bool already_in;
|
||||
|
||||
OBIDMS_p dms = NULL;
|
||||
OBIDMS_p ref_dms = NULL;
|
||||
@ -488,10 +492,11 @@ int obi_ecotag(const char* dms_name,
|
||||
|
||||
for (i=0; i < query_count; i++)
|
||||
{
|
||||
if (i%1000 == 0)
|
||||
if (i%10 == 0)
|
||||
fprintf(stderr,"\rDone : %f %% ", (i / (float) query_count)*100);
|
||||
|
||||
best_match_count = 0;
|
||||
best_match_taxid_count = 0;
|
||||
best_match_ids_length = 0;
|
||||
threshold = ecotag_threshold;
|
||||
best_score = 0.0;
|
||||
@ -543,6 +548,7 @@ int obi_ecotag(const char* dms_name,
|
||||
// Reset the array with that match
|
||||
best_match_ids_length = 0;
|
||||
best_match_count = 0;
|
||||
best_match_taxid_count = 0;
|
||||
}
|
||||
|
||||
// Store in best match array
|
||||
@ -585,8 +591,27 @@ int obi_ecotag(const char* dms_name,
|
||||
|
||||
// Save match
|
||||
best_match_array[best_match_count] = j;
|
||||
best_match_taxids[best_match_count] = obi_get_int_with_elt_idx_and_col_p_in_view(ref_view, ref_taxid_column, j, 0);
|
||||
best_match_count++;
|
||||
|
||||
// Save best match taxid only if not already in array
|
||||
taxid_to_store = obi_get_int_with_elt_idx_and_col_p_in_view(ref_view, ref_taxid_column, j, 0);
|
||||
already_in = false;
|
||||
for (t=0; t<best_match_taxid_count; t++)
|
||||
{
|
||||
taxid = best_match_taxids[t];
|
||||
//fprintf(stderr, "\ntaxid %d, taxid_to_store %d\n", taxid, taxid_to_store);
|
||||
if (taxid == taxid_to_store)
|
||||
{
|
||||
already_in = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (! already_in)
|
||||
{
|
||||
best_match_taxids[best_match_taxid_count] = taxid_to_store;
|
||||
best_match_taxid_count++;
|
||||
}
|
||||
|
||||
strcpy(best_match_ids+best_match_ids_length, id);
|
||||
best_match_ids_length = best_match_ids_length + id_len + 1;
|
||||
}
|
||||
@ -693,7 +718,7 @@ int obi_ecotag(const char* dms_name,
|
||||
assigned_name_column, lca_name,
|
||||
assigned_status_column, assigned,
|
||||
best_match_ids_column, best_match_ids_to_store, best_match_ids_length,
|
||||
best_match_taxids_column, best_match_taxids_to_store, best_match_count,
|
||||
best_match_taxids_column, best_match_taxids_to_store, best_match_taxid_count,
|
||||
score_column, best_score
|
||||
) < 0)
|
||||
return -1;
|
||||
|
@ -2463,6 +2463,32 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
|
||||
return -1;
|
||||
}
|
||||
|
||||
// Write the rest of the taxa from the current taxa list
|
||||
while (nT < (tax->taxa)->count)
|
||||
{
|
||||
// Add element from taxa list
|
||||
// Enlarge structure if needed
|
||||
if (n == buffer_size)
|
||||
{
|
||||
buffer_size = buffer_size * 2;
|
||||
tax->merged_idx = (ecomergedidx_t*) realloc(tax->merged_idx, sizeof(ecomergedidx_t) + sizeof(ecomerged_t) * buffer_size);
|
||||
if (tax->merged_idx == NULL)
|
||||
{
|
||||
obi_set_errno(OBI_MALLOC_ERROR);
|
||||
obidebug(1, "\nError reallocating memory for a taxonomy structure");
|
||||
closedir(tax_dir);
|
||||
fclose(file);
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
|
||||
(tax->merged_idx)->merged[n].taxid = (tax->taxa)->taxon[nT].taxid;
|
||||
(tax->merged_idx)->merged[n].idx = nT;
|
||||
|
||||
nT++;
|
||||
n++;
|
||||
}
|
||||
|
||||
// Store count
|
||||
(tax->merged_idx)->count = n;
|
||||
|
||||
|
@ -75,7 +75,7 @@ typedef struct {
|
||||
*/
|
||||
int32_t max_taxid; /**< Maximum taxid existing in the taxon index.
|
||||
*/
|
||||
int32_t buffer_size; /**< Number of taxa. // TODO kept this but not sure of its use
|
||||
int32_t buffer_size; /**< . // TODO kept this but not sure of its use
|
||||
*/
|
||||
ecotx_t taxon[]; /**< Taxon array.
|
||||
*/
|
||||
|
@ -36,7 +36,7 @@
|
||||
*/
|
||||
#define COLUMN_GROWTH_FACTOR (2) /**< The growth factor when a column is enlarged.
|
||||
*/
|
||||
#define MAXIMUM_LINE_COUNT (1000000000) /**< The maximum line count for the data of a column (1E9). //TODO
|
||||
#define MAXIMUM_LINE_COUNT (1000000000000) /**< The maximum line count for the data of a column (1E12). //TODO
|
||||
*/
|
||||
#define COMMENTS_MAX_LENGTH (4096) /**< The maximum length for comments.
|
||||
*/
|
||||
|
@ -29,6 +29,8 @@
|
||||
#define OBIQual_int_NA (NULL) /**< NA value for the type OBI_QUAL if the quality is in integer format */
|
||||
#define OBITuple_NA (NULL) /**< NA value for tuples of any type */
|
||||
|
||||
#define OBI_INT_MAX (INT32_MAX) /**< Maximum value for the type OBI_INT */
|
||||
|
||||
|
||||
/**
|
||||
* @brief enum for the boolean OBIType.
|
||||
|
Reference in New Issue
Block a user