Compare commits

...

140 Commits

Author SHA1 Message Date
d27a5b9115 Switch to version 3.0.0b37 2020-10-30 10:47:13 +01:00
20bd3350b4 New command: obi addtaxids to add NCBI taxids to sequences from their
taxon name.
2020-10-30 10:46:55 +01:00
2e191372d7 Now handling sequences with Uracil (U) nucleotides by converting to
Thymine (T)
2020-10-30 10:46:17 +01:00
112e12cab0 Taxonomy: new functions to find taxa by name 2020-10-30 10:45:20 +01:00
b9b4cec5b5 import: now can import SILVA fasta files 2020-10-30 10:43:04 +01:00
199f3772e8 Small fixes (potential compilation problems) 2020-10-30 10:41:58 +01:00
422a6450fa ecotag: clarified similarity circle documentation 2020-09-29 17:57:29 +02:00
137c109f86 obi ls: now done in C (preparing things for R packages to read DMS) and
switch to version 3.0.0b36
2020-09-29 17:51:39 +02:00
b6648ae81e Revert "Fixed version numbering mistake (should be b34 not b35)"
This reverts commit f6dffbecfe
2020-09-25 16:25:39 +02:00
f6dffbecfe Fixed version numbering mistake (should be b34 not b35) 2020-09-25 16:24:23 +02:00
c4696ac865 ecotag: added separate threshold for minimum circle identity (and switch
to version 3.0.0b35
2020-09-25 16:22:09 +02:00
11a0945a9b obi cat: fixed open file descriptor leak and switch to version 3.0.0b34 2020-08-28 10:41:22 +02:00
f23c40c905 obi cat: fixed a bug introduced in 3.0.0b28 and switch to version
3.0.0b33
2020-08-27 18:38:16 +02:00
f99fc13b75 switch to version 3.0.0b32 2020-08-13 18:17:09 +02:00
1da6aac1b8 C: patch for failed creation of AVL with errno EEXIST 2020-08-12 17:55:08 +02:00
159803b40a export: now automatically sorts dictionary keys alphabetically for
tab/csv output
2020-07-31 16:43:35 +02:00
7dcbc34017 import: fixed entry count estimation when importing fastq files 2020-07-30 16:56:36 +02:00
db2202c8b4 uniq: added a check to make sure that there is more than one element for
one tag when merging its information
2020-07-30 16:14:37 +02:00
d33ff97846 switch to version 3.0.0b31 2020-07-28 09:31:19 +02:00
1dcdf69f1f export: fixed a bug introduced in version 3.0.0b28 2020-07-28 09:31:05 +02:00
dec114eed6 Python: added "date created" information in view representation 2020-07-27 17:38:45 +02:00
f36691053b Python: added the OBITools3 version that generated the view in view
comments
2020-07-27 16:50:00 +02:00
f2aa5fcf8b alignpairedend: fixed division by 0 bug and switch to version 3.0.0b30 2020-07-27 10:15:59 +02:00
bccb3e6874 switch to version 3.0.0b29 2020-07-26 17:40:26 +02:00
f5a17bea68 C: added a missing error check 2020-07-26 17:39:55 +02:00
e28507639a C and Cython: fixed and improved the associated columns system 2020-07-26 17:39:29 +02:00
e6feac93fe obi test: made less heavy to be faster 2020-07-26 17:37:21 +02:00
50b292b489 obi import: added --space-priority option to import a view line by line 2020-07-26 17:36:52 +02:00
24a737aa55 switch to version 3.0.0b28 2020-07-24 16:10:10 +02:00
8aa455ad8a Python: made all commands handle output to buffer object (e.g. stdout) 2020-07-24 16:09:48 +02:00
46ca693ca9 Cython: View: new method to print a view to a buffer (e.g. stdout) 2020-07-24 16:03:23 +02:00
9a9afde113 Cython: progress bar: set default refresh rate to 5 seconds 2020-07-24 15:29:12 +02:00
8dd403a118 grep: now prints the number of entries grepped 2020-07-13 17:08:13 +02:00
9672f01c6a alignpairedend: improved/fixed the output tags for the alignment score
and lengths. Removed minimum score option
2020-07-13 15:59:50 +02:00
ed9549acfb ngsfilter: unidentified sequences are now stored untrimmed 2020-07-13 15:56:40 +02:00
9ace9989c4 Switch to version 3.0.0b27 2020-07-07 16:47:21 +02:00
a3ebe5f118 C: AVL trees: fixed a bug where storing the difference between 2 crc64
values in an int64 would mess trees up resulting in failed data
dereplication
2020-07-07 16:47:00 +02:00
9100e14899 obi uniq: quick fix for bug where some sequences are not correctly
dereplicated
2020-07-03 17:36:57 +02:00
ccda0661ce small help documentation improvement 2020-07-01 18:20:38 +02:00
aab59f2214 obi clean: fixed a memory bug, fixed the behaviour when no sample info,
and added checks warnings and error handling when sample info not
dereplicated
2020-07-01 18:17:47 +02:00
ade1107b42 switch to version 3.0.0b26 2020-06-17 18:56:07 +02:00
9c7d24406f export: dictionaries are now formatted like in the original OBITools
when exporting in tabular format and tuple formatting is cleaner
2020-06-17 18:55:46 +02:00
03bc9915f2 Cython: utils: added handling of tuples to bytes2str_object function 2020-06-17 18:54:14 +02:00
24b1dab573 Cython: Columns: added a keys() method that returns all element names 2020-06-17 18:53:41 +02:00
7593673f3f ngsfilter: now setting 'reversed' tag to False instead of None when
false
2020-06-17 18:52:35 +02:00
aa01236cae switch to version 3.0.0b25 2020-06-13 21:48:49 +02:00
49b8810a76 C: made indexer opening/closing cleaner 2020-06-13 21:47:03 +02:00
7a39df54c0 ls: fixed an issue where big DMS couldn't be read by ls 2020-06-13 21:45:22 +02:00
09e483b0d6 switch to temporary version 3.0.0b24a 2020-06-10 17:47:56 +02:00
14a2579173 uniq: now outputs an empty view if input view is empty instead of
displaying an error
2020-06-10 17:47:26 +02:00
36a8aaa92e grep: now creating empty views instead of displaying an error when
selecting on an unexisting column/tag
2020-06-10 16:57:42 +02:00
a17eb445c2 ngsfilter: made one of the tag error messages more accurate 2020-06-10 16:27:36 +02:00
e4a32788c2 Switch to version 3.0.0b24 2020-06-09 14:36:58 +02:00
2442cc80bf Cython: View: fixed bash history display 2020-06-09 14:36:37 +02:00
aa836b2ace uniq: improved progress bar of second browsing 2020-06-09 14:36:02 +02:00
8776ce22e6 C: fixed a bug where indexers referring to tuples of certain types were
not properly closed and imported
2020-06-09 14:34:43 +02:00
4aa772c405 ecotag: Added list of taxids for all best matches (closes #80) 2020-06-09 14:33:14 +02:00
b0b96ac37a version 3.0.0b23a 2020-06-05 16:10:24 +02:00
687e42ad22 C: kmer alignment: fixed a bug where scores of 0 were at
(0+kmer_length-1) (and now setting alignment direction to None if score
is 0
2020-06-05 16:09:33 +02:00
5fbbb6d304 alignpairedend: fixed a bug when rebuilding joined (unaligned) sequences
where only the forward sequence was kept
2020-06-05 16:06:43 +02:00
359a9fe237 Switch to version 3.0.0b23 2020-06-04 15:35:03 +02:00
f9b6851f75 Python: correctly flagged some mandatory options as required 2020-06-04 15:34:24 +02:00
29a2652bbf Fixed installation on Ubuntu without pip 2020-06-04 15:06:35 +02:00
2a2c233936 obi import: fixed a bug when skipping an entry 2020-05-29 21:19:42 +02:00
faf8ea9d86 Switch to version 3.0.0b21 2020-05-28 20:42:09 +02:00
ffe2485e94 Genbank parser: now reading ORIGIN lines with comments without
triggering error
2020-05-28 20:41:34 +02:00
6094ce2bbc obi import: skip on error more robust 2020-05-28 20:40:36 +02:00
a7dcf16c06 Minor changes for pip release 2020-05-20 15:59:04 +02:00
f13f8f6165 obi import: minor doc/display improvements 2020-05-20 11:46:29 +02:00
b5a29ac413 Switch to version 3.0.0b19 2020-05-20 10:29:36 +02:00
efd2b9d338 Cleaner installation 2020-05-20 10:29:12 +02:00
ca6e3e7aad obi import: fixed to work with seq genbank extension 2020-05-20 10:28:14 +02:00
76ed8e18e5 Switch to version 3.0.0b18 with version formatting that fits setuptools 2020-05-18 17:08:55 +02:00
1d17f28aec setup: now using setuptools instead of distutils to work with pip 2020-05-18 17:08:09 +02:00
fa834e4b8b obi import: small bug fix 2020-05-18 17:06:58 +02:00
a72fea3cc9 Python: fasta parser: fixed a bug stopping the program when the last
line contained a single nucleotide
2020-05-12 11:24:12 +02:00
e9a37d8a6e Switch to version 3.0.0-beta16 2020-05-07 17:09:26 +02:00
ef074f8455 typo 2020-05-07 17:08:59 +02:00
aec5e69f2c C, views: no more automatic COUNT column if MERGED_sample column exists 2020-05-07 17:08:07 +02:00
170ef3f1ba Views: added obi prefix to commands in bash history 2020-05-07 17:07:01 +02:00
f999946582 obi uniq: fixed the remerging of already merged informations, and
efficiency improvements
2020-05-07 17:05:54 +02:00
773b36ec37 obi import: fixed the import of old obitools files with premerged
informations, and other minor improvements
2020-05-07 17:03:04 +02:00
69cb434a6c version 3.0.0-beta15c 2020-04-29 14:25:33 +02:00
55d4f98d60 obi annotate: fixed annotation at ranks 2020-04-29 14:24:40 +02:00
0bec2631e8 ecotag: fixed a bug where all the full DMS path weren't properly sent to
the C layer
2020-04-29 10:35:55 +02:00
e6b6c6fa84 AVLs: Made an error message more informative 2020-04-29 10:14:04 +02:00
974528b2e6 build_ref_db: fixed bug erasing some of the higher LCAs (i.e. lowest
similarities)
2020-04-28 15:56:06 +02:00
1b346b54f9 ecotag: better specificity by now correctly looking for similarities
within refs above best score instead of ecotag threshold
2020-04-28 15:10:07 +02:00
058f2ad8b3 ecopcr: fixed a bug where sequences were considered circular (generating
false positives)
2020-04-27 14:44:35 +02:00
60bfd3ae8d obi annotate: now defaults to setting str if expression is not valid 2020-04-24 11:35:20 +02:00
67bdee105a C: build_ref_db: added progress display for each step 2020-04-18 14:24:08 +02:00
0f745e0113 C: Columns: optimizing column file growth 2020-04-18 13:55:47 +02:00
da8de52ba4 export: fixed progress bar bug 2020-04-17 15:09:10 +02:00
4d36538c6e C: SSE lcs alignment: band-aid for memory bug I don't understand
(triggered on specific db on ubuntu)
2020-04-17 15:07:52 +02:00
8d0b17d87d Switch to version 3.0.0-beta14 2020-04-15 17:47:26 +02:00
343999a627 Taxonomy: fixed a critical memory bug when building the list of merged
taxids
2020-04-15 17:46:13 +02:00
e9a40630e9 C: Columns: rounding column growth to ceil to avoid looping on small
values
2020-04-13 19:02:10 +02:00
8dbcd3025a C: Columns: reduced column growth factor from 2 to 1.3 to avoid errno28 2020-04-13 14:47:56 +02:00
4cf635d001 Switch to version 3.0.0-beta13 2020-04-12 17:42:58 +02:00
b7e7cc232a Made completion script cleaner 2020-04-12 17:41:59 +02:00
b6ab792ceb C: made error message more detailed when checking that sequences and
qualities match
2020-04-12 17:40:24 +02:00
ddea5a2964 obi import: fixed inconsequential error when precomputing number of
entries in some formats
2020-04-12 17:38:42 +02:00
30852ab7d5 View bash history: removed useless shebang 2020-04-12 17:36:04 +02:00
4d0299904e all commands (almost): cleaner DMS closing at the end 2020-04-12 17:31:58 +02:00
eef5156d95 obi stats: fixed error when printing bool keys 2020-04-12 17:12:04 +02:00
e62c991bbc goes with previous commit 2020-04-10 11:22:26 +02:00
1218eed7fd ecopcr: now printing a warning instead of interrupting with an error
when a taxid is not found
2020-04-10 11:22:04 +02:00
cd9cea8c97 obi import: fixed critical bug where the last entry of embl and genbank
files was not imported
2020-04-09 19:26:27 +02:00
98cfb70d73 ecopcr: made some errors more informative 2020-04-09 09:15:28 +02:00
b9f68c76c8 ecopcr: added warnings and check of primer length (related to #75) 2020-04-05 18:40:56 +02:00
0b98371688 ngsfilter: added warning about primer length in -h (#75) 2020-04-05 18:39:20 +02:00
f0d152fcbd ngsfilter: now checking primer length (fixes #75) 2020-04-05 18:29:10 +02:00
8019dee68e ecotag: now closing all DMS properly 2020-04-05 13:20:49 +02:00
0b4a234671 Swich to version 3.0.0-beta11 2020-02-12 14:23:42 +01:00
d32cfdcce5 ecotag: fixed the generated column comments formatting that would
generate errors
2020-02-12 14:23:17 +01:00
219c0d6fdc obi cat: Fixed the handling when concatenating views with dictionaries
having different key sets
2020-02-12 14:21:39 +01:00
dc9f897917 switch to version 3.0.0-beta10 2020-02-02 21:15:27 +01:00
bb72682f7d obi import: new option --preread to do a first readthrough of the
dataset if it contains huge dictionaries for a much faster import.
2020-02-02 21:12:34 +01:00
52920c3c71 URI decoding: dirty temp fix for bug where default dms makes a mess when
should guess file
2020-02-02 21:11:05 +01:00
18c22cecf9 switch to version 3.0.0-beta9 2020-02-01 15:48:55 +01:00
1bfb96023c obi import: rewriting a column now deletes the old one to save disk
space
2020-02-01 15:31:14 +01:00
c67d668989 obi import: fixed a bug when the first entry would contain a dictionary
with one key. Switch to beta8
2020-01-29 20:23:39 +01:00
db0ac37d41 switch to version 3.0.0-beta7 2020-01-29 16:18:53 +01:00
d0c21ecd39 Removed an OpenMP clause that was not obligatory and triggered a known
gcc bug involving macros
2020-01-24 16:00:53 +01:00
53212168a2 History: added 'obi' in bash history for practical reasons 2020-01-23 16:51:49 +01:00
b4b2e62195 Cleaner handling of reverse quality columns 2020-01-18 19:28:12 +01:00
ced82c4242 Switching to version 3.0-beta6 2020-01-18 17:29:23 +01:00
a524f8829e New command: obi cat to concatenate views (not optimized yet) 2020-01-18 17:28:31 +01:00
5c9091e9eb C: closing DMS after cleaning it instead of counting on upper layer 2020-01-18 17:27:35 +01:00
822000cb70 Fixes in documentation 2020-01-18 17:26:18 +01:00
b9cd9bee9a C: Changed obibool definitions because of conflict with R 2020-01-06 15:11:31 +01:00
b1f3e082f9 ngsfilter: fixed a bug when there is only one tag introduced in latest
edit
2020-01-06 13:53:38 +01:00
6c018b403c ecopcr: fixed and improved the options to keep nuclotides around the
amplicon
2019-12-26 20:45:54 +01:00
694d1934a8 Tagging version beta3 2019-12-12 17:03:13 +01:00
fc3ac03630 clean_dms: now works with extension 2019-12-12 17:02:50 +01:00
d75e54a078 uniq: added forced deletion of reverse sequence quality 2019-12-12 17:02:36 +01:00
6bfd7441f3 ngsfilter: fixed sequence cutting when dealing with unaligned sequences.
Could use optimization
2019-12-12 17:01:31 +01:00
81a179239c ngsfilter: fixed sequence cut bug on aligned sequences. Still exists for
unaligned sequences
2019-12-10 18:13:27 +01:00
35ce37c0f7 ngsfilter: fixed a bug with unaligned chimeras (unpaired primers) and
made error annotations more explicit
2019-12-10 13:43:32 +01:00
53f18316b0 ngsfilter: made more robust and practical to use with empty tags 2019-11-29 15:21:08 +01:00
78 changed files with 2972 additions and 940 deletions

View File

@ -6,7 +6,7 @@ recursive-include doc/sphinx/source *.txt *.rst *.py
recursive-include doc/sphinx/sphinxext *.py
include doc/sphinx/Makefile
include doc/sphinx/Doxyfile
include README.txt
include README.md
include requirements.txt
include scripts/obi

View File

@ -1,4 +1,3 @@
#/usr/bin/env bash
_obi_comp ()
{

View File

@ -39,6 +39,12 @@ def __addImportInputOption(optionManager):
const=b'fastq',
help="Input file is in fastq format")
group.add_argument('--silva-input',
action="store_const", dest="obi:inputformat",
default=None,
const=b'silva',
help="Input file is in SILVA fasta format")
group.add_argument('--embl-input',
action="store_const", dest="obi:inputformat",
default=None,
@ -55,7 +61,7 @@ def __addImportInputOption(optionManager):
action="store_const", dest="obi:inputformat",
default=None,
const=b'ngsfilter',
help="Input file is an ngsfilter file")
help="Input file is an ngsfilter file. If not using tags, use ':' or 'None:None' or '-:-' or any combination")
group.add_argument('--ecopcr-result-input',
action="store_const", dest="obi:inputformat",
@ -222,7 +228,7 @@ def __addDMSOutputOption(optionManager):
group.add_argument('--no-create-dms',
action="store_true", dest="obi:nocreatedms",
default=False,
help="Don't create an output DMS it does not already exist")
help="Don't create an output DMS if it does not already exist")
def __addEltLimitOption(optionManager):

View File

@ -30,12 +30,12 @@ cdef class ProgressBar:
off_t maxi,
dict config={},
str head="",
double seconde=0.1,
double seconds=5,
cut=False):
self.starttime = self.clock()
self.lasttime = self.starttime
self.tickcount = <clock_t> (seconde * CLOCKS_PER_SEC)
self.tickcount = <clock_t> (seconds * CLOCKS_PER_SEC)
self.freq = 1
self.cycle = 0
self.arrow = 0

View File

@ -0,0 +1,230 @@
#cython: language_level=3
from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.dms.column.column cimport Column
from functools import reduce
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes, tostr
from io import BufferedWriter
from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
ID_COLUMN, \
DEFINITION_COLUMN, \
QUALITY_COLUMN, \
COUNT_COLUMN, \
TAXID_COLUMN
from obitools3.dms.capi.obitypes cimport OBI_INT
from obitools3.dms.capi.obitaxonomy cimport MIN_LOCAL_TAXID
import time
import math
import sys
from cpython.exc cimport PyErr_CheckSignals
__title__="Annotate sequences with their corresponding NCBI taxid found from the taxon scientific name."
def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi addtaxids specific options')
group.add_argument('-t', '--taxid-tag',
action="store",
dest="addtaxids:taxid_tag",
metavar="<TAXID_TAG>",
default=b"TAXID",
help="Name of the tag to store the found taxid "
"(default: 'TAXID'.")
group.add_argument('-n', '--taxon-name-tag',
action="store",
dest="addtaxids:taxon_name_tag",
metavar="<SCIENTIFIC_NAME_TAG>",
default=b"SCIENTIFIC_NAME",
help="Name of the tag giving the scientific name of the taxon "
"(default: 'SCIENTIFIC_NAME'.")
group.add_argument('-g', '--try-genus-match',
action="store_true", dest="addtaxids:try_genus_match",
default=False,
help="Try matching the first word of <SCIENTIFIC_NAME_TAG> when can't find corresponding taxid for a taxon. "
"If there is a match it is added in the 'parent_taxid' tag. (Can be used by 'obi taxonomy' to add the taxon under that taxid).")
group.add_argument('-a', '--restricting-ancestor',
action="store",
dest="addtaxids:restricting_ancestor",
metavar="<RESTRICTING_ANCESTOR>",
default=None,
help="Enables to restrict the search of taxids under an ancestor specified by its taxid.")
group.add_argument('-l', '--log-file',
action="store",
dest="addtaxids:log_file",
metavar="<LOG_FILE>",
default='',
help="Path to a log file to write informations about not found taxids.")
def run(config):
DMS.obi_atexit()
logger("info", "obi addtaxids")
# Open the input
input = open_uri(config['obi']['inputURI'])
if input is None:
raise Exception("Could not read input view")
i_dms = input[0]
i_view = input[1]
i_view_name = input[1].name
# Open the output: only the DMS, as the output view is going to be created by cloning the input view
# (could eventually be done via an open_uri() argument)
output = open_uri(config['obi']['outputURI'],
input=False,
dms_only=True)
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
output_0 = output[0]
o_view_name = output[1]
# stdout output: create temporary view
if type(output_0)==BufferedWriter:
o_dms = i_dms
i=0
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in output DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
imported_view_name = o_view_name
# If the input and output DMS are not the same, import the input view in the output DMS before cloning it to modify it
# (could be the other way around: clone and modify in the input DMS then import the new view in the output DMS)
if i_dms != o_dms:
imported_view_name = i_view_name
i=0
while imported_view_name in o_dms: # Making sure view name is unique in output DMS
imported_view_name = i_view_name+b"_"+str2bytes(str(i))
i+=1
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], i_view_name, imported_view_name)
i_view = o_dms[imported_view_name]
# Clone output view from input view
o_view = i_view.clone(o_view_name)
if o_view is None:
raise Exception("Couldn't create output view")
i_view.close()
# Open taxonomy
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]
# Initialize the progress bar
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(o_view), config)
else:
pb = None
try:
if config['addtaxids']['log_file']:
logfile = open(config['addtaxids']['log_file'], 'w')
else:
logfile = None
if config['addtaxids']['try_genus_match']:
try_genus = True
else:
try_genus = False
if 'restricting_ancestor' in config['addtaxids']:
res_anc = int(config['addtaxids']['restricting_ancestor'])
else:
res_anc = None
taxid_column_name = config['addtaxids']['taxid_tag']
parent_taxid_column_name = "PARENT_TAXID" # TODO macro
taxon_name_column_name = config['addtaxids']['taxon_name_tag']
taxid_column = Column.new_column(o_view, taxid_column_name, OBI_INT)
parent_taxid_column = Column.new_column(o_view, parent_taxid_column_name, OBI_INT)
taxon_name_column = o_view[taxon_name_column_name]
found_count = 0
not_found_count = 0
parent_found_count = 0
for i in range(len(o_view)):
PyErr_CheckSignals()
if pb is not None:
pb(i)
taxon_name = taxon_name_column[i]
taxon = taxo.get_taxon_by_name(taxon_name, res_anc)
if taxon is not None:
taxid_column[i] = taxon.taxid
found_count+=1
elif try_genus: # try finding genus or other parent taxon from the first word
taxon_name_sp = taxon_name.split(b" ")
taxon = taxo.get_taxon_by_name(taxon_name_sp[0], res_anc)
if taxon is not None:
parent_taxid_column[i] = taxon.taxid
parent_found_count+=1
if logfile:
print("Found parent taxon for", tostr(taxon_name), file=logfile)
else:
not_found_count+=1
if logfile:
print("No taxid found for", tostr(taxon_name), file=logfile)
else:
not_found_count+=1
if logfile:
print("No taxid found for", tostr(taxon_name), file=logfile)
except Exception, e:
raise RollbackException("obi addtaxids error, rollbacking view: "+str(e), o_view)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
logger("info", "\nTaxids found: "+str(found_count)+"/"+str(len(o_view))+" ("+str(round(found_count*100.0/len(o_view), 2))+"%)")
if config['addtaxids']['try_genus_match']:
logger("info", "\nParent taxids found: "+str(parent_found_count)+"/"+str(len(o_view))+" ("+str(round(parent_found_count*100.0/len(o_view), 2))+"%)")
logger("info", "\nTaxids not found: "+str(not_found_count)+"/"+str(len(o_view))+" ("+str(round(not_found_count*100.0/len(o_view), 2))+"%)")
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
input_dms_name=[input[0].name]
input_view_name=[i_view_name]
if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
o_view.write_config(config, "addtaxids", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(o_dms, imported_view_name)
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,7 +4,7 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
@ -12,6 +12,9 @@ from obitools3.utils cimport tobytes, str2bytes
from obitools3.dms.capi.obilcsalign cimport obi_lcs_align_one_column, \
obi_lcs_align_two_columns
from io import BufferedWriter
from cpython.exc cimport PyErr_CheckSignals
import time
import sys
@ -23,6 +26,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi align specific options')
@ -201,20 +205,20 @@ def run(config):
if output is None:
raise Exception("Could not create output")
o_dms = output[0]
output_0 = output[0]
o_dms_name = o_dms.name
final_o_view_name = output[1]
o_view_name = final_o_view_name
# If the input and output DMS are not the same, align creating a temporary view in the input dms that will be exported to
# If stdout output or the input and output DMS are not the same, align creating a temporary view in the input dms that will be exported to
# the right DMS and deleted in the other afterwards.
if i_dms != o_dms:
temporary_view_name = final_o_view_name
i=0
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = final_o_view_name+b"_"+str2bytes(str(i))
if i_dms != o_dms or type(output_0)==BufferedWriter:
if type(output_0)==BufferedWriter:
o_dms = i_dms
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in input DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
else:
o_view_name = final_o_view_name
# Save command config in View comments
command_line = " ".join(sys.argv[1:])
@ -263,12 +267,19 @@ def run(config):
View.delete_view(i_dms, i_view_name_2)
i_dms_2.close()
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view = o_dms[o_view_name]
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# If the input and the output DMS are different, delete the temporary result view in the input DMS
if i_dms != o_dms:
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
o_dms.close(force=True)
i_dms.close()
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -6,7 +6,7 @@ from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from obitools3.dms.column.column cimport Column
from obitools3.dms.capi.obiview cimport QUALITY_COLUMN
from obitools3.dms.capi.obitypes cimport OBI_QUAL
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.libalign._qsassemble import QSolexaReverseAssemble
@ -14,8 +14,10 @@ from obitools3.libalign._qsrassemble import QSolexaRightReverseAssemble
from obitools3.libalign._solexapairend import buildConsensus, buildJoinedSequence
from obitools3.dms.obiseq cimport Nuc_Seq
from obitools3.libalign.shifted_ali cimport Kmer_similarity, Ali_shifted
from obitools3.commands.ngsfilter import REVERSE_SEQ_COLUMN_NAME, REVERSE_QUALITY_COLUMN_NAME
from obitools3.dms.capi.obiview cimport REVERSE_SEQUENCE_COLUMN, REVERSE_QUALITY_COLUMN
from obitools3.utils cimport str2bytes
from io import BufferedWriter
import sys
import os
@ -29,6 +31,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi alignpairedend specific options')
@ -39,12 +42,13 @@ def addOptions(parser):
type=str,
help="URI to the reverse reads if they are in a different view than the forward reads")
group.add_argument('--score-min',
action="store", dest="alignpairedend:smin",
metavar="#.###",
default=None,
type=float,
help="Minimum score for keeping alignments")
# group.add_argument('--score-min',
# action="store", dest="alignpairedend:smin",
# metavar="#.###",
# default=None,
# type=float,
# help="Minimum score for keeping alignments. "
# "(for kmer alignment) The score is an approximation of the number of nucleotides matching in the overlap of the alignment.")
# group.add_argument('-A', '--true-ali',
# action="store_true", dest="alignpairedend:trueali",
@ -102,7 +106,7 @@ def alignmentIterator(entries, aligner):
seqR = reverse[i]
else:
seqF = Nuc_Seq.new_from_stored(entries[i])
seqR = Nuc_Seq(seqF.id, seqF[REVERSE_SEQ_COLUMN_NAME], quality=seqF[REVERSE_QUALITY_COLUMN_NAME])
seqR = Nuc_Seq(seqF.id, seqF[REVERSE_SEQUENCE_COLUMN], quality=seqF[REVERSE_QUALITY_COLUMN])
seqR.index = i
ali = aligner(seqF, seqR)
@ -170,18 +174,29 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
view = output[1]
output_0 = output[0]
o_dms = output[0]
Column.new_column(view, QUALITY_COLUMN, OBI_QUAL) #TODO output URI quality option?
if 'smin' in config['alignpairedend']:
smin = config['alignpairedend']['smin']
# stdout output: create temporary view
if type(output_0)==BufferedWriter:
i_dms = forward.dms # using any dms
o_dms = i_dms
i=0
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in input DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
o_view = View_NUC_SEQS.new(o_dms, o_view_name, quality=True)
else:
smin = 0
o_view = output[1]
Column.new_column(o_view, QUALITY_COLUMN, OBI_QUAL)
# Initialize the progress bar
pb = ProgressBar(entries_len, config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(entries_len, config)
else:
pb = None
#if config['alignpairedend']['trueali']:
# kmer_ali = False
# aligner = buildAlignment
@ -196,8 +211,8 @@ def run(config):
reversed_column=None)
else:
aligner = Kmer_similarity(entries, \
column2=entries[REVERSE_SEQ_COLUMN_NAME], \
qual_column2=entries[REVERSE_QUALITY_COLUMN_NAME], \
column2=entries[REVERSE_SEQUENCE_COLUMN], \
qual_column2=entries[REVERSE_QUALITY_COLUMN], \
kmer_size=config['alignpairedend']['kmersize'], \
reversed_column=entries[b'reversed']) # column created by the ngsfilter tool
@ -206,51 +221,63 @@ def run(config):
i = 0
for ali in ba:
pb(i)
if pb is not None:
pb(i)
PyErr_CheckSignals()
consensus = view[i]
consensus = o_view[i]
if not two_views:
seqF = entries[i]
else:
seqF = forward[i]
if ali.score > smin and ali.consensus_len > 0 :
if ali.overlap_len > 0 :
buildConsensus(ali, consensus, seqF)
else:
if not two_views:
seqR = Nuc_Seq(seqF.id, seqF[REVERSE_SEQ_COLUMN_NAME], quality = seqF[REVERSE_QUALITY_COLUMN_NAME])
seqR = Nuc_Seq(seqF.id, seqF[REVERSE_SEQUENCE_COLUMN], quality = seqF[REVERSE_QUALITY_COLUMN])
else:
seqR = reverse[i]
buildJoinedSequence(ali, seqR, consensus, forward=seqF)
consensus[b"smin"] = smin
if kmer_ali :
ali.free()
i+=1
pb(i, force=True)
print("", file=sys.stderr)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
if kmer_ali :
aligner.free()
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
view.write_config(config, "alignpairedend", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
output[0].record_command_line(command_line)
o_view.write_config(config, "alignpairedend", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(view), file=sys.stderr)
input[0].close()
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# If stdout output, delete the temporary imported view used to create the final file
if type(output_0)==BufferedWriter:
View_NUC_SEQS.delete_view(o_dms, o_view_name)
output_0.close()
# Close all DMS
input[0].close(force=True)
if two_views:
rinput[0].close()
output[0].close()
rinput[0].close(force=True)
o_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,16 +4,18 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from functools import reduce
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
from io import BufferedWriter
from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
ID_COLUMN, \
DEFINITION_COLUMN, \
QUALITY_COLUMN, \
COUNT_COLUMN
COUNT_COLUMN, \
TAXID_COLUMN
import time
import math
@ -33,6 +35,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi annotate specific options')
@ -175,8 +178,8 @@ def sequenceTaggerGenerator(config, taxo=None):
counter[0]+=1
for rank in annoteRank:
if 'taxid' in seq:
taxid = seq['taxid']
if TAXID_COLUMN in seq:
taxid = seq[TAXID_COLUMN]
if taxid is not None:
rtaxid = taxo.get_taxon_at_rank(taxid, rank)
if rtaxid is not None:
@ -190,58 +193,50 @@ def sequenceTaggerGenerator(config, taxo=None):
seq['seq_rank']=counter[0]
for i,v in toSet:
#try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(v, environ, seq)
#except Exception,e: # TODO discuss usefulness of this
# if options.onlyValid:
# raise e
# val = v
try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(v, environ, seq)
except Exception: # set string if not a valid expression
val = v
seq[i]=val
if length:
seq['seq_length']=len(seq)
if newId is not None:
# try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newId, environ, seq)
# except Exception,e:
# if options.onlyValid:
# raise e
# val = newId
try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newId, environ, seq)
except Exception: # set string if not a valid expression
val = newId
seq.id=val
if newDef is not None:
# try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newDef, environ, seq)
# except Exception,e:
# if options.onlyValid:
# raise e
# val = newDef
try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newDef, environ, seq)
except Exception: # set string if not a valid expression
val = newDef
seq.definition=val
#
if newSeq is not None:
# try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newSeq, environ, seq)
# except Exception,e:
# if options.onlyValid:
# raise e
# val = newSeq
try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
val = eval(newSeq, environ, seq)
except Exception: # set string if not a valid expression
val = newSeq
seq.seq=val
if 'seq_length' in seq:
seq['seq_length']=len(seq)
@ -251,15 +246,14 @@ def sequenceTaggerGenerator(config, taxo=None):
seq.view.delete_column(QUALITY_COLUMN)
if run is not None:
# try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
eval(run, environ, seq)
# except Exception,e:
# if options.onlyValid:
# raise e
try:
if taxo is not None:
environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
else:
environ = {'sequence':seq, 'counter':counter[0], 'math':math}
eval(run, environ, seq)
except Exception,e:
raise e
return sequenceTagger
@ -286,8 +280,19 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
output_0 = output[0]
o_view_name = output[1]
# stdout output: create temporary view
if type(output_0)==BufferedWriter:
o_dms = i_dms
i=0
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in output DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
imported_view_name = o_view_name
# If the input and output DMS are not the same, import the input view in the output DMS before cloning it to modify it
# (could be the other way around: clone and modify in the input DMS then import the new view in the output DMS)
if i_dms != o_dms:
@ -298,7 +303,7 @@ def run(config):
i+=1
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], i_view_name, imported_view_name)
i_view = o_dms[imported_view_name]
# Clone output view from input view
o_view = i_view.clone(o_view_name)
if o_view is None:
@ -315,7 +320,10 @@ def run(config):
taxo = None
# Initialize the progress bar
pb = ProgressBar(len(o_view), config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(o_view), config)
else:
pb = None
try:
@ -354,14 +362,16 @@ def run(config):
sequenceTagger = sequenceTaggerGenerator(config, taxo=taxo)
for i in range(len(o_view)):
PyErr_CheckSignals()
pb(i)
if pb is not None:
pb(i)
sequenceTagger(o_view[i])
except Exception, e:
raise RollbackException("obi annotate error, rollbacking view: "+str(e), o_view)
pb(i, force=True)
print("", file=sys.stderr)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
@ -371,15 +381,21 @@ def run(config):
input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
o_view.write_config(config, "annotate", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
output[0].record_command_line(command_line)
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary imported view used to create the final view
if i_dms != o_dms:
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(o_dms, imported_view_name)
o_dms.close()
i_dms.close()
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,14 +4,16 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms.dms cimport DMS
from obitools3.dms.view import RollbackException
from obitools3.dms.capi.build_reference_db cimport build_reference_db
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
from obitools3.dms.view.view cimport View
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from io import BufferedWriter
import sys
from cpython.exc cimport PyErr_CheckSignals
__title__="Tag a set of sequences for PCR and sequencing errors identification"
@ -22,6 +24,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi build_ref_db specific options')
@ -56,17 +59,20 @@ def run(config):
if output is None:
raise Exception("Could not create output")
o_dms = output[0]
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, build the database creating a temporary view that will be exported to
# If stdout output or the input and output DMS are not the same, build the database creating a temporary view that will be exported to
# the right DMS and deleted in the other afterwards.
if i_dms != o_dms:
temporary_view_name = final_o_view_name
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = final_o_view_name+b"_"+str2bytes(str(i))
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
@ -80,26 +86,33 @@ def run(config):
input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
comments = View.print_config(config, "build_ref_db", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
if build_reference_db(i_dms.name_with_full_path, tobytes(i_view_name), tobytes(taxonomy_name), tobytes(o_view_name), comments, config['build_ref_db']['threshold']) < 0:
raise Exception("Error building a reference database")
# If the input and output DMS are not the same, export result view to output DMS
if i_dms != o_dms:
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view = o_dms[o_view_name]
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# Save command config in DMS comments
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_dms[final_o_view_name]), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary result view in the input DMS
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
o_dms.close(force=True)
i_dms.close()
i_dms.close(force=True)
logger("info", "Done.")

164
python/obitools3/commands/cat.pyx Executable file
View File

@ -0,0 +1,164 @@
#cython: language_level=3
from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport str2bytes
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from obitools3.dms.view.view cimport View
from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, REVERSE_SEQUENCE_COLUMN, \
QUALITY_COLUMN, REVERSE_QUALITY_COLUMN
from obitools3.dms.capi.obitypes cimport OBI_SEQ, OBI_QUAL
from obitools3.dms.column.column cimport Column
from io import BufferedWriter
import time
import sys
from cpython.exc cimport PyErr_CheckSignals
__title__="Concatenate views."
def addOptions(parser):
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi cat specific options')
group.add_argument("-c",
action="append", dest="cat:views_to_cat",
metavar="<VIEW_NAME>",
default=[],
type=str,
help="URI of a view to concatenate. (e.g. 'my_dms/my_view'). "
"Several -c options can be used on the same "
"command line.")
def run(config):
DMS.obi_atexit()
logger("info", "obi cat")
# Check the views to concatenate
idms_list = []
iview_list = []
total_len = 0
remove_qual = False
remove_rev_qual = False
v_type = View_NUC_SEQS
for v_uri in config["cat"]["views_to_cat"]:
input = open_uri(v_uri)
if input is None:
raise Exception("Could not read input view")
i_dms = input[0]
i_view = input[1]
if input[2] != View_NUC_SEQS: # Check view type (output view is nuc_seqs view if all input view are nuc_seqs view)
v_type = View
if QUALITY_COLUMN not in i_view: # Check if keep quality column in output view (if all input views have it)
remove_qual = True
if REVERSE_QUALITY_COLUMN not in i_view: # same as above for reverse quality
remove_rev_qual = True
total_len += len(i_view)
idms_list.append(i_dms)
iview_list.append(i_view.name)
i_view.close()
# Open the output: only the DMS
output = open_uri(config['obi']['outputURI'],
input=False,
newviewtype=v_type)
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
output_0 = output[0]
o_view = output[1]
# stdout output
if type(output_0)==BufferedWriter:
o_dms = i_dms
# Initialize quality columns and their associated sequence columns if needed
if type(output_0) != BufferedWriter:
if not remove_qual:
if NUC_SEQUENCE_COLUMN not in o_view:
Column.new_column(o_view, NUC_SEQUENCE_COLUMN, OBI_SEQ)
Column.new_column(o_view, QUALITY_COLUMN, OBI_QUAL, associated_column_name=NUC_SEQUENCE_COLUMN, associated_column_version=o_view[NUC_SEQUENCE_COLUMN].version)
if not remove_rev_qual:
Column.new_column(o_view, REVERSE_SEQUENCE_COLUMN, OBI_SEQ)
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:
dict_cols = {}
for v_uri in config["cat"]["views_to_cat"]:
v = open_uri(v_uri)[1]
for coln in v.keys():
col = v[coln]
if v[coln].nb_elements_per_line > 1:
if coln not in dict_cols:
dict_cols[coln] = {}
dict_cols[coln]['eltnames'] = set(v[coln].elements_names)
dict_cols[coln]['nbelts'] = v[coln].nb_elements_per_line
dict_cols[coln]['obitype'] = v[coln].data_type_int
else:
dict_cols[coln]['eltnames'] = set(v[coln].elements_names + list(dict_cols[coln]['eltnames']))
dict_cols[coln]['nbelts'] = len(dict_cols[coln]['eltnames'])
v.close()
for coln in dict_cols:
Column.new_column(o_view, coln, dict_cols[coln]['obitype'],
nb_elements_per_line=dict_cols[coln]['nbelts'], elements_names=list(dict_cols[coln]['eltnames']))
# Initialize the progress bar
if not config['obi']['noprogressbar']:
pb = ProgressBar(total_len, config)
else:
pb = None
i = 0
for v_uri in config["cat"]["views_to_cat"]:
v = open_uri(v_uri)[1]
for entry in v:
PyErr_CheckSignals()
if pb is not None:
pb(i)
if type(output_0)==BufferedWriter:
rep = repr(entry)
output_0.write(str2bytes(rep)+b"\n")
else:
o_view[i] = entry
i+=1
v.close()
# Deletes quality columns if needed
if type(output_0)!=BufferedWriter:
if QUALITY_COLUMN in o_view and remove_qual :
o_view.delete_column(QUALITY_COLUMN)
if REVERSE_QUALITY_COLUMN in o_view and remove_rev_qual :
o_view.delete_column(REVERSE_QUALITY_COLUMN)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
# Save command config in DMS comments
command_line = " ".join(sys.argv[1:])
o_view.write_config(config, "cat", command_line, input_dms_name=[d.name for d in idms_list], input_view_name=[vname for vname in iview_list])
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(view), file=sys.stderr)
for d in idms_list:
d.close(force=True)
o_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,13 +4,14 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms.dms cimport DMS
from obitools3.dms.view import RollbackException
from obitools3.dms.capi.obiclean cimport obi_clean
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
from obitools3.dms.view.view cimport View
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from io import BufferedWriter
import sys
@ -21,7 +22,8 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi clean specific options')
group.add_argument('--distance', '-d',
@ -36,8 +38,7 @@ def addOptions(parser):
dest="clean:sample-tag-name",
metavar="<SAMPLE TAG NAME>",
type=str,
default="merged_sample",
help="Name of the tag where sample counts are kept.")
help="Name of the tag where merged sample count informations are kept (typically generated by obi uniq, usually MERGED_sample, default: None).")
group.add_argument('--ratio', '-r',
action="store", dest="clean:ratio",
@ -89,17 +90,20 @@ def run(config):
if output is None:
raise Exception("Could not create output")
o_dms = output[0]
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, run obiclean creating a temporary view that will be exported to
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported to
# the right DMS and deleted in the other afterwards.
if i_dms != o_dms:
temporary_view_name = final_o_view_name
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = final_o_view_name+b"_"+str2bytes(str(i))
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
@ -107,6 +111,9 @@ def run(config):
command_line = " ".join(sys.argv[1:])
comments = View.print_config(config, "clean", command_line, input_dms_name=[i_dms_name], input_view_name=[i_view_name])
if 'sample-tag-name' not in config['clean']:
config['clean']['sample-tag-name'] = ""
if obi_clean(i_dms.name_with_full_path, tobytes(i_view_name), tobytes(config['clean']['sample-tag-name']), tobytes(o_view_name), comments, \
config['clean']['distance'], config['clean']['ratio'], config['clean']['heads-only'], config['clean']['thread-count']) < 0:
raise Exception("Error running obiclean")
@ -114,18 +121,25 @@ def run(config):
# If the input and output DMS are not the same, export result view to output DMS
if i_dms != o_dms:
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view = o_dms[o_view_name]
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
# Save command config in DMS comments
o_dms.record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_dms[final_o_view_name]), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary result view in the input DMS
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
o_dms.close(force=True)
i_dms.close()
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -21,7 +21,10 @@ def run(config):
logger("info", "obi clean_dms")
if obi_clean_dms(tobytes(config['obi']['inputURI'])) < 0 :
dms_path = tobytes(config['obi']['inputURI'])
if b'.obidms' in dms_path:
dms_path = dms_path.split(b'.obidms')[0]
if obi_clean_dms(dms_path) < 0 :
raise Exception("Error cleaning DMS", config['obi']['inputURI'])
logger("info", "Done.")

View File

@ -22,7 +22,7 @@ def addOptions(parser):
group.add_argument('-s','--sequence',
action="store_true", dest="count:sequence",
default=False,
help="Prints only the number of sequence records.")
help="Prints only the number of sequence records (much faster, default: False).")
group.add_argument('-a','--all',
action="store_true", dest="count:all",
@ -56,3 +56,5 @@ def run(config):
print(count2)
else:
print(count1)
input[0].close(force=True)

View File

@ -5,10 +5,10 @@ from obitools3.dms.dms cimport DMS
from obitools3.dms.capi.obidms cimport OBIDMS_p
from obitools3.dms.view import RollbackException
from obitools3.dms.capi.obiecopcr cimport obi_ecopcr
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addTaxonomyOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addTaxonomyOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes
from obitools3.utils cimport tobytes, str2bytes
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from obitools3.dms.view import View
@ -16,6 +16,7 @@ from libc.stdlib cimport malloc, free
from libc.stdint cimport int32_t
import sys
from io import BufferedWriter
__title__="in silico PCR"
@ -27,6 +28,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi ecopcr specific options')
@ -35,13 +37,15 @@ def addOptions(parser):
action="store", dest="ecopcr:primer1",
metavar='<PRIMER>',
type=str,
help="Forward primer.")
required=True,
help="Forward primer, length must be less than or equal to 32")
group.add_argument('--primer2', '-R',
action="store", dest="ecopcr:primer2",
metavar='<PRIMER>',
type=str,
help="Reverse primer.")
required=True,
help="Reverse primer, length must be less than or equal to 32")
group.add_argument('--error', '-e',
action="store", dest="ecopcr:error",
@ -107,14 +111,20 @@ def addOptions(parser):
help="Defines the method used for estimating the Tm (melting temperature) between the primers and their corresponding "
"target sequences. SANTALUCIA: 1, or OWCZARZY: 2. Default: 1.")
group.add_argument('--keep-primers', '-p',
action="store_true",
dest="ecopcr:keep-primers",
default=False,
help="Whether to keep the primers attached to the output sequences (default: the primers are cut out).")
group.add_argument('--keep-nucs', '-D',
action="store",
dest="ecopcr:keep-nucs",
metavar="<INTEGER>",
metavar="<N>",
type=int,
default=0,
help="Keeps the specified number of nucleotides on each side of the in silico amplified sequences, "
"(already including the amplified DNA fragment plus the two target sequences of the primers).")
help="Keeps N nucleotides on each side of the in silico amplified sequences, "
"not including the primers (implying that primers are automatically kept if N > 0).")
group.add_argument('--kingdom-mode', '-k',
action="store_true",
@ -161,11 +171,20 @@ def run(config):
if output is None:
raise Exception("Could not create output")
o_dms = output[0]
output_0 = output[0]
o_dms_name = output[0].name
o_view_name = output[1]
# Read taxonomy name
taxonomy_name = config['obi']['taxoURI'].split("/")[-1] # Robust in theory
# If stdout output create a temporary view in the input dms that will be deleted afterwards.
if type(output_0)==BufferedWriter:
o_dms = i_dms
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in input DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
# Save command config in View comments
command_line = " ".join(sys.argv[1:])
@ -185,7 +204,7 @@ def run(config):
config['ecopcr']['min-length'], config['ecopcr']['max-length'], \
restrict_to_taxids_p, ignore_taxids_p, \
config['ecopcr']['circular'], config['ecopcr']['salt-concentration'], config['ecopcr']['salt-correction-method'], \
config['ecopcr']['keep-nucs'], config['ecopcr']['kingdom-mode']) < 0:
config['ecopcr']['keep-nucs'], config['ecopcr']['keep-primers'], config['ecopcr']['kingdom-mode']) < 0:
raise Exception("Error running ecopcr")
# Save command config in DMS comments
@ -193,10 +212,22 @@ def run(config):
free(restrict_to_taxids_p)
free(ignore_taxids_p)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view = o_dms[o_view_name]
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_dms[o_view_name]), file=sys.stderr)
o_dms.close()
# If stdout output, delete the temporary result view in the input DMS
if type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
i_dms.close(force=True)
o_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,7 +4,7 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms.dms cimport DMS
from obitools3.dms.view import RollbackException
from obitools3.dms.capi.obiecotag cimport obi_ecotag
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
@ -12,6 +12,7 @@ from obitools3.dms.view.view cimport View
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
import sys
from io import BufferedWriter
__title__="Taxonomic assignment of sequences"
@ -22,6 +23,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi ecotag specific options')
@ -39,6 +41,17 @@ def addOptions(parser):
help="Minimum identity to consider for assignment, as a normalized identity, e.g. 0.95 for an identity of 95%%. "
"Default: 0.00 (no threshold).")
group.add_argument('--minimum-circle','-c',
action="store", dest="ecotag:bubble_threshold",
metavar='<CIRCLE_THRESHOLD>',
default=0.99,
type=float,
help="Minimum identity considered for the assignment circle "
"(sequence is assigned to the LCA of all sequences within a similarity circle of the best matches; "
"the threshold for this circle is the highest value between <CIRCLE_THRESHOLD> and the best assignment score found for the query sequence). "
"Give value as a normalized identity, e.g. 0.95 for an identity of 95%%. "
"Default: 0.99.")
def run(config):
DMS.obi_atexit()
@ -64,10 +77,9 @@ def run(config):
ref_view_name = ref[1]
# Check that the threshold demanded is greater than or equal to the threshold used to build the reference database
if config['ecotag']['threshold'] < eval(i_dms[ref_view_name].comments["ref_db_threshold"]) :
print("Error: The threshold demanded (%f) is lower than the threshold used to build the reference database (%f).",
config['ecotag']['threshold'], i_dms[ref_view_name].comments["ref_db_threshold"])
if config['ecotag']['bubble_threshold'] < eval(ref_dms[ref_view_name].comments["ref_db_threshold"]) :
raise Exception(f"Error: The threshold demanded ({config['ecotag']['bubble_threshold']}) is lower than the threshold used to build the reference database ({float(ref_dms[ref_view_name].comments['ref_db_threshold'])}).")
# Open the output: only the DMS
output = open_uri(config['obi']['outputURI'],
input=False,
@ -75,17 +87,19 @@ def run(config):
if output is None:
raise Exception("Could not create output")
o_dms = output[0]
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, run ecotag creating a temporary view that will be exported to
# the right DMS and deleted in the other afterwards.
if i_dms != o_dms:
temporary_view_name = final_o_view_name
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported and deleted.
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = final_o_view_name+b"_"+str2bytes(str(i))
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
@ -107,10 +121,11 @@ def run(config):
comments = View.print_config(config, "ecotag", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
if obi_ecotag(i_dms.name_with_full_path, tobytes(i_view_name), \
tobytes(ref_dms_name), tobytes(ref_view_name), \
tobytes(taxo_dms_name), tobytes(taxonomy_name), \
tobytes(o_view_name), comments,
config['ecotag']['threshold']) < 0:
ref_dms.name_with_full_path, tobytes(ref_view_name), \
taxo_dms.name_with_full_path, tobytes(taxonomy_name), \
tobytes(o_view_name), comments, \
config['ecotag']['threshold'], \
config['ecotag']['bubble_threshold']) < 0:
raise Exception("Error running ecotag")
# If the input and output DMS are not the same, export result view to output DMS
@ -120,15 +135,24 @@ def run(config):
# Save command config in DMS comments
o_dms.record_command_line(command_line)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view = o_dms[o_view_name]
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_dms[final_o_view_name]), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary result view in the input DMS
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
o_dms.close(force=True)
i_dms.close()
taxo_dms.close(force=True)
ref_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -59,13 +59,23 @@ def run(config):
# Check that the input view has the type NUC_SEQS if needed # TODO discuss, maybe bool property
if (output[2] == Nuc_Seq) and (iview.type != b"NUC_SEQS_VIEW") : # Nuc_Seq_Stored? TODO
raise Exception("Error: the view to export in fasta or fastq format is not a NUC_SEQS view")
if config['obi']['only'] is not None:
withoutskip = min(input[4], config['obi']['only'])
else:
withoutskip = input[4]
if config['obi']['skip'] is not None:
skip = min(input[4], config['obi']['skip'])
else:
skip = 0
# Initialize the progress bar
if config['obi']['noprogressbar']:
pb = None
else:
pb = ProgressBar(len(iview), config, seconde=5)
pb = ProgressBar(withoutskip - skip, config)
i=0
for seq in iview :
PyErr_CheckSignals()
@ -79,14 +89,14 @@ def run(config):
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
print("", file=sys.stderr)
# TODO save command in input dms?
if not BrokenPipeError and not IOError:
output_object.close()
iview.close()
input[0].close()
input[0].close(force=True)
logger("info", "Done.")

View File

@ -4,7 +4,7 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, str2bytes
@ -14,6 +14,7 @@ import time
import re
import sys
import ast
from io import BufferedWriter
from cpython.exc cimport PyErr_CheckSignals
@ -28,6 +29,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group("obi grep specific options")
@ -36,14 +38,13 @@ def addOptions(parser):
metavar="<PREDICATE>",
default=[],
type=str,
help="Warning: use bytes for character strings (b'text' instead of 'text'). "
"Python boolean expression to be evaluated in the "
help="Python boolean expression to be evaluated in the "
"sequence/line context. The attribute name can be "
"used in the expression as a variable name. "
"An extra variable named 'sequence' or 'line' refers "
"to the sequence or line object itself. "
"Several -p options can be used on the same "
"commande line.")
"command line.")
group.add_argument("-S", "--sequence",
action="store", dest="grep:seq_pattern",
@ -162,8 +163,7 @@ def obi_eval(compiled_expr, loc_env, line):
return obi_eval_result
def Filter_generator(options, tax_filter):
#taxfilter = taxonomyFilterGenerator(options)
def Filter_generator(options, tax_filter, i_view):
# Initialize conditions
predicates = None
@ -172,6 +172,9 @@ def Filter_generator(options, tax_filter):
attributes = None
if "attributes" in options and len(options["attributes"]) > 0:
attributes = options["attributes"]
for attribute in attributes:
if attribute not in i_view:
return None
lmax = None
if "lmax" in options:
lmax = options["lmax"]
@ -197,6 +200,8 @@ def Filter_generator(options, tax_filter):
if "attribute_patterns" in options and len(options["attribute_patterns"]) > 0:
for p in options["attribute_patterns"]:
attribute, pattern = p.split(":", 1)
if attribute not in i_view:
return None
attribute_patterns[tobytes(attribute)] = re.compile(tobytes(pattern))
def filter(line, loc_env):
@ -301,16 +306,21 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
o_view_name_final = output[1]
o_view_name = o_view_name_final
# If the input and output DMS are not the same, create output view in input DMS first, then export it
# to output DMS, making sure the temporary view name is unique in the input DMS
if i_dms != o_dms:
output_0 = output[0]
final_o_view_name = output[1]
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported and deleted afterwards.
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while o_view_name in i_dms:
o_view_name = o_view_name_final+b"_"+str2bytes(str(i))
i+=1
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
taxo_uri = open_uri(config["obi"]["taxoURI"])
@ -321,33 +331,49 @@ def run(config):
taxo = None
# Initialize the progress bar
pb = ProgressBar(len(i_view), config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(i_view), config)
else:
pb = None
# Apply filter
tax_filter = Taxonomy_filter_generator(taxo, config["grep"])
filter = Filter_generator(config["grep"], tax_filter)
filter = Filter_generator(config["grep"], tax_filter, i_view)
selection = Line_selection(i_view)
for i in range(len(i_view)):
PyErr_CheckSignals()
pb(i)
line = i_view[i]
loc_env = {"sequence": line, "line": line, "taxonomy": taxo, "obi_eval_result": False}
good = filter(line, loc_env)
if good :
if filter is None and config["grep"]["invert_selection"]: # all sequences are selected: filter is None if no line will be selected because some columns don't exist
for i in range(len(i_view)):
PyErr_CheckSignals()
if pb is not None:
pb(i)
selection.append(i)
pb(i, force=True)
print("", file=sys.stderr)
elif filter is not None : # filter is None if no line will be selected because some columns don't exist
for i in range(len(i_view)):
PyErr_CheckSignals()
if pb is not None:
pb(i)
line = i_view[i]
loc_env = {"sequence": line, "line": line, "taxonomy": taxo, "obi_eval_result": False}
good = filter(line, loc_env)
if good :
selection.append(i)
if pb is not None:
pb(len(i_view), force=True)
print("", file=sys.stderr)
# Create output view with the line selection
try:
o_view = selection.materialize(o_view_name)
except Exception, e:
raise RollbackException("obi grep error, rollbacking view: "+str(e), o_view)
logger("info", "Grepped %d entries" % len(o_view))
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
input_dms_name=[input[0].name]
@ -362,16 +388,22 @@ def run(config):
# and delete the temporary view in the input DMS
if i_dms != o_dms:
o_view.close()
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, o_view_name_final)
o_view = o_dms[o_view_name_final]
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
o_view = o_dms[final_o_view_name]
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary imported view used to create the final view
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
i_dms.close()
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -4,14 +4,15 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport str2bytes
import time
import sys
from io import BufferedWriter
from cpython.exc cimport PyErr_CheckSignals
@ -22,6 +23,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi head specific options')
@ -53,31 +55,41 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
o_view_name_final = output[1]
o_view_name = o_view_name_final
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, create output view in input DMS first, then export it
# to output DMS, making sure the temporary view name is unique in the input DMS
if i_dms != o_dms:
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported and deleted.
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while o_view_name in i_dms:
o_view_name = o_view_name_final+b"_"+str2bytes(str(i))
i+=1
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
n = min(config['head']['count'], len(i_view))
# Initialize the progress bar
pb = ProgressBar(n, config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(n, config)
else:
pb = None
selection = Line_selection(i_view)
for i in range(n):
PyErr_CheckSignals()
pb(i)
if pb is not None:
pb(i)
selection.append(i)
pb(i, force=True)
print("", file=sys.stderr)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
# Create output view with the line selection
try:
@ -94,16 +106,22 @@ def run(config):
# and delete the temporary view in the input DMS
if i_dms != o_dms:
o_view.close()
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, o_view_name_final)
o_view = o_dms[o_view_name_final]
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
o_view = o_dms[final_o_view_name]
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(view), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary imported view used to create the final view
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
i_dms.close()
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -54,4 +54,5 @@ def run(config):
print(bytes2str(entries.ascii_history))
else:
raise Exception("ASCII history only available for views")
input[0].close(force=True)

View File

@ -11,6 +11,7 @@ from obitools3.dms.column.column cimport Column
from obitools3.dms.obiseq cimport Nuc_Seq
from obitools3.dms import DMS
from obitools3.dms.taxo.taxo cimport Taxonomy
from obitools3.files.uncompress cimport CompressedFile
from obitools3.utils cimport tobytes, \
@ -24,13 +25,16 @@ from obitools3.dms.capi.obiview cimport VIEW_TYPE_NUC_SEQS, \
DEFINITION_COLUMN, \
QUALITY_COLUMN, \
COUNT_COLUMN, \
TAXID_COLUMN
TAXID_COLUMN, \
MERGED_PREFIX, \
SCIENTIFIC_NAME_COLUMN
from obitools3.dms.capi.obidms cimport obi_import_view
from obitools3.dms.capi.obitypes cimport obitype_t, \
OBI_VOID, \
OBI_QUAL
OBI_QUAL, \
OBI_STR
from obitools3.dms.capi.obierrno cimport obi_errno
@ -45,6 +49,8 @@ from obitools3.apps.config import logger
from cpython.exc cimport PyErr_CheckSignals
from io import BufferedWriter
__title__="Imports sequences from different formats into a DMS"
@ -65,6 +71,19 @@ def addOptions(parser):
addTaxdumpInputOption(parser)
addMinimalOutputOption(parser)
group = parser.add_argument_group('obi import specific options')
group.add_argument('--preread',
action="store_true", dest="import:preread",
default=False,
help="Do a first readthrough of the dataset if it contains huge dictionaries (more than 100 keys) for "
"a much faster import. This option is not recommended and will slow down the import in any other case.")
group.add_argument('--space-priority',
action="store_true", dest="import:space_priority",
default=False,
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.")
def run(config):
@ -77,6 +96,7 @@ def run(config):
cdef obitype_t new_type
cdef bint get_quality
cdef bint NUC_SEQS_view
cdef bint silva
cdef int nb_elts
cdef object d
cdef View view
@ -87,6 +107,8 @@ def run(config):
cdef Column seq_col
cdef Column qual_col
cdef Column old_column
cdef Column sci_name_col
cdef bytes sci_name
cdef bint rewrite
cdef dict dcols
cdef int skipping
@ -120,7 +142,7 @@ def run(config):
if entry_count > 0:
logger("info", "Importing %d entries", entry_count)
else:
logger("info", "Importing an unknow number of entries")
logger("info", "Importing an unknown number of entries")
# TODO a bit dirty?
if input[2]==Nuc_Seq or input[2]==View_NUC_SEQS:
@ -130,7 +152,7 @@ def run(config):
else:
v = None
if config['obi']['taxdump'] or isinstance(input[1], View):
if config['obi']['taxdump'] or (isinstance(input[1], View) and not config['import']['space_priority']):
dms_only=True
else:
dms_only=False
@ -154,23 +176,24 @@ def run(config):
taxo.write(taxo_name)
taxo.close()
o_dms.record_command_line(" ".join(sys.argv[1:]))
o_dms.close()
o_dms.close(force=True)
logger("info", "Done.")
return
# If importing a view between two DMS, use C API
if isinstance(input[1], View):
# 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 :
input[0].close(force=True)
output[0].close(force=True)
raise Exception("Error importing a view in a DMS")
o_dms.record_command_line(" ".join(sys.argv[1:]))
o_dms.close()
input[0].close(force=True)
output[0].close(force=True)
logger("info", "Done.")
return
if entry_count >= 0:
pb = ProgressBar(entry_count, config, seconde=5)
entries = input[1]
pb = ProgressBar(entry_count, config)
NUC_SEQS_view = False
if isinstance(output[1], View) :
@ -185,9 +208,73 @@ def run(config):
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":
silva = True
sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
else:
silva = False
dcols = {}
# First read through the entries to prepare columns with dictionaries as they are very time-expensive to rewrite
if config['import']['preread']:
logger("info", "First readthrough...")
entries = input[1]
i = 0
dict_dict = {}
for entry in entries:
PyErr_CheckSignals()
if entry is None: # error or exception handled at lower level, not raised because Python generators can't resume after any exception is raised
if config['obi']['skiperror']:
i-=1
continue
else:
raise Exception("obi import error in first readthrough")
if pb is not None:
pb(i)
elif not i%50000:
logger("info", "Read %d entries", i)
for tag in entry :
newtag = tag
if tag[:7] == b"merged_":
newtag = MERGED_PREFIX+tag[7:]
if type(entry[tag]) == dict :
if tag in dict_dict:
dict_dict[newtag][0].update(entry[tag].keys())
else:
dict_dict[newtag] = [set(entry[tag].keys()), get_obitype(entry[tag])]
i+=1
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
for tag in dict_dict:
dcols[tag] = (Column.new_column(view, tag, dict_dict[tag][1], \
nb_elements_per_line=len(dict_dict[tag][0]), \
elements_names=list(dict_dict[tag][0])), \
dict_dict[tag][1])
# Reinitialize the input
if isinstance(input[0], CompressedFile):
input_is_file = True
if entry_count >= 0:
pb = ProgressBar(entry_count, config)
try:
input[0].close()
except AttributeError:
pass
input = open_uri(config['obi']['inputURI'], force_file=input_is_file)
if input is None:
raise Exception("Could not open input URI")
entries = input[1]
i = 0
for entry in entries :
@ -195,7 +282,6 @@ def run(config):
if entry is None: # error or exception handled at lower level, not raised because Python generators can't resume after any exception is raised
if config['obi']['skiperror']:
i-=1
continue
else:
raise RollbackException("obi import error, rollbacking view", view)
@ -204,115 +290,139 @@ def run(config):
pb(i)
elif not i%50000:
logger("info", "Imported %d entries", i)
if NUC_SEQS_view:
id_col[i] = entry.id
def_col[i] = entry.definition
seq_col[i] = entry.seq
# Check if there is a sequencing quality associated by checking the first entry # TODO haven't found a more robust solution yet
if i == 0:
get_quality = QUALITY_COLUMN in entry
try:
if NUC_SEQS_view:
id_col[i] = entry.id
def_col[i] = entry.definition
seq_col[i] = entry.seq
# Check if there is a sequencing quality associated by checking the first entry # TODO haven't found a more robust solution yet
if i == 0:
get_quality = QUALITY_COLUMN in entry
if get_quality:
Column.new_column(view, QUALITY_COLUMN, OBI_QUAL)
qual_col = view[QUALITY_COLUMN]
if get_quality:
Column.new_column(view, QUALITY_COLUMN, OBI_QUAL)
qual_col = view[QUALITY_COLUMN]
if get_quality:
qual_col[i] = entry.quality
for tag in entry :
if tag != ID_COLUMN and tag != DEFINITION_COLUMN and tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN : # TODO dirty
value = entry[tag]
if tag == b"taxid":
tag = TAXID_COLUMN
if tag == b"count":
tag = COUNT_COLUMN
qual_col[i] = entry.quality
if tag not in dcols :
value_type = type(value)
nb_elts = 1
value_obitype = OBI_VOID
if value_type == dict or value_type == list :
nb_elts = len(value)
elt_names = list(value)
else :
nb_elts = 1
elt_names = None
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), value_obitype)
# Fill value
dcols[tag][0][i] = value
# TODO else log error?
else :
rewrite = False
# Check type adequation
old_type = dcols[tag][1]
new_type = OBI_VOID
new_type = update_obitype(old_type, value)
if old_type != new_type :
rewrite = True
try:
# Fill value
dcols[tag][0][i] = value
except IndexError :
# Parse taxon scientific name if SILVA file
if silva:
sci_name = entry.definition.split(b";")[-1]
sci_name_col[i] = sci_name
for tag in entry :
if tag != ID_COLUMN and tag != DEFINITION_COLUMN and tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN : # TODO dirty
value = entry[tag]
if tag == b"taxid":
tag = TAXID_COLUMN
if tag == b"count":
tag = COUNT_COLUMN
if tag[:7] == b"merged_":
tag = MERGED_PREFIX+tag[7:]
if tag not in dcols :
value_type = type(value)
old_column = dcols[tag][0]
old_nb_elements_per_line = old_column.nb_elements_per_line
new_nb_elements_per_line = 0
old_elements_names = old_column.elements_names
new_elements_names = None
nb_elts = 1
value_obitype = OBI_VOID
if value_type == dict or value_type == list :
nb_elts = len(value)
elt_names = list(value)
else :
nb_elts = 1
elt_names = None
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), value_obitype)
# Fill value
if value_type == dict and nb_elts == 1: # special case that makes the OBI3 create a 1 elt/line column which won't read a dict value
value = value[list(value.keys())[0]] # The solution is to transform the value in a simple atomic one acceptable by the column
dcols[tag][0][i] = value
# TODO else log error?
#####################################################################
# Check the length and keys of column lines if needed
if value_type == dict : # Check dictionary keys
for k in value :
if k not in old_elements_names :
new_elements_names = list(set(old_elements_names+[tobytes(k) for k in value]))
rewrite = True
break
elif value_type == list or value_type == tuple : # Check vector length
if old_nb_elements_per_line < len(value) :
new_nb_elements_per_line = len(value)
rewrite = True
#####################################################################
if rewrite :
if new_nb_elements_per_line == 0 and new_elements_names is not None :
new_nb_elements_per_line = len(new_elements_names)
# Reset obierrno
obi_errno = 0
dcols[tag] = (view.rewrite_column_with_diff_attributes(old_column.name,
new_data_type=new_type,
new_nb_elements_per_line=new_nb_elements_per_line,
new_elements_names=new_elements_names,
rewrite_last_line=False),
new_type)
# Update the dictionary:
for t in dcols :
dcols[t] = (view[t], dcols[t][1])
else :
rewrite = False
# Check type adequation
old_type = dcols[tag][1]
new_type = OBI_VOID
new_type = update_obitype(old_type, value)
if old_type != new_type :
rewrite = True
try:
# Check that it's not the case where the first entry contained a dict of length 1 and now there is a new key
if type(value) == dict and \
dcols[tag][0].nb_elements_per_line == 1 \
and set(dcols[tag][0].elements_names) != set(value.keys()) :
raise IndexError # trigger column rewrite
# Fill value
dcols[tag][0][i] = value
except IndexError :
value_type = type(value)
old_column = dcols[tag][0]
old_nb_elements_per_line = old_column.nb_elements_per_line
new_nb_elements_per_line = 0
old_elements_names = old_column.elements_names
new_elements_names = None
#####################################################################
# Check the length and keys of column lines if needed
if value_type == dict : # Check dictionary keys
for k in value :
if k not in old_elements_names :
new_elements_names = list(set(old_elements_names+[tobytes(k) for k in value]))
rewrite = True
break
elif value_type == list or value_type == tuple : # Check vector length
if old_nb_elements_per_line < len(value) :
new_nb_elements_per_line = len(value)
rewrite = True
#####################################################################
if rewrite :
if new_nb_elements_per_line == 0 and new_elements_names is not None :
new_nb_elements_per_line = len(new_elements_names)
# Reset obierrno
obi_errno = 0
dcols[tag] = (view.rewrite_column_with_diff_attributes(old_column.name,
new_data_type=new_type,
new_nb_elements_per_line=new_nb_elements_per_line,
new_elements_names=new_elements_names,
rewrite_last_line=False),
new_type)
# Update the dictionary:
for t in dcols :
dcols[t] = (view[t], dcols[t][1])
# Fill value
dcols[tag][0][i] = value
except Exception as e:
print("\nCould not import sequence id:", entry.id, "(error raised:", e, ")")
if 'skiperror' in config['obi'] and not config['obi']['skiperror']:
raise e
else:
pass
i+=1
if pb is not None:
@ -333,7 +443,7 @@ def run(config):
except AttributeError:
pass
try:
output[0].close()
output[0].close(force=True)
except AttributeError:
pass

View File

@ -46,5 +46,5 @@ def run(config):
process.wait()
iview.close()
input[0].close()
input[0].close(force=True)

View File

@ -31,24 +31,11 @@ def run(config):
input = open_uri(config['obi']['inputURI'])
if input is None:
raise Exception("Could not read input")
if input[2] == DMS and not config['ls']['longformat']:
dms = input[0]
l = []
for view in input[0]:
l.append(tostr(view) + "\t(Date created: " + str(bytes2str_object(dms[view].comments["Date created"]))+")")
l.sort()
for v in l:
print(v)
# Print representation
if config['ls']['longformat']:
print(input[1].repr_longformat())
else:
print(repr(input[1]))
if input[2] == DMS:
taxolist = ["\n### Taxonomies:"]
for t in Taxonomy.list_taxos(input[0]):
taxolist.append("\t"+tostr(t))
if len(taxolist) > 1:
for t in taxolist:
print(t)
if config['ls']['longformat'] and len(input[1].comments) > 0:
print("\n### Comments:")
print(str(input[1].comments))
input[0].close(force=True)

357
python/obitools3/commands/ngsfilter.pyx Executable file → Normal file
View File

@ -2,10 +2,10 @@
from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view import RollbackException
from obitools3.dms.view import RollbackException, View
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from obitools3.dms.column.column cimport Column, Column_line
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.libalign._freeendgapfm import FreeEndGapFullMatch
@ -13,17 +13,19 @@ from obitools3.libalign.apat_pattern import Primer_search
from obitools3.dms.obiseq cimport Nuc_Seq
from obitools3.dms.capi.obitypes cimport OBI_SEQ, OBI_QUAL
from obitools3.dms.capi.apat cimport MAX_PATTERN
from obitools3.utils cimport tobytes
from obitools3.dms.capi.obiview cimport REVERSE_SEQUENCE_COLUMN, REVERSE_QUALITY_COLUMN
from obitools3.utils cimport tobytes, str2bytes
from libc.stdint cimport INT32_MAX
from functools import reduce
import math
import sys
from cpython.exc cimport PyErr_CheckSignals
from io import BufferedWriter
REVERSE_SEQ_COLUMN_NAME = b"REVERSE_SEQUENCE" # used by alignpairedend tool
REVERSE_QUALITY_COLUMN_NAME = b"REVERSE_QUALITY" # used by alignpairedend tool
#REVERSE_SEQ_COLUMN_NAME = b"REVERSE_SEQUENCE" # used by alignpairedend tool
#REVERSE_QUALITY_COLUMN_NAME = b"REVERSE_QUALITY" # used by alignpairedend tool
__title__="Assigns sequence records to the corresponding experiment/sample based on DNA tags and primers"
@ -33,7 +35,8 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi ngsfilter specific options')
group.add_argument('-t','--info-view',
@ -41,7 +44,9 @@ def addOptions(parser):
metavar="<URI>",
type=str,
default=None,
help="URI to the view containing the samples definition (with tags, primers, sample names,...)")
required=True,
help="URI to the view containing the samples definition (with tags, primers, sample names,...).\n"
"\nWarning: primer lengths must be less than or equal to 32")
group.add_argument('-R', '--reverse-reads',
action="store", dest="ngsfilter:reverse",
@ -55,7 +60,12 @@ def addOptions(parser):
metavar="<URI>",
type=str,
default=None,
help="URI to the view used to store the sequences unassigned to any sample")
help="URI to the view used to store the sequences unassigned to any sample. Those sequences are untrimmed.")
group.add_argument('--no-tags',
action="store_true", dest="ngsfilter:notags",
default=False,
help="Use this option if your experiment does not use tags to identify samples")
group.add_argument('-e','--error',
action="store", dest="ngsfilter:error",
@ -166,8 +176,15 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
primer_list = []
i=0
for p in info_view:
# Check primer length: should not be longer than 32, the max allowed by the apat lib
if len(p[b'forward_primer']) > 32:
raise RollbackException("Error: primers can not be longer than 32bp, rollbacking views")
if len(p[b'reverse_primer']) > 32:
raise RollbackException("Error: primers can not be longer than 32bp, rollbacking views")
forward=Primer(p[b'forward_primer'],
len(p[b'forward_tag']) if p[b'forward_tag']!=b'-' else None,
len(p[b'forward_tag']) if (b'forward_tag' in p and p[b'forward_tag']!=None) else None,
True,
max_errors=max_errors,
verbose=verbose,
@ -178,7 +195,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
infos[forward]=fp
reverse=Primer(p[b'reverse_primer'],
len(p[b'reverse_tag']) if p[b'reverse_tag']!=b'-' else None,
len(p[b'reverse_tag']) if (b'reverse_tag' in p and p[b'reverse_tag']!=None) else None,
False,
max_errors=max_errors,
verbose=verbose,
@ -213,10 +230,11 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
rpp=rp.get(cf,{})
rp[cf]=rpp
tags = (p[b'forward_tag'] if p[b'forward_tag']!=b'-' else None,
p[b'reverse_tag'] if p[b'reverse_tag']!=b'-' else None)
tags = (p[b'forward_tag'] if (b'forward_tag' in p and p[b'forward_tag']!=None) else None,
p[b'reverse_tag'] if (b'reverse_tag' in p and p[b'reverse_tag']!=None) else None)
assert tags not in dpp, \
if tags != (None, None):
assert tags not in dpp, \
"Tag pair %s is already used with primer pairs: (%s,%s)" % (str(tags),forward,reverse)
# Save additional data
@ -234,7 +252,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
return infos, primer_list
cdef tuple annotate(sequences, infos, verbose=False):
cdef tuple annotate(sequences, infos, no_tags, verbose=False):
def sortMatch(match):
if match[1] is None:
@ -249,17 +267,12 @@ cdef tuple annotate(sequences, infos, verbose=False):
return match[1][1]
not_aligned = len(sequences) > 1
sequenceF = sequences[0]
sequenceR = None
if not not_aligned:
final_sequence = sequenceF
else:
final_sequence = sequenceF.clone() # TODO maybe not cloning and then deleting quality tags is more efficient
sequences[0] = sequences[0].clone()
if not_aligned:
sequenceR = sequences[1]
final_sequence[REVERSE_SEQ_COLUMN_NAME] = sequenceR.seq # used by alignpairedend tool
final_sequence[REVERSE_QUALITY_COLUMN_NAME] = sequenceR.quality # used by alignpairedend tool
sequences[1] = sequences[1].clone()
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
for seq in sequences:
if hasattr(seq, "quality_array"):
@ -275,8 +288,6 @@ cdef tuple annotate(sequences, infos, verbose=False):
# Try direct matching:
directmatch = []
first_matched_seq = None
second_matched_seq = None
for seq in sequences:
new_seq = True
pattern = 0
@ -295,60 +306,96 @@ cdef tuple annotate(sequences, infos, verbose=False):
directmatch = directmatch[0] if directmatch[0][1] is not None else None
if directmatch is None:
final_sequence[b'error']=b'No primer match'
return False, final_sequence
if not_aligned:
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
sequences[0][b'error']=b'No primer match'
return False, sequences[0]
first_matched_seq = directmatch[2]
if id(first_matched_seq) == id(sequenceF) and not_aligned:
second_matched_seq = sequenceR
if id(directmatch[2]) == id(sequences[0]):
first_match_first_seq = True
else:
second_matched_seq = sequenceF
match = first_matched_seq[directmatch[1][1]:directmatch[1][2]]
first_match_first_seq = False
match = directmatch[2][directmatch[1][1]:directmatch[1][2]]
if not not_aligned:
final_sequence[b'seq_length_ori']=len(final_sequence)
sequences[0][b'seq_length_ori']=len(sequences[0])
if not not_aligned or id(first_matched_seq) == id(sequenceF):
final_sequence = final_sequence[directmatch[1][2]:]
if not not_aligned or first_match_first_seq:
sequences[0] = sequences[0][directmatch[1][2]:]
else:
cut_seq = sequenceR[directmatch[1][2]:]
final_sequence[REVERSE_SEQ_COLUMN_NAME] = cut_seq.seq # used by alignpairedend tool
final_sequence[REVERSE_QUALITY_COLUMN_NAME] = cut_seq.quality # used by alignpairedend tool
sequences[1] = sequences[1][directmatch[1][2]:]
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
if directmatch[0].forward:
final_sequence[b'direction']=b'forward'
final_sequence[b'forward_errors']=directmatch[1][0]
final_sequence[b'forward_primer']=directmatch[0].raw
final_sequence[b'forward_match']=match.seq
sequences[0][b'direction']=b'forward'
sequences[0][b'forward_errors']=directmatch[1][0]
sequences[0][b'forward_primer']=directmatch[0].raw
sequences[0][b'forward_match']=match.seq
else:
final_sequence[b'direction']=b'reverse'
final_sequence[b'reverse_errors']=directmatch[1][0]
final_sequence[b'reverse_primer']=directmatch[0].raw
final_sequence[b'reverse_match']=match.seq
sequences[0][b'direction']=b'reverse'
sequences[0][b'reverse_errors']=directmatch[1][0]
sequences[0][b'reverse_primer']=directmatch[0].raw
sequences[0][b'reverse_match']=match.seq
# Keep only paired reverse primer
infos = infos[directmatch[0]]
infos = infos[directmatch[0]]
reverse_primer = list(infos.keys())[0]
direct_primer = directmatch[0]
# If not aligned, look for other match in already computed matches (choose the one that makes the biggest amplicon)
if not_aligned:
i=1
# TODO comment
while i<len(all_direct_matches) and (all_direct_matches[i][1] is None or all_direct_matches[i][0].forward == directmatch[0].forward or all_direct_matches[i][0] == directmatch[0]):
while i<len(all_direct_matches) and \
(all_direct_matches[i][1] is None or \
all_direct_matches[i][0].forward == directmatch[0].forward or \
all_direct_matches[i][0] == directmatch[0] or \
reverse_primer != all_direct_matches[i][0]) :
i+=1
if i < len(all_direct_matches):
reversematch = all_direct_matches[i]
else:
reversematch = None
# Cut reverse primer out of 1st matched seq if it contains it, because if it's also in the other sequence, the next step will "choose" only the one on the other sequence
if not_aligned:
# do it on same seq
if first_match_first_seq:
r = reverse_primer.revcomp(sequences[0])
else:
r = reverse_primer.revcomp(sequences[1])
if r is not None: # found
if first_match_first_seq :
sequences[0] = sequences[0][:r[1]]
else:
sequences[1] = sequences[1][:r[1]]
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
# do the same on the other seq
if first_match_first_seq:
r = direct_primer.revcomp(sequences[1])
else:
r = direct_primer.revcomp(sequences[0])
if r is not None: # found
if first_match_first_seq:
sequences[1] = sequences[1][:r[1]]
else:
sequences[0] = sequences[0][:r[1]]
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality
# Look for other primer in the other direction on the sequence, or
# If sequences are not already aligned and reverse primer not found in most likely sequence (the one without the forward primer), try matching on the same sequence than the first match (primer in the other direction)
if not not_aligned or (not_aligned and (reversematch is None or reversematch[1] is None)):
if not not_aligned:
sequence_to_match = second_matched_seq
if not_aligned and first_match_first_seq:
seq_to_match = sequences[1]
else:
sequence_to_match = first_matched_seq
seq_to_match = sequences[0]
reversematch = []
# Compute begin
begin=directmatch[1][2]+1 # end of match + 1 on the same sequence
@ -365,7 +412,7 @@ cdef tuple annotate(sequences, infos, verbose=False):
primer=p
# Saving original primer as 4th member of the tuple to serve as correct key in infos dict even if it might have been reversed complemented
# (3rd member already used by directmatch)
reversematch.append((primer, primer(sequence_to_match, same_sequence=not new_seq, pattern=pattern, begin=begin), None, p))
reversematch.append((primer, primer(seq_to_match, same_sequence=not new_seq, pattern=pattern, begin=begin), None, p))
new_seq = False
pattern+=1
# Choose match closer to the end of the sequence
@ -378,11 +425,11 @@ cdef tuple annotate(sequences, infos, verbose=False):
message = b'No reverse primer match'
else:
message = b'No direct primer match'
final_sequence[b'error']=message
return False, final_sequence
sequences[0][b'error']=message
return False, sequences[0]
if reversematch is None:
final_sequence[b'status']=b'partial'
sequences[0][b'status']=b'partial'
if directmatch[0].forward:
tags=(directmatch[1][3],None)
@ -392,78 +439,86 @@ cdef tuple annotate(sequences, infos, verbose=False):
samples = infos[None]
else:
final_sequence[b'status']=b'full'
sequences[0][b'status']=b'full'
if not not_aligned or first_match_first_seq:
match = sequences[0][reversematch[1][1]:reversematch[1][2]]
else:
match = sequences[1][reversematch[1][1]:reversematch[1][2]]
match = second_matched_seq[reversematch[1][1]:reversematch[1][2]]
match = match.reverse_complement
if not not_aligned or id(second_matched_seq) == id(sequenceF):
final_sequence = final_sequence[0:reversematch[1][1]]
else:
cut_seq = sequenceR[reversematch[1][2]:]
if not not_aligned:
sequences[0] = sequences[0][0:reversematch[1][1]]
elif first_match_first_seq:
sequences[1] = sequences[1][reversematch[1][2]:]
if not directmatch[0].forward:
cut_seq = cut_seq.reverse_complement
final_sequence[REVERSE_SEQ_COLUMN_NAME] = cut_seq.seq # used by alignpairedend tool
final_sequence[REVERSE_QUALITY_COLUMN_NAME] = cut_seq.quality # used by alignpairedend tool
sequences[1] = sequences[1].reverse_complement
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]:]
if directmatch[0].forward:
tags=(directmatch[1][3], reversematch[1][3])
final_sequence[b'reverse_errors'] = reversematch[1][0]
final_sequence[b'reverse_primer'] = reversematch[0].raw
final_sequence[b'reverse_match'] = match.seq
sequences[0][b'reverse_errors'] = reversematch[1][0]
sequences[0][b'reverse_primer'] = reversematch[0].raw
sequences[0][b'reverse_match'] = match.seq
else:
tags=(reversematch[1][3], directmatch[1][3])
final_sequence[b'forward_errors'] = reversematch[1][0]
final_sequence[b'forward_primer'] = reversematch[0].raw
final_sequence[b'forward_match'] = match.seq
sequences[0][b'forward_errors'] = reversematch[1][0]
sequences[0][b'forward_primer'] = reversematch[0].raw
sequences[0][b'forward_match'] = match.seq
if tags[0] is not None:
final_sequence[b'forward_tag'] = tags[0]
sequences[0][b'forward_tag'] = tags[0]
if tags[1] is not None:
final_sequence[b'reverse_tag'] = tags[1]
sequences[0][b'reverse_tag'] = tags[1]
samples = infos[reversematch[3]]
if not directmatch[0].forward:
final_sequence = final_sequence.reverse_complement
final_sequence[b'reversed'] = True # used by the alignpairedend tool (in kmer_similarity.c)
sequences[0] = sequences[0].reverse_complement
sequences[0][b'reversed'] = True # used by the alignpairedend tool (in kmer_similarity.c)
else:
sequences[0][b'reversed'] = False # used by the alignpairedend tool (in kmer_similarity.c)
sample=None
if tags[0] is not None: # Direct tag known
if tags[1] is not None: # Reverse tag known
sample = samples.get(tags, None)
else: # Only direct tag known
s=[samples[x] for x in samples if x[0]==tags[0]]
if len(s)==1:
sample=s[0]
elif len(s)>1:
final_sequence[b'error']=b'multiple samples match tags'
return False, final_sequence
else:
sample=None
else:
if tags[1] is not None: # Only reverse tag known
s=[samples[x] for x in samples if x[1]==tags[1]]
if len(s)==1:
sample=s[0]
elif len(s)>1:
final_sequence[b'error']=b'multiple samples match tags'
return False, final_sequence
else:
sample=None
if not no_tags:
if tags[0] is not None: # Direct tag known
if tags[1] is not None: # Reverse tag known
sample = samples.get(tags, None)
else: # Only direct tag known
s=[samples[x] for x in samples if x[0]==tags[0]]
if len(s)==1:
sample=s[0]
elif len(s)>1:
sequences[0][b'error']=b'Did not found reverse tag'
return False, sequences[0]
else:
sample=None
else:
if tags[1] is not None: # Only reverse tag known
s=[samples[x] for x in samples if x[1]==tags[1]]
if len(s)==1:
sample=s[0]
elif len(s)>1:
sequences[0][b'error']=b'Did not found forward tag'
return False, sequences[0]
else:
sample=None
if sample is None:
sequences[0][b'error']=b"No sample with that tag combination"
return False, sequences[0]
if sample is None:
final_sequence[b'error']=b"Cannot assign sequence to a sample"
return False, final_sequence
final_sequence.update(sample)
sequences[0].update(sample)
if not not_aligned:
final_sequence[b'seq_length']=len(final_sequence)
sequences[0][b'seq_length']=len(sequences[0])
return True, final_sequence
return True, sequences[0]
def run(config):
@ -486,7 +541,8 @@ def run(config):
raise Exception("Could not open input reads")
if input[2] != View_NUC_SEQS:
raise NotImplementedError('obi ngsfilter only works on NUC_SEQS views')
i_dms = input[0]
if "reverse" in config["ngsfilter"]:
forward = input[1]
@ -527,8 +583,19 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
output_0 = output[0]
o_view = output[1]
# If stdout output, create a temporary view in the input dms that will be deleted afterwards.
if type(output_0)==BufferedWriter:
o_dms = i_dms
o_view_name = b"temp"
while o_view_name in i_dms: # Making sure view name is unique in input DMS
o_view_name = o_view_name+b"_"+str2bytes(str(i))
i+=1
o_view = View_NUC_SEQS.new(i_dms, o_view_name)
# Open the view containing the informations about the tags and the primers
info_input = open_uri(config['ngsfilter']['info_view'])
if info_input is None:
@ -549,10 +616,19 @@ def run(config):
unidentified = None
# Initialize the progress bar
pb = ProgressBar(entries_len, config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(entries_len, config)
else:
pb = None
# Check and store primers and tags
infos, primer_list = read_info_view(info_view, max_errors=config['ngsfilter']['error'], verbose=False, not_aligned=not_aligned) # TODO obi verbose option
try:
infos, primer_list = read_info_view(info_view, max_errors=config['ngsfilter']['error'], verbose=False, not_aligned=not_aligned) # TODO obi verbose option
except RollbackException, e:
if unidentified is not None:
raise RollbackException("obi ngsfilter error, rollbacking views: "+str(e), o_view, unidentified)
else:
raise RollbackException("obi ngsfilter error, rollbacking view: "+str(e), o_view)
aligner = Primer_search(primer_list, config['ngsfilter']['error'])
@ -564,50 +640,75 @@ def run(config):
paired_p.revcomp.aligner = aligner
if not_aligned: # create columns used by alignpairedend tool
Column.new_column(o_view, REVERSE_SEQ_COLUMN_NAME, OBI_SEQ)
Column.new_column(o_view, REVERSE_QUALITY_COLUMN_NAME, OBI_QUAL, associated_column_name=REVERSE_SEQ_COLUMN_NAME, associated_column_version=o_view[REVERSE_SEQ_COLUMN_NAME].version)
Column.new_column(o_view, REVERSE_SEQUENCE_COLUMN, OBI_SEQ)
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)
Column.new_column(unidentified, REVERSE_SEQ_COLUMN_NAME, OBI_SEQ)
Column.new_column(unidentified, REVERSE_QUALITY_COLUMN_NAME, OBI_QUAL, associated_column_name=REVERSE_SEQ_COLUMN_NAME, associated_column_version=unidentified[REVERSE_SEQ_COLUMN_NAME].version)
if unidentified is not None:
Column.new_column(unidentified, REVERSE_SEQUENCE_COLUMN, OBI_SEQ)
Column.new_column(unidentified, REVERSE_QUALITY_COLUMN, OBI_QUAL, associated_column_name=REVERSE_SEQUENCE_COLUMN, associated_column_version=unidentified[REVERSE_SEQUENCE_COLUMN].version)
g = 0
u = 0
no_tags = config['ngsfilter']['notags']
try:
for i in range(entries_len):
PyErr_CheckSignals()
pb(i)
if pb is not None:
pb(i)
if not_aligned:
modseq = [Nuc_Seq.new_from_stored(forward[i]), Nuc_Seq.new_from_stored(reverse[i])]
else:
modseq = [Nuc_Seq.new_from_stored(entries[i])]
good, oseq = annotate(modseq, infos)
good, oseq = annotate(modseq, infos, no_tags)
if good:
o_view[g].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq)
g+=1
elif unidentified is not None:
unidentified[u].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq)
# Untrim sequences (put original back)
if len(modseq) > 1:
oseq[REVERSE_SEQUENCE_COLUMN] = reverse[i].seq
oseq[REVERSE_QUALITY_COLUMN] = reverse[i].quality
unidentified[u].set(oseq.id, forward[i].seq, definition=oseq.definition, quality=forward[i].quality, tags=oseq)
else:
unidentified[u].set(oseq.id, entries[i].seq, definition=oseq.definition, quality=entries[i].quality, tags=oseq)
u+=1
except Exception, e:
raise RollbackException("obi ngsfilter error, rollbacking views: "+str(e), o_view, unidentified)
if unidentified is not None:
raise RollbackException("obi ngsfilter error, rollbacking views: "+str(e), o_view, unidentified)
else:
raise RollbackException("obi ngsfilter error, rollbacking view: "+str(e), o_view)
pb(i, force=True)
print("", file=sys.stderr)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
o_view.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
unidentified.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
# Add comment about unidentified seqs
unidentified.comments["info"] = "View containing sequences categorized as unidentified by the ngsfilter command"
output[0].record_command_line(command_line)
if unidentified is not None:
unidentified.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
# Add comment about unidentified seqs
unidentified.comments["info"] = "View containing sequences categorized as unidentified by the ngsfilter command"
o_dms.record_command_line(command_line)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
input[0].close()
output[0].close()
info_input[0].close()
unidentified_input[0].close()
# If stdout output, delete the temporary result view in the input DMS
if type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
i_dms.close(force=True)
o_dms.close(force=True)
info_input[0].close(force=True)
if unidentified is not None:
unidentified_input[0].close(force=True)
aligner.free()
logger("info", "Done.")

View File

@ -4,7 +4,7 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport str2bytes
@ -24,6 +24,7 @@ from obitools3.dms.capi.obitypes cimport OBI_BOOL, \
import time
import sys
from cpython.exc cimport PyErr_CheckSignals
from io import BufferedWriter
NULL_VALUE = {OBI_BOOL: OBIBool_NA,
@ -42,6 +43,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi sort specific options')
@ -59,7 +61,7 @@ def addOptions(parser):
def line_cmp(line, key, pb):
pb
pb
if line[key] is None:
return NULL_VALUE[line.view[key].data_type_int]
else:
@ -86,20 +88,28 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
o_view_name_final = output[1]
o_view_name = o_view_name_final
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, create output view in input DMS first, then export it
# to output DMS, making sure the temporary view name is unique in the input DMS
if i_dms != o_dms:
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported and deleted.
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while o_view_name in i_dms:
o_view_name = o_view_name_final+b"_"+str2bytes(str(i))
i+=1
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
# Initialize the progress bar
pb = ProgressBar(len(i_view), config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(i_view), config)
else:
pb = None
keys = config['sort']['keys']
selection = Line_selection(i_view)
@ -110,10 +120,14 @@ def run(config):
for k in keys: # TODO order?
PyErr_CheckSignals()
selection.sort(key=lambda line_idx: line_cmp(i_view[line_idx], k, pb(line_idx)), reverse=config['sort']['reverse'])
pb(len(i_view), force=True)
print("", file=sys.stderr)
if pb is not None:
selection.sort(key=lambda line_idx: line_cmp(i_view[line_idx], k, pb(line_idx)), reverse=config['sort']['reverse'])
else:
selection.sort(key=lambda line_idx: line_cmp(i_view[line_idx], k, None), reverse=config['sort']['reverse'])
if pb is not None:
pb(len(i_view), force=True)
print("", file=sys.stderr)
# Create output view with the sorted line selection
try:
@ -132,16 +146,23 @@ def run(config):
# and delete the temporary view in the input DMS
if i_dms != o_dms:
o_view.close()
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, o_view_name_final)
o_view = o_dms[o_view_name_final]
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
o_view = o_dms[final_o_view_name]
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary imported view used to create the final view
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
i_dms.close()
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -162,7 +162,7 @@ def run(config):
lcat=0
# Initialize the progress bar
pb = ProgressBar(len(i_view), config, seconde=5)
pb = ProgressBar(len(i_view), config)
for i in range(len(i_view)):
PyErr_CheckSignals()
@ -251,7 +251,7 @@ def run(config):
for i in range(len(sorted_stats)):
c = sorted_stats[i][0]
for v in c:
if v is not None:
if type(v) == bytes:
print(pcat % tostr(v)+"\t", end="")
else:
print(pcat % str(v)+"\t", end="")
@ -268,6 +268,6 @@ def run(config):
print("%7d" %catcount[c], end="")
print("%9d" %totcount[c])
input[0].close()
input[0].close(force=True)
logger("info", "Done.")

View File

@ -4,7 +4,7 @@ from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.dms import DMS
from obitools3.dms.view.view cimport View, Line_selection
from obitools3.uri.decode import open_uri
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption
from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputOption, addNoProgressBarOption
from obitools3.dms.view import RollbackException
from obitools3.apps.config import logger
from obitools3.utils cimport str2bytes
@ -12,6 +12,7 @@ from obitools3.utils cimport str2bytes
import time
import sys
from cpython.exc cimport PyErr_CheckSignals
from io import BufferedWriter
__title__="Keep the N last lines of a view."
@ -21,6 +22,7 @@ def addOptions(parser):
addMinimalInputOption(parser)
addMinimalOutputOption(parser)
addNoProgressBarOption(parser)
group=parser.add_argument_group('obi tail specific options')
@ -52,31 +54,41 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
o_dms = output[0]
o_view_name_final = output[1]
o_view_name = o_view_name_final
output_0 = output[0]
final_o_view_name = output[1]
# If the input and output DMS are not the same, create output view in input DMS first, then export it
# to output DMS, making sure the temporary view name is unique in the input DMS
if i_dms != o_dms:
# If stdout output or the input and output DMS are not the same, create a temporary view that will be exported and deleted.
if i_dms != o_dms or type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while o_view_name in i_dms:
o_view_name = o_view_name_final+b"_"+str2bytes(str(i))
i+=1
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
if type(output_0)==BufferedWriter:
o_dms = i_dms
else:
o_view_name = final_o_view_name
start = max(len(i_view) - config['tail']['count'], 0)
# Initialize the progress bar
pb = ProgressBar(len(i_view) - start, config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(i_view) - start, config)
else:
pb = None
selection = Line_selection(i_view)
for i in range(start, len(i_view)):
PyErr_CheckSignals()
pb(i)
if pb is not None:
pb(i)
selection.append(i)
pb(i, force=True)
print("", file=sys.stderr)
if pb is not None:
pb(i, force=True)
print("", file=sys.stderr)
# Save command config in View comments
command_line = " ".join(sys.argv[1:])
@ -97,16 +109,22 @@ def run(config):
# and delete the temporary view in the input DMS
if i_dms != o_dms:
o_view.close()
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, o_view_name_final)
o_view = o_dms[o_view_name_final]
View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], o_view_name, final_o_view_name)
o_view = o_dms[final_o_view_name]
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
# If the input and the output DMS are different, delete the temporary imported view used to create the final view
if i_dms != o_dms:
# If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
if i_dms != o_dms or type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
o_dms.close()
i_dms.close()
o_dms.close(force=True)
i_dms.close(force=True)
logger("info", "Done.")

View File

@ -23,6 +23,7 @@ from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
import shutil
import string
import random
import sys
from cpython.exc cimport PyErr_CheckSignals
@ -366,7 +367,7 @@ def random_new_view(config, infos, first=False):
infos['view'] = View_NUC_SEQS.new(infos['dms'], random_unique_name(infos), comments=random_comments(config)) # TODO quality column
else :
infos['view'] = View.new(infos['dms'], random_unique_name(infos), comments=random_comments(config)) # TODO quality column
infos['view'].write_config(config, "test", infos["command_line"], input_dms_name=[infos['dms'].name], input_view_name=["random"])
print_test(config, repr(infos['view']))
if v_to_clone is not None :
if line_selection is None:
@ -441,7 +442,7 @@ def addOptions(parser):
default=20,
type=int,
help="Maximum length of tuples. "
"Default: 200")
"Default: 50")
group.add_argument('--max_ini_col_count','-o',
action="store", dest="test:maxinicolcount",
@ -457,7 +458,7 @@ def addOptions(parser):
default=10000,
type=int,
help="Maximum number of lines in a column. "
"Default: 10000")
"Default: 1000")
group.add_argument('--max_elts_per_line','-e',
action="store", dest="test:maxelts",
@ -497,7 +498,8 @@ def run(config):
(b"OBI_SEQ", False): random_seq, (b"OBI_SEQ", True): random_seq_tuples,
(b"OBI_STR", False): random_bytes, (b"OBI_STR", True): random_bytes_tuples
},
'tests': [test_set_and_get, test_add_col, test_delete_col, test_col_alias, test_new_view]
'tests': [test_set_and_get, test_add_col, test_delete_col, test_col_alias, test_new_view],
'command_line': " ".join(sys.argv[1:])
}
# TODO ???
@ -529,7 +531,7 @@ def run(config):
test_taxo(config, infos)
infos['view'].close()
infos['dms'].close()
infos['dms'].close(force=True)
shutil.rmtree(config['obi']['defaultdms']+'.obidms', ignore_errors=True)
print("Done.")

View File

@ -5,5 +5,5 @@ from obitools3.dms.taxo.taxo cimport Taxonomy
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy)
cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, list mergedKeys_list=*, Taxonomy taxonomy=*, bint mergeIds=*, list categories=*, int max_elts=*)
cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy, dict config)
cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, dict config, list mergedKeys_list=*, Taxonomy taxonomy=*, bint mergeIds=*, list categories=*, int max_elts=*)

View File

@ -8,22 +8,24 @@ from obitools3.dms.view import RollbackException
from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
from obitools3.dms.column.column cimport Column, Column_line
from obitools3.dms.capi.obiview cimport QUALITY_COLUMN, COUNT_COLUMN, NUC_SEQUENCE_COLUMN, ID_COLUMN, TAXID_COLUMN, \
TAXID_DIST_COLUMN, MERGED_TAXID_COLUMN, MERGED_COLUMN, MERGED_PREFIX
TAXID_DIST_COLUMN, MERGED_TAXID_COLUMN, MERGED_COLUMN, MERGED_PREFIX, \
REVERSE_QUALITY_COLUMN
from obitools3.dms.capi.obitypes cimport OBI_INT, OBI_STR, index_t
from obitools3.apps.optiongroups import addMinimalInputOption, \
addMinimalOutputOption, \
addTaxonomyOption, \
addEltLimitOption
addEltLimitOption, \
addNoProgressBarOption
from obitools3.uri.decode import open_uri
from obitools3.apps.config import logger
from obitools3.utils cimport tobytes, tostr
from obitools3.utils cimport tobytes, tostr, str2bytes
import sys
from cpython.exc cimport PyErr_CheckSignals
from io import BufferedWriter
__title__="Group sequence records together"
def addOptions(parser):
@ -32,6 +34,7 @@ def addOptions(parser):
addTaxonomyOption(parser)
addMinimalOutputOption(parser)
addEltLimitOption(parser)
addNoProgressBarOption(parser)
group = parser.add_argument_group('obi uniq specific options')
@ -56,7 +59,7 @@ def addOptions(parser):
"(option can be used several times).")
cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy, dict config) :
cdef int taxid
cdef Nuc_Seq_Stored seq
@ -69,7 +72,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
cdef object gn_sn
cdef object fa_sn
# Create columns
# Create columns and save them for efficiency
if b"species" in o_view and o_view[b"species"].data_type_int != OBI_INT :
o_view.delete_column(b"species")
if b"species" not in o_view:
@ -77,6 +80,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"species",
OBI_INT
)
species_column = o_view[b"species"]
if b"genus" in o_view and o_view[b"genus"].data_type_int != OBI_INT :
o_view.delete_column(b"genus")
@ -85,6 +89,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"genus",
OBI_INT
)
genus_column = o_view[b"genus"]
if b"family" in o_view and o_view[b"family"].data_type_int != OBI_INT :
o_view.delete_column(b"family")
@ -93,6 +98,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"family",
OBI_INT
)
family_column = o_view[b"family"]
if b"species_name" in o_view and o_view[b"species_name"].data_type_int != OBI_STR :
o_view.delete_column(b"species_name")
@ -101,6 +107,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"species_name",
OBI_STR
)
species_name_column = o_view[b"species_name"]
if b"genus_name" in o_view and o_view[b"genus_name"].data_type_int != OBI_STR :
o_view.delete_column(b"genus_name")
@ -109,6 +116,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"genus_name",
OBI_STR
)
genus_name_column = o_view[b"genus_name"]
if b"family_name" in o_view and o_view[b"family_name"].data_type_int != OBI_STR :
o_view.delete_column(b"family_name")
@ -117,6 +125,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"family_name",
OBI_STR
)
family_name_column = o_view[b"family_name"]
if b"rank" in o_view and o_view[b"rank"].data_type_int != OBI_STR :
o_view.delete_column(b"rank")
@ -125,6 +134,7 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"rank",
OBI_STR
)
rank_column = o_view[b"rank"]
if b"scientific_name" in o_view and o_view[b"scientific_name"].data_type_int != OBI_STR :
o_view.delete_column(b"scientific_name")
@ -133,9 +143,19 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
b"scientific_name",
OBI_STR
)
for seq in o_view:
PyErr_CheckSignals()
scientific_name_column = o_view[b"scientific_name"]
# Initialize the progress bar
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(o_view), config)
else:
pb = None
i=0
for seq in o_view:
PyErr_CheckSignals()
if pb is not None:
pb(i)
if MERGED_TAXID_COLUMN in seq :
m_taxids = []
m_taxids_dict = seq[MERGED_TAXID_COLUMN]
@ -165,20 +185,24 @@ cdef merge_taxonomy_classification(View_NUC_SEQS o_view, Taxonomy taxonomy) :
else:
fa_sn = None
tfa = None
seq[b"species"] = tsp
seq[b"genus"] = tgn
seq[b"family"] = tfa
seq[b"species_name"] = sp_sn
seq[b"genus_name"] = gn_sn
seq[b"family_name"] = fa_sn
seq[b"rank"] = taxonomy.get_rank(taxid)
seq[b"scientific_name"] = taxonomy.get_scientific_name(taxid)
species_column[i] = tsp
genus_column[i] = tgn
family_column[i] = tfa
species_name_column[i] = sp_sn
genus_name_column[i] = gn_sn
family_name_column[i] = fa_sn
rank_column[i] = taxonomy.get_rank(taxid)
scientific_name_column[i] = taxonomy.get_scientific_name(taxid)
i+=1
if pb is not None:
pb(len(o_view), force=True)
cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, list mergedKeys_list=None, Taxonomy taxonomy=None, bint mergeIds=False, list categories=None, int max_elts=1000000) :
cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, dict config, list mergedKeys_list=None, Taxonomy taxonomy=None, bint mergeIds=False, list categories=None, int max_elts=1000000) :
cdef int i
cdef int k
@ -187,6 +211,7 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
cdef int u_idx
cdef int i_idx
cdef int i_count
cdef int o_count
cdef str key_str
cdef bytes key
cdef bytes mkey
@ -209,7 +234,6 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
cdef Nuc_Seq_Stored i_seq
cdef Nuc_Seq_Stored o_seq
cdef Nuc_Seq_Stored u_seq
cdef Column i_col
cdef Column i_seq_col
cdef Column i_id_col
cdef Column i_taxid_col
@ -217,6 +241,8 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
cdef Column o_id_col
cdef Column o_taxid_dist_col
cdef Column o_merged_col
cdef Column o_count_col
cdef Column i_count_col
cdef Column_line i_mcol
cdef object taxid_dist_dict
cdef object iter_view
@ -252,7 +278,12 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
mergedKeys_m = []
for k in range(k_count):
mergedKeys_m.append(MERGED_PREFIX + mergedKeys[k])
# Check that not trying to remerge without total count information
for key in mergedKeys_m:
if key in view and COUNT_COLUMN not in view:
raise Exception("\n>>>>\nError: trying to re-merge tags without total count tag. Run obi annotate to add the count tag from the relevant merged tag, i.e.: \nobi annotate --set-tag COUNT:'sum([value for key,value in sequence['MERGED_sample'].items()])' dms/input dms/output\n")
if categories is None:
categories = []
@ -274,7 +305,8 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
iter_view = iter(view)
for i_seq in iter_view :
PyErr_CheckSignals()
pb(i)
if pb is not None:
pb(i)
# This can't be done in the same line as the unique_id tuple creation because it generates a bug
# where Cython (version 0.25.2) does not detect the reference to the categs_list variable and deallocates
@ -284,6 +316,7 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
for x in categories :
catl.append(i_seq[x])
#unique_id = tuple(catl) + (i_seq_col[i],)
unique_id = tuple(catl) + (i_seq_col.get_line_idx(i),)
#unique_id = tuple(i_seq[x] for x in categories) + (seq_col.get_line_idx(i),) # The line that cython can't read properly
@ -320,7 +353,14 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
for k in range(k_count):
key = mergedKeys[k]
merged_col_name = mergedKeys_m[k]
i_col = view[key]
if merged_infos[merged_col_name]['nb_elts'] == 1:
raise Exception("Can't merge information from a tag with only one element (e.g. one sample ; don't use -m option)")
if merged_col_name in view:
i_col = view[merged_col_name]
else:
i_col = view[key]
if merged_infos[merged_col_name]['nb_elts'] > max_elts:
str_merged_cols.append(merged_col_name)
@ -374,23 +414,35 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
alias=MERGED_COLUMN
)
# Keep columns that are going to be used a lot in variables
# Keep columns in variables for efficiency
o_id_col = o_view[ID_COLUMN]
if TAXID_DIST_COLUMN in o_view:
o_taxid_dist_col = o_view[TAXID_DIST_COLUMN]
if MERGED_COLUMN in o_view:
o_merged_col = o_view[MERGED_COLUMN]
if COUNT_COLUMN not in o_view:
Column.new_column(o_view,
COUNT_COLUMN,
OBI_INT)
o_count_col = o_view[COUNT_COLUMN]
if COUNT_COLUMN in view:
i_count_col = view[COUNT_COLUMN]
if pb is not None:
pb(len(view), force=True)
print("")
pb(len(view), force=True)
print("")
logger("info", "Second browsing through the input")
# Initialize the progress bar
pb = ProgressBar(len(uniques), seconde=5)
if pb is not None:
pb = ProgressBar(len(view))
o_idx = 0
total_treated = 0
for unique_id in uniques :
PyErr_CheckSignals()
pb(o_idx)
merged_sequences = uniques[unique_id]
@ -407,7 +459,7 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
merged_list = list(set(merged_list)) # deduplicate the list
o_merged_col[o_idx] = merged_list
o_seq[COUNT_COLUMN] = 0
o_count = 0
if TAXID_DIST_COLUMN in u_seq and i_taxid_dist_col[u_idx] is not None:
taxid_dist_dict = i_taxid_dist_col[u_idx]
@ -419,16 +471,20 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
merged_dict[mkey] = {}
for i_idx in merged_sequences:
PyErr_CheckSignals()
if pb is not None:
pb(total_treated)
i_id = i_id_col[i_idx]
i_seq = view[i_idx]
if COUNT_COLUMN not in i_seq or i_seq[COUNT_COLUMN] is None:
if COUNT_COLUMN not in i_seq or i_count_col[i_idx] is None:
i_count = 1
else:
i_count = i_seq[COUNT_COLUMN]
i_count = i_count_col[i_idx]
o_seq[COUNT_COLUMN] += i_count
o_count += i_count
for k in range(k_count):
@ -463,42 +519,55 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
mcol[key2] = i_mcol[key2]
else:
mcol[key2] = mcol[key2] + i_mcol[key2]
# Write taxid_dist
if mergeIds and TAXID_COLUMN in mergedKeys:
if TAXID_DIST_COLUMN in str_merged_cols:
o_taxid_dist_col[o_idx] = str(taxid_dist_dict)
else:
o_taxid_dist_col[o_idx] = taxid_dist_dict
# Write merged dicts
for mkey in merged_dict:
if mkey in str_merged_cols:
mkey_cols[mkey][o_idx] = str(merged_dict[mkey])
else:
mkey_cols[mkey][o_idx] = merged_dict[mkey]
# Sets NA values to 0 # TODO discuss, for now keep as None and test for None instead of testing for 0 in tools
#for key in mkey_cols[mkey][o_idx]:
# if mkey_cols[mkey][o_idx][key] is None:
# mkey_cols[mkey][o_idx][key] = 0
for key in i_seq.keys():
# Delete informations that differ between the merged sequences
# TODO make special columns list?
# TODO make special columns list? // could be more efficient
if key != COUNT_COLUMN and key != ID_COLUMN and key != NUC_SEQUENCE_COLUMN and key in o_seq and o_seq[key] != i_seq[key] \
and key not in merged_dict :
o_seq[key] = None
total_treated += 1
# Write merged dicts
for mkey in merged_dict:
if mkey in str_merged_cols:
mkey_cols[mkey][o_idx] = str(merged_dict[mkey])
else:
mkey_cols[mkey][o_idx] = merged_dict[mkey]
# Sets NA values to 0 # TODO discuss, for now keep as None and test for None instead of testing for 0 in tools
#for key in mkey_cols[mkey][o_idx]:
# if mkey_cols[mkey][o_idx][key] is None:
# mkey_cols[mkey][o_idx][key] = 0
# Write taxid_dist
if mergeIds and TAXID_COLUMN in mergedKeys:
if TAXID_DIST_COLUMN in str_merged_cols:
o_taxid_dist_col[o_idx] = str(taxid_dist_dict)
else:
o_taxid_dist_col[o_idx] = taxid_dist_dict
o_count_col[o_idx] = o_count
o_idx += 1
# Deletes quality column if there is one because the matching between sequence and quality will be broken (quality set to NA when sequence not)
if pb is not None:
pb(len(view), force=True)
# Deletes quality columns if there is one because the matching between sequence and quality will be broken (quality set to NA when sequence not)
if QUALITY_COLUMN in view:
o_view.delete_column(QUALITY_COLUMN)
if REVERSE_QUALITY_COLUMN in view:
o_view.delete_column(REVERSE_QUALITY_COLUMN)
# Delete old columns that are now merged
for k in range(k_count):
if mergedKeys[k] in o_view:
o_view.delete_column(mergedKeys[k])
if taxonomy is not None:
print("") # TODO because in the middle of progress bar. Better solution?
logger("info", "Merging taxonomy classification")
merge_taxonomy_classification(o_view, taxonomy)
merge_taxonomy_classification(o_view, taxonomy, config)
@ -530,8 +599,23 @@ def run(config):
if output is None:
raise Exception("Could not create output view")
i_dms = input[0]
entries = input[1]
o_view = output[1]
o_dms = output[0]
output_0 = output[0]
# If stdout output create a temporary view that will be exported and deleted.
if type(output_0)==BufferedWriter:
temporary_view_name = b"temp"
i=0
while temporary_view_name in i_dms: # Making sure view name is unique in input DMS
temporary_view_name = temporary_view_name+b"_"+str2bytes(str(i))
i+=1
o_view_name = temporary_view_name
o_dms = i_dms
o_view = View_NUC_SEQS.new(i_dms, o_view_name)
else:
o_view = output[1]
if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
taxo_uri = open_uri(config['obi']['taxoURI'])
@ -542,15 +626,19 @@ def run(config):
taxo = None
# Initialize the progress bar
pb = ProgressBar(len(entries), config, seconde=5)
if config['obi']['noprogressbar'] == False:
pb = ProgressBar(len(entries), config)
else:
pb = None
try:
uniq_sequences(entries, o_view, pb, mergedKeys_list=config['uniq']['merge'], taxonomy=taxo, mergeIds=config['uniq']['mergeids'], categories=config['uniq']['categories'], max_elts=config['obi']['maxelts'])
except Exception, e:
raise RollbackException("obi uniq error, rollbacking view: "+str(e), o_view)
if len(entries) > 0:
try:
uniq_sequences(entries, o_view, pb, config, mergedKeys_list=config['uniq']['merge'], taxonomy=taxo, mergeIds=config['uniq']['mergeids'], categories=config['uniq']['categories'], max_elts=config['obi']['maxelts'])
except Exception, e:
raise RollbackException("obi uniq error, rollbacking view: "+str(e), o_view)
pb(len(entries), force=True)
print("", file=sys.stderr)
if pb is not None:
print("", file=sys.stderr)
# Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:])
@ -560,13 +648,23 @@ def run(config):
input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
o_view.write_config(config, "uniq", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
output[0].record_command_line(command_line)
o_dms.record_command_line(command_line)
# stdout output: write to buffer
if type(output_0)==BufferedWriter:
logger("info", "Printing to output...")
o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
o_view.close()
#print("\n\nOutput view:\n````````````", file=sys.stderr)
#print(repr(o_view), file=sys.stderr)
input[0].close()
output[0].close()
# If stdout output, delete the temporary result view in the input DMS
if type(output_0)==BufferedWriter:
View.delete_view(i_dms, o_view_name)
i_dms.close(force=True)
o_dms.close(force=True)
logger("info", "Done.")

View File

@ -34,6 +34,7 @@ cdef extern from "obidms.h" nogil:
int obi_close_dms(OBIDMS_p dms, bint force)
char* obi_dms_get_dms_path(OBIDMS_p dms)
char* obi_dms_get_full_path(OBIDMS_p dms, const_char_p path_name)
char* obi_dms_formatted_infos(OBIDMS_p dms, bint detailed)
void obi_close_atexit()
obiversion_t obi_import_column(const char* dms_path_1, const char* dms_path_2, const char* column_name, obiversion_t version_number)

View File

@ -68,3 +68,6 @@ cdef extern from "obidmscolumn.h" nogil:
int obi_column_write_comments(OBIDMS_column_p column, const char* comments)
int obi_column_add_comment(OBIDMS_column_p column, const char* key, const char* value)
char* obi_column_formatted_infos(OBIDMS_column_p column, bint detailed)

View File

@ -23,6 +23,7 @@ cdef extern from "obi_ecopcr.h" nogil:
double salt_concentration,
int salt_correction_method,
int keep_nucleotides,
bint keep_primers,
bint kingdom_mode)

View File

@ -11,4 +11,5 @@ cdef extern from "obi_ecotag.h" nogil:
const char* taxonomy_name,
const char* output_view_name,
const char* output_view_comments,
double ecotag_threshold)
double ecotag_threshold,
double bubble_threshold)

View File

@ -7,6 +7,8 @@ from libc.stdint cimport int32_t
cdef extern from "obidms_taxonomy.h" nogil:
extern int MIN_LOCAL_TAXID
struct ecotxnode :
int32_t taxid
int32_t rank
@ -18,6 +20,13 @@ cdef extern from "obidms_taxonomy.h" nogil:
ctypedef ecotxnode ecotx_t
struct econame_t : # can't get this struct to be accepted by Cython ('unknown size')
char* name
char* class_name
int32_t is_scientific_name
ecotxnode* taxon
struct ecotxidx_t :
int32_t count
int32_t max_taxid
@ -30,9 +39,14 @@ cdef extern from "obidms_taxonomy.h" nogil:
char** label
struct econameidx_t :
int32_t count
econame_t* names
struct OBIDMS_taxonomy_t :
ecorankidx_t* ranks
# econameidx_t* names
econameidx_t* names
ecotxidx_t* taxa
ctypedef OBIDMS_taxonomy_t* OBIDMS_taxonomy_p
@ -51,7 +65,11 @@ cdef extern from "obidms_taxonomy.h" nogil:
ecotx_t* obi_taxo_get_parent_at_rank(ecotx_t* taxon, int32_t rankidx)
ecotx_t* obi_taxo_get_taxon_with_taxid(OBIDMS_taxonomy_p taxonomy, int32_t taxid)
char* obi_taxo_get_name_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx)
ecotx_t* obi_taxo_get_taxon_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx)
bint obi_taxo_is_taxon_under_taxid(ecotx_t* taxon, int32_t other_taxid)
ecotx_t* obi_taxo_get_species(ecotx_t* taxon, OBIDMS_taxonomy_p taxonomy)
@ -71,4 +89,4 @@ cdef extern from "obidms_taxonomy.h" nogil:
int obi_taxo_add_preferred_name_with_taxon(OBIDMS_taxonomy_p tax, ecotx_t* taxon, const char* preferred_name)
const char* obi_taxo_rank_index_to_label(int32_t rank_idx, ecorankidx_t* ranks)

View File

@ -24,7 +24,10 @@ cdef extern from "obiview.h" nogil:
extern const_char_p ID_COLUMN
extern const_char_p DEFINITION_COLUMN
extern const_char_p QUALITY_COLUMN
extern const_char_p REVERSE_QUALITY_COLUMN
extern const_char_p REVERSE_SEQUENCE_COLUMN
extern const_char_p COUNT_COLUMN
extern const_char_p SCIENTIFIC_NAME_COLUMN
extern const_char_p TAXID_COLUMN
extern const_char_p MERGED_TAXID_COLUMN
extern const_char_p MERGED_PREFIX
@ -100,14 +103,18 @@ cdef extern from "obiview.h" nogil:
const_char_p comments,
bint create)
int obi_view_delete_column(Obiview_p view, const_char_p column_name)
int obi_view_delete_column(Obiview_p view, const_char_p column_name, bint delete_file)
OBIDMS_column_p obi_view_get_column(Obiview_p view, const_char_p column_name)
OBIDMS_column_p* obi_view_get_pointer_on_column_in_view(Obiview_p view, const_char_p column_name)
int obi_view_create_column_alias(Obiview_p view, const_char_p current_name, const_char_p alias)
char* obi_view_formatted_infos(Obiview_p view, bint detailed)
char* obi_view_formatted_infos_one_line(Obiview_p view)
int obi_view_write_comments(Obiview_p view, const_char_p comments)
int obi_view_add_comment(Obiview_p view, const_char_p key, const_char_p value)

View File

@ -22,6 +22,7 @@ cdef class Column(OBIWrapper) :
cdef inline OBIDMS_column_p pointer(self)
cdef read_elements_names(self)
cpdef list keys(self)
@staticmethod
cdef type get_column_class(obitype_t obitype, bint multi_elts, bint tuples)

View File

@ -14,6 +14,7 @@ from ..capi.obidms cimport obi_import_column
from ..capi.obidmscolumn cimport OBIDMS_column_header_p, \
obi_close_column, \
obi_get_elements_names, \
obi_column_formatted_infos, \
obi_column_write_comments
from ..capi.obiutils cimport obi_format_date
@ -21,7 +22,11 @@ from ..capi.obiutils cimport obi_format_date
from ..capi.obiview cimport obi_view_add_column, \
obi_view_get_pointer_on_column_in_view, \
Obiview_p, \
NUC_SEQUENCE_COLUMN
NUC_SEQUENCE_COLUMN, \
QUALITY_COLUMN, \
REVERSE_SEQUENCE_COLUMN, \
REVERSE_QUALITY_COLUMN
from ..object cimport OBIDeactivatedInstanceError
@ -34,8 +39,9 @@ from obitools3.utils cimport tobytes, \
from obitools3.dms.column import typed_column
from libc.stdlib cimport free
from libc.stdlib cimport free
from libc.string cimport strcpy
import importlib
import inspect
import pkgutil
@ -92,6 +98,7 @@ cdef class Column(OBIWrapper) :
object alias=b""):
# TODO indexer_name?
cdef Column column
cdef bytes column_name_b = tobytes(column_name)
cdef bytes alias_b = tobytes(alias)
cdef bytes comments_b = str2bytes(json.dumps(bytes2str_object(comments)))
@ -122,11 +129,19 @@ cdef class Column(OBIWrapper) :
if data_type == OBI_QUAL:
if associated_column_name_b == b"":
if NUC_SEQUENCE_COLUMN not in view:
raise RuntimeError("Cannot create column %s in view %s: trying to create quality column but no NUC_SEQ column to associate it with in the view" % (bytes2str(column_name_b),
bytes2str(view.name)))
associated_column_name_b = NUC_SEQUENCE_COLUMN
associated_column_version = view[NUC_SEQUENCE_COLUMN].version
if column_name == QUALITY_COLUMN:
if NUC_SEQUENCE_COLUMN not in view:
raise RuntimeError("Cannot create column %s in view %s: trying to create quality column but no NUC_SEQ column to associate it with in the view" % (bytes2str(column_name_b),
bytes2str(view.name)))
associated_column_name_b = NUC_SEQUENCE_COLUMN
associated_column_version = view[NUC_SEQUENCE_COLUMN].version
elif column_name == REVERSE_QUALITY_COLUMN:
if REVERSE_SEQUENCE_COLUMN not in view:
raise RuntimeError("Cannot create column %s in view %s: trying to create reverse quality column but no REVERSE_SEQUENCE column to associate it with in the view" % (bytes2str(column_name_b),
bytes2str(view.name)))
associated_column_name_b = REVERSE_SEQUENCE_COLUMN
associated_column_version = view[REVERSE_SEQUENCE_COLUMN].version
if (obi_view_add_column(view = view.pointer(),
column_name = column_name_b,
@ -146,8 +161,19 @@ cdef class Column(OBIWrapper) :
create = True)<0):
raise RuntimeError("Cannot create column %s in view %s" % (bytes2str(column_name_b),
bytes2str(view.name)))
return Column.open(view, alias_b)
column = Column.open(view, alias_b)
# Automatically associate nuc sequence column to quality column if necessary
if data_type == OBI_QUAL:
if column_name == QUALITY_COLUMN:
view[NUC_SEQUENCE_COLUMN].associated_column_name = column.name
view[NUC_SEQUENCE_COLUMN].associated_column_version = column.version
elif column_name == REVERSE_QUALITY_COLUMN:
view[REVERSE_SEQUENCE_COLUMN].associated_column_name = column.name
view[REVERSE_SEQUENCE_COLUMN].associated_column_version = column.version
return column
@staticmethod
@ -276,10 +302,25 @@ cdef class Column(OBIWrapper) :
@OBIWrapper.checkIsActive
def __repr__(self) :
cdef bytes s
s = self._alias + b", data type: " + self.data_type
return bytes2str(s)
cdef str s
cdef char* sc
cdef OBIDMS_column_p pointer = self.pointer()
sc = obi_column_formatted_infos(pointer, False)
s = bytes2str(sc)
free(sc)
return s
@OBIWrapper.checkIsActive
def repr_longformat(self) :
cdef str s
cdef char* sc
cdef OBIDMS_column_p pointer = self.pointer()
sc = obi_column_formatted_infos(pointer, True)
s = bytes2str(sc)
free(sc)
return s
def close(self): # TODO discuss, can't be called bc then bug when closing view that tries to close it in C
@ -305,7 +346,10 @@ cdef class Column(OBIWrapper) :
free(elts_names_b)
return elts_names_list
cpdef list keys(self):
return self._elements_names
# Column alias property getter and setter
@property
def name(self):
@ -322,7 +366,7 @@ cdef class Column(OBIWrapper) :
@property
def elements_names(self):
return self._elements_names
# nb_elements_per_line property getter
@property
def nb_elements_per_line(self):
@ -386,6 +430,31 @@ cdef class Column(OBIWrapper) :
raise OBIDeactivatedInstanceError()
return obi_format_date(self.pointer().header.creation_date)
# associated_column name property getter and setter
@property
def associated_column_name(self):
if not self.active() :
raise OBIDeactivatedInstanceError()
return self.pointer().header.associated_column.column_name
@associated_column_name.setter
def associated_column_name(self, object new_name):
strcpy(self.pointer().header.associated_column.column_name, tobytes(new_name))
# associated_column version property getter and setter
@property
def associated_column_version(self):
if not self.active() :
raise OBIDeactivatedInstanceError()
return self.pointer().header.associated_column.version
@associated_column_version.setter
def associated_column_version(self, int new_version):
self.pointer().header.associated_column.version = new_version
# comments property getter
@property
def comments(self):

View File

@ -10,7 +10,8 @@ from .capi.obidms cimport obi_open_dms, \
obi_dms_exists, \
obi_dms_get_full_path, \
obi_close_atexit, \
obi_dms_write_comments
obi_dms_write_comments, \
obi_dms_formatted_infos
from .capi.obitypes cimport const_char_p
@ -32,6 +33,8 @@ from .object import OBIWrapper
import json
import time
from libc.stdlib cimport free
cdef class DMS(OBIWrapper):
@ -94,16 +97,16 @@ cdef class DMS(OBIWrapper):
return dms
def close(self) :
def close(self, force=False) :
'''
Closes the DMS instance and free the associated memory
Closes the DMS instance and free the associated memory (no counter, closing is final)
The `close` method is automatically called by the object destructor.
'''
cdef OBIDMS_p pointer = self.pointer()
if self.active() :
OBIWrapper.close(self)
if (obi_close_dms(pointer, False)) < 0 :
if (obi_close_dms(pointer, force=force)) < 0 :
raise Exception("Problem closing an OBIDMS")
@ -223,11 +226,24 @@ cdef class DMS(OBIWrapper):
@OBIWrapper.checkIsActive
def __repr__(self):
cdef str s
s=""
for view_name in self.keys():
s = s + repr(self.get_view(view_name)) + "\n"
def __repr__(self) :
cdef str s
cdef char* sc
cdef OBIDMS_p pointer = self.pointer()
sc = obi_dms_formatted_infos(pointer, False)
s = bytes2str(sc)
free(sc)
return s
@OBIWrapper.checkIsActive
def repr_longformat(self) :
cdef str s
cdef char* sc
cdef OBIDMS_p pointer = self.pointer()
sc = obi_dms_formatted_infos(pointer, True)
s = bytes2str(sc)
free(sc)
return s
@ -254,12 +270,13 @@ cdef class DMS(OBIWrapper):
# bash command history property getter
@property
def bash_history(self):
s = b"#!/bin/bash\n\n"
#s = b"#!${bash}/bin/bash\n\n"
s = b""
first = True
for command in self.command_line_history:
s+=b"#"
s+=command[b"time"]
s+=b"\n"
s+=b"\nobi "
s+=command[b"command"]
s+=b"\n"
return s

View File

@ -39,4 +39,6 @@ cdef class Nuc_Seq_Stored(Seq_Stored) :
cpdef set_quality_char(self, object new_qual, int offset=*)
cpdef object build_quality_array(self, list quality)
cpdef bytes build_reverse_complement(self)
cpdef str get_str(self)
cpdef str get_str(self)
cpdef repr_bytes(self)

View File

@ -431,9 +431,12 @@ cdef class Nuc_Seq_Stored(Seq_Stored) :
return len(self._view.get_column(NUC_SEQUENCE_COLUMN).get_line(self.index))
def __repr__(self):
return bytes2str(self.repr_bytes())
cpdef repr_bytes(self):
if self.quality is None:
formatter = FastaFormat()
else:
formatter = FastqFormat()
return bytes2str(formatter(self))
return formatter(self)

View File

@ -11,11 +11,14 @@ cdef class Taxonomy(OBIWrapper) :
cdef bytes _name
cdef DMS _dms
cdef list _ranks
cdef dict _name_dict
cdef inline OBIDMS_taxonomy_p pointer(self)
cdef fill_name_dict(self)
cpdef Taxon get_taxon_by_idx(self, int idx)
cpdef Taxon get_taxon_by_taxid(self, int taxid)
cpdef Taxon get_taxon_by_name(self, object taxon_name, object restricting_taxid=*)
cpdef write(self, object prefix)
cpdef int add_taxon(self, str name, str rank_name, int parent_taxid, int min_taxid=*)
cpdef object get_species(self, int taxid)

View File

@ -15,7 +15,11 @@ from ..capi.obitaxonomy cimport obi_taxonomy_exists, \
obi_taxo_get_species, \
obi_taxo_get_genus, \
obi_taxo_get_family, \
ecotx_t
ecotx_t, \
econame_t, \
obi_taxo_get_name_from_name_idx, \
obi_taxo_get_taxon_from_name_idx
from cpython.pycapsule cimport PyCapsule_New, PyCapsule_GetPointer
import tarfile
@ -24,11 +28,29 @@ from libc.stdlib cimport free
cdef class Taxonomy(OBIWrapper) :
# TODO function to import taxonomy?
# TODO function to import taxonomy?
cdef inline OBIDMS_taxonomy_p pointer(self) :
return <OBIDMS_taxonomy_p>(self._pointer)
cdef fill_name_dict(self):
print("Indexing taxon names...")
cdef OBIDMS_taxonomy_p pointer = self.pointer()
cdef ecotx_t* taxon_p
cdef object taxon_capsule
cdef bytes name
cdef int count
cdef int n
count = (<OBIDMS_taxonomy_p>pointer).names.count
for n in range(count) :
name = obi_taxo_get_name_from_name_idx(pointer, n)
taxon_p = obi_taxo_get_taxon_from_name_idx(pointer, n)
taxon_capsule = PyCapsule_New(taxon_p, NULL, NULL)
self._name_dict[name] = Taxon(taxon_capsule, self)
@staticmethod
def exists(DMS dms, object name) :
@ -75,7 +97,8 @@ cdef class Taxonomy(OBIWrapper) :
taxo._dms = dms
taxo._name = tobytes(name)
taxo._name_dict = {}
taxo.fill_name_dict()
taxo._ranks = []
for r in range((<OBIDMS_taxonomy_p>pointer).ranks.count) :
taxo._ranks.append(obi_taxo_rank_index_to_label(r, (<OBIDMS_taxonomy_p>pointer).ranks))
@ -118,7 +141,8 @@ cdef class Taxonomy(OBIWrapper) :
taxo._dms = dms
taxo._name = folder_path
taxo._name_dict = {}
taxo.fill_name_dict()
taxo._ranks = []
for r in range((<OBIDMS_taxonomy_p>pointer).ranks.count) :
taxo._ranks.append(obi_taxo_rank_index_to_label(r, (<OBIDMS_taxonomy_p>pointer).ranks))
@ -129,8 +153,8 @@ cdef class Taxonomy(OBIWrapper) :
def __getitem__(self, object ref):
if type(ref) == int :
return self.get_taxon_by_taxid(ref)
else :
raise NotImplementedError()
elif type(ref) == str or type(ref) == bytes :
return self.get_taxon_by_name(ref)
cpdef Taxon get_taxon_by_taxid(self, int taxid):
@ -143,6 +167,19 @@ cdef class Taxonomy(OBIWrapper) :
return Taxon(taxon_capsule, self)
cpdef Taxon get_taxon_by_name(self, object taxon_name, object restricting_taxid=None):
taxon = self._name_dict.get(tobytes(taxon_name), None)
if not taxon:
return None
elif restricting_taxid:
if self.is_ancestor(restricting_taxid, taxon.taxid):
return taxon
else:
return None
else:
return taxon
cpdef Taxon get_taxon_by_idx(self, int idx):
cdef ecotx_t* taxa
cdef ecotx_t* taxon_p
@ -232,7 +269,7 @@ cdef class Taxonomy(OBIWrapper) :
taxa = self.pointer().taxa.taxon
# Yield each taxid
# Yield each taxon
for t in range(self.pointer().taxa.count):
taxon_p = <ecotx_t*> (taxa+t)
taxon_capsule = PyCapsule_New(taxon_p, NULL, NULL)

View File

@ -20,9 +20,14 @@ cdef class View(OBIWrapper):
cdef DMS _dms
cdef inline Obiview_p pointer(self)
cpdef print_to_output(self,
object output,
bint noprogressbar=*)
cpdef delete_column(self,
object column_name)
object column_name,
bint delete_file=*)
cpdef rename_column(self,
object current_name,
@ -60,6 +65,8 @@ cdef class Line :
cdef index_t _index
cdef View _view
cpdef repr_bytes(self)
cdef register_view_class(bytes view_type_name,
type view_class)

View File

@ -6,6 +6,9 @@ cdef dict __VIEW_CLASS__= {}
from libc.stdlib cimport malloc
from obitools3.apps.progress cimport ProgressBar # @UnresolvedImport
from obitools3.version import version
from ..capi.obiview cimport Alias_column_pair_p, \
obi_new_view, \
obi_open_view, \
@ -16,7 +19,9 @@ from ..capi.obiview cimport Alias_column_pair_p, \
obi_view_delete_column, \
obi_view_create_column_alias, \
obi_view_write_comments, \
obi_delete_view
obi_delete_view, \
obi_view_formatted_infos, \
obi_view_formatted_infos_one_line
from ..capi.obidmscolumn cimport OBIDMS_column_p
from ..capi.obidms cimport OBIDMS_p
@ -48,10 +53,15 @@ from ..capi.obidms cimport obi_import_view
from obitools3.format.tab import TabFormat
from cpython.exc cimport PyErr_CheckSignals
import importlib
import inspect
import pkgutil
import json
import sys
from libc.stdlib cimport free
cdef class View(OBIWrapper) :
@ -178,13 +188,50 @@ cdef class View(OBIWrapper) :
@OBIWrapper.checkIsActive
def __repr__(self) :
cdef str s = "#View name:\n{name:s}\n#Line count:\n{line_count:d}\n#Columns:\n".format(name = bytes2str(self.name),
line_count = self.line_count)
for column_name in self.keys() :
s = s + repr(self[column_name]) + '\n'
def __repr__(self) :
cdef str s
cdef char* sc
cdef Obiview_p pointer = self.pointer()
sc = obi_view_formatted_infos(pointer, False)
s = bytes2str(sc)
free(sc)
return s
@OBIWrapper.checkIsActive
def repr_longformat(self) :
cdef str s
cdef char* sc
cdef Obiview_p pointer = self.pointer()
sc = obi_view_formatted_infos(pointer, True)
s = bytes2str(sc)
free(sc)
return s
cpdef print_to_output(self, object output, bint noprogressbar=False):
cdef int i
cdef Line entry
self.checkIsActive(self)
# Initialize the progress bar
if noprogressbar == False:
pb = ProgressBar(len(self))
else:
pb = None
i=0
for entry in self:
PyErr_CheckSignals()
if pb is not None:
pb(i)
output.write(entry.repr_bytes()+b"\n")
i+=1
if pb is not None:
pb(len(self), force=True)
print("", file=sys.stderr)
def keys(self):
@ -227,7 +274,8 @@ cdef class View(OBIWrapper) :
cpdef delete_column(self,
object column_name) :
object column_name,
bint delete_file=False) :
cdef bytes column_name_b = tobytes(column_name)
@ -239,7 +287,7 @@ cdef class View(OBIWrapper) :
col.close()
# Remove the column from the view which closes the C structure
if obi_view_delete_column(self.pointer(), column_name_b) < 0 :
if obi_view_delete_column(self.pointer(), column_name_b, delete_file) < 0 :
raise RollbackException("Problem deleting column %s from a view",
bytes2str(column_name_b), self)
@ -297,11 +345,17 @@ cdef class View(OBIWrapper) :
nb_elements_per_line=new_nb_elements_per_line, elements_names=new_elements_names,
comments=old_column.comments, alias=column_name_b+tobytes('___new___'))
switch_to_dict = old_column.nb_elements_per_line == 1 and new_nb_elements_per_line > 1
ori_key = old_column._elements_names[0]
for i in range(length) :
new_column[i] = old_column[i]
if switch_to_dict :
new_column[i] = {ori_key: old_column[i]}
else:
new_column[i] = old_column[i]
# Remove old column from view
self.delete_column(column_name_b)
self.delete_column(column_name_b, delete_file=True)
# Rename new
new_column.name = column_name_b
@ -398,6 +452,7 @@ cdef class View(OBIWrapper) :
for i in range(len(input_view_name)):
input_str.append(tostr(input_dms_name[i])+"/"+tostr(input_view_name[i]))
comments["input_str"] = input_str
comments["version"] = version
return bytes2str_object(comments)
@ -519,11 +574,12 @@ cdef class View(OBIWrapper) :
# bash command history property getter
@property
def bash_history(self):
s = b"#!/bin/bash\n\n"
s = b""
first = True
for level in self.view_history:
command_list = [level[input][b"command_line"] for input in level.keys()]
for command in command_list:
s+=b"obi "
s+=command
s+=b"\n"
return s
@ -749,8 +805,12 @@ cdef class Line :
def __repr__(self):
return bytes2str(self).repr_bytes()
cpdef repr_bytes(self):
formatter = TabFormat(header=False)
return bytes2str(formatter(self))
return formatter(self)
# View property getter

View File

@ -3,8 +3,9 @@
cimport cython
from obitools3.dms.view.view cimport Line
from obitools3.utils cimport bytes2str_object, str2bytes, tobytes
from obitools3.dms.column.column cimport Column_line
from obitools3.dms.column.column cimport Column_line, Column_multi_elts
import sys
cdef class TabFormat:
@ -25,19 +26,33 @@ cdef class TabFormat:
for k in self.tags:
if self.header and self.first_line:
value = tobytes(k)
if isinstance(data.view[k], Column_multi_elts):
keys = data.view[k].keys()
keys.sort()
for k2 in keys:
line.append(tobytes(k)+b':'+tobytes(k2))
else:
line.append(tobytes(k))
else:
value = data[k]
if value is not None:
if type(value) == Column_line:
value = value.bytes()
if isinstance(data.view[k], Column_multi_elts):
keys = data.view[k].keys()
keys.sort()
if value is None: # all keys at None
for k2 in keys: # TODO could be much more efficient
line.append(self.NAString)
else:
value = str2bytes(str(bytes2str_object(value))) # genius programming
if value is None:
value = self.NAString
line.append(value)
for k2 in keys: # TODO could be much more efficient
if value[k2] is not None:
line.append(str2bytes(str(bytes2str_object(value[k2])))) # genius programming
else:
line.append(self.NAString)
else:
if value is not None:
line.append(str2bytes(str(bytes2str_object(value))))
else:
line.append(self.NAString)
if self.first_line:
self.first_line = False

View File

@ -6,6 +6,7 @@ from .solexapairend import iterOnAligment
from .shifted_ali cimport Ali_shifted
from obitools3.dms.capi.obiview cimport Obiview_p, QUALITY_COLUMN, NUC_SEQUENCE_COLUMN, \
REVERSE_SEQUENCE_COLUMN, REVERSE_QUALITY_COLUMN, \
obi_set_qual_int_with_elt_idx_and_col_p_in_view, \
obi_set_str_with_elt_idx_and_col_p_in_view
@ -13,7 +14,6 @@ from obitools3.dms.capi.obidmscolumn cimport OBIDMS_column_p
from obitools3.dms.view.view cimport View
from obitools3.dms.column.column cimport Column
from obitools3.commands.ngsfilter import REVERSE_SEQ_COLUMN_NAME, REVERSE_QUALITY_COLUMN_NAME
from math import log10
@ -183,12 +183,13 @@ def buildConsensus(ali, seq, ref_tags=None):
# doesn't work because uint8_t* are forced into bytes by cython (nothing read/kept beyond 0 values)
#obi_set_qual_int_with_elt_idx_and_col_p_in_view(view_p, col_p, seq.index, 0, ali.consensus_qual, ali.consensus_len)
seq.set(ref_tags.id+b"_CONS", ali.consensus_seq, quality=ali.consensus_qual)
seq[b'ali_length'] = ali.consensus_len
seq[b'score_norm']=float(ali.score)/ali.consensus_len
seq[b"seq_length"] = ali.consensus_len
seq[b"overlap_length"] = ali.overlap_len
seq[b'score_norm']=round(float(ali.score)/ali.overlap_len, 3)
seq[b'shift']=ali.shift
else:
if len(ali[0])>999: # TODO why?
raise AssertionError,"Too long alignemnt"
raise AssertionError,"Too long alignment"
ic=IterOnConsensus(ali)
@ -233,7 +234,7 @@ def buildConsensus(ali, seq, ref_tags=None):
seq[b'mode']=b'alignment'
for tag in ref_tags:
if tag != REVERSE_SEQ_COLUMN_NAME and tag != REVERSE_QUALITY_COLUMN_NAME and \
if tag != REVERSE_SEQUENCE_COLUMN and tag != REVERSE_QUALITY_COLUMN and \
tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN:
seq[tag] = ref_tags[tag]
@ -250,11 +251,22 @@ def buildJoinedSequence(ali, reverse, seq, forward=None):
quality.extend(reverse.quality)
seq.set(forward.id +b"_PairedEnd", s, definition=forward.definition, quality=quality)
seq[b"score"]=ali.score
seq[b"ali_direction"]=ali.direction
if len(ali.direction) > 0:
seq[b"ali_direction"]=ali.direction
else:
seq[b"ali_direction"]=None
seq[b"mode"]=b"joined"
seq[b"pairedend_limit"]=len(forward)
seq[b"pairedend_limit"]=len(forward)
seq[b"seq_length"] = ali.consensus_len
seq[b"overlap_length"] = ali.overlap_len
if ali.overlap_len > 0:
seq[b'score_norm']=round(float(ali.score)/ali.overlap_len, 3)
else:
seq[b"score_norm"]=0.0
for tag in forward:
if tag != REVERSE_SEQ_COLUMN_NAME and tag != REVERSE_QUALITY_COLUMN_NAME:
if tag != REVERSE_SEQUENCE_COLUMN and tag != REVERSE_QUALITY_COLUMN and \
tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN:
seq[tag] = forward[tag]
return seq

View File

@ -156,6 +156,9 @@ def emblIterator_file(lineiterator,
yield seq
read+=1
# Last sequence
seq = emblParser(entry)
yield seq
free(entry)
@ -174,7 +177,7 @@ def emblIterator_dir(dir_path,
for filename in files:
if read==only:
return
print("Parsing file %s (%d/%d)" % (tostr(filename), read_files, len(files)))
print("Parsing file %s (%d/%d)" % (tostr(filename), read_files+1, len(files)))
f = uopen(filename)
if only is not None:
only_f = only-read

View File

@ -104,6 +104,7 @@ def fastaNucIterator(lineiterator,
cdef bytes sequence
cdef int skipped, ionly, read
cdef Nuc_Seq seq
cdef bint stop
if only is None:
ionly = -1
@ -130,7 +131,8 @@ def fastaNucIterator(lineiterator,
else:
line = firstline
while True:
stop=False
while not stop:
if ionly >= 0 and read >= ionly:
break
@ -153,7 +155,7 @@ def fastaNucIterator(lineiterator,
s.append(line[0:-1])
line = next(iterator)
except StopIteration:
pass
stop=True
sequence = b"".join(s)

View File

@ -25,8 +25,9 @@ from libc.string cimport strcpy, strlen
_featureMatcher = re.compile(b'^FEATURES.+\n(?=ORIGIN)',re.DOTALL + re.M)
_headerMatcher = re.compile(b'^LOCUS.+(?=\nFEATURES)', re.DOTALL + re.M)
_seqMatcher = re.compile(b'(?<=ORIGIN).+(?=//\n)', re.DOTALL + re.M)
_cleanSeq = re.compile(b'[ \n0-9]+')
_seqMatcher = re.compile(b'ORIGIN.+(?=//\n)', re.DOTALL + re.M)
_cleanSeq1 = re.compile(b'ORIGIN.+\n')
_cleanSeq2 = re.compile(b'[ \n0-9]+')
_acMatcher = re.compile(b'(?<=^ACCESSION ).+',re.M)
_deMatcher = re.compile(b'(?<=^DEFINITION ).+\n( .+\n)*',re.M)
_cleanDe = re.compile(b'\n *')
@ -42,7 +43,8 @@ def genbankParser(bytes text):
ft = _featureMatcher.search(text).group()
s = _seqMatcher.search(text).group()
s = _cleanSeq.sub(b'', s).upper()
s = _cleanSeq1.sub(b'', s)
s = _cleanSeq2.sub(b'', s)
acs = _acMatcher.search(text).group()
acs = acs.split()
@ -51,23 +53,23 @@ def genbankParser(bytes text):
de = _deMatcher.search(header).group()
de = _cleanDe.sub(b' ',de).strip().strip(b'.')
tags = {}
extractTaxon(ft, tags)
seq = Nuc_Seq(ac,
s,
definition=de,
quality=None,
offset=-1,
tags=tags)
except Exception as e:
print("\nCould not import sequence id:", text.split()[1], "(error raised:", e, ")")
# Do not raise any Exception if you need the possibility to resume the generator
# (Python generators can't resume after any exception is raised)
return None
tags = {}
extractTaxon(ft, tags)
seq = Nuc_Seq(ac,
s,
definition=de,
quality=None,
offset=-1,
tags=tags)
return seq
@ -153,6 +155,9 @@ def genbankIterator_file(lineiterator,
yield seq
read+=1
# Last sequence
seq = genbankParser(entry)
yield seq
free(entry)
@ -168,10 +173,12 @@ def genbankIterator_dir(dir_path,
read = 0
read_files = 0
files = [filename for filename in glob.glob(os.path.join(path, b'*.gbff*'))]
files.extend([filename for filename in glob.glob(os.path.join(path, b'*.seq*'))]) # new genbank extension
files = list(set(files))
for filename in files:
if read==only:
return
print("Parsing file %s (%d/%d)" % (tostr(filename), read_files, len(files)))
print("Parsing file %s (%d/%d)" % (tostr(filename), read_files+1, len(files)))
f = uopen(filename)
if only is not None:
only_f = only-read

3
python/obitools3/parsers/ngsfilter.pyx Executable file → Normal file
View File

@ -57,6 +57,9 @@ def ngsfilterIterator(lineiterator,
split_line = line.split()
tags = split_line.pop(2)
tags = tags.split(b":")
for t_idx in range(len(tags)):
if tags[t_idx]==b"-" or tags[t_idx]==b"None" or tags[t_idx]==b"":
tags[t_idx] = nastring
if len(tags) == 1: # Forward and reverse tags are the same
tags.append(tags[0])
split_line.insert(2, tags[0])

26
python/obitools3/uri/decode.pyx Executable file → Normal file
View File

@ -171,7 +171,8 @@ Reads an URI and returns a tuple containing:
def open_uri(uri,
bint input=True,
type newviewtype=View,
dms_only=False):
dms_only=False,
force_file=False):
cdef bytes urib = tobytes(uri)
cdef bytes scheme
@ -195,9 +196,9 @@ def open_uri(uri,
if 'obi' not in config:
config['obi']={}
try:
if not force_file and "defaultdms" in config["obi"]:
default_dms=config["obi"]["defaultdms"]
except KeyError:
else:
default_dms=None
try:
@ -209,10 +210,11 @@ def open_uri(uri,
error = None
if 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
if 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
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)
@ -274,11 +276,11 @@ def open_uri(uri,
iseq = urib
objclass = bytes
else: # TODO update uopen to be able to write?
if urip.path:
file = open(urip.path, 'wb')
else:
if not urip.path or urip.path == b'-':
file = sys.stdout.buffer
else:
file = open(urip.path, 'wb')
if file is not None:
qualifiers=parse_qs(urip.query)
@ -462,7 +464,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":
if format==b"fasta" or format==b"silva":
if input:
iseq = fastaNucIterator(file,
skip=skip,

View File

@ -2,7 +2,7 @@
from obitools3.dms.capi.obitypes cimport obitype_t, index_t
cpdef bytes format_separator(bytes format)
cpdef bytes format_uniq_pattern(bytes format)
cpdef int count_entries(file, bytes format)
cdef obi_errno_to_exception(index_t line_nb=*, object elt_id=*, str error_message=*)

View File

@ -24,11 +24,11 @@ import glob
import gzip
cpdef bytes format_separator(bytes format):
cpdef bytes format_uniq_pattern(bytes format):
if format == b"fasta":
return b"\n>"
elif format == b"fastq":
return b"\n@"
return b"\n\+\n"
elif format == b"ngsfilter" or format == b"tabular":
return b"\n"
elif format == b"genbank" or format == b"embl":
@ -42,7 +42,7 @@ cpdef bytes format_separator(bytes format):
cpdef int count_entries(file, bytes format):
try:
sep = format_separator(format)
sep = format_uniq_pattern(format)
if sep is None:
return -1
sep = re.compile(sep)
@ -72,7 +72,7 @@ cpdef int count_entries(file, bytes format):
return -1
mmapped_file = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
total_count += len(re.findall(sep, mmapped_file))
if format != b"ngsfilter" and format != b"tabular":
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)
except:
@ -166,7 +166,9 @@ cdef object bytes2str_object(object value): # Only works if complex types are d
value[k] = bytes2str(v)
if type(k) == bytes:
value[bytes2str(k)] = value.pop(k)
elif isinstance(value, list):
elif isinstance(value, list) or isinstance(value, tuple):
if isinstance(value, tuple):
value = list(value)
for i in range(len(value)):
if isinstance(value[i], list) or isinstance(value[i], dict):
value[i] = bytes2str_object(value[i])

View File

@ -1,5 +1,5 @@
major = 3
minor = 0
serial= '0-beta1'
serial= '0b37'
version ="%d.%02d.%s" % (major,minor,serial)
version ="%d.%d.%s" % (major,minor,serial)

5
requirements.txt Executable file
View File

@ -0,0 +1,5 @@
--extra-index-url https://pypi.python.org/simple/
Cython>=0.24
Sphinx>=1.2.0
ipython>=3.0.0
breathe>=4.0.0

View File

@ -5,8 +5,9 @@ import re
import subprocess
from distutils import log
from distutils.core import setup
#from distutils.core import setup
from setuptools import setup # to work with pip
from distutils.core import Extension
from distutils.sysconfig import get_python_lib
@ -16,6 +17,8 @@ from distutils.extension import Extension
from distutils.dist import Distribution as ori_Distribution
from python.obitools3.version import version
class Distribution(ori_Distribution):
@ -24,10 +27,11 @@ class Distribution(ori_Distribution):
ori_Distribution.__init__(self, attrs)
self.global_options.insert(0,('cobitools3', None, "intall location of the C library"
self.global_options.insert(0,('cobitools3', None, "install location of the C library"
))
from distutils.command.build import build as build_ori
from setuptools.command.bdist_egg import bdist_egg as bdist_egg_ori
from distutils.core import Command
@ -68,6 +72,12 @@ class build(build_ori):
build_ori.run(self)
class bdist_egg(bdist_egg_ori):
def run(self):
self.run_command('build_clib')
bdist_egg_ori.run(self)
sys.path.append(os.path.abspath("python"))
@ -83,12 +93,13 @@ def findPackage(root,base=None):
PACKAGE = "OBITools3"
VERSION = "3.0.0-beta1"
VERSION = version
AUTHOR = 'Celine Mercier'
EMAIL = 'celine.mercier@metabarcoding.org'
URL = "http://metabarcoding.org/obitools3"
URL = "https://metabarcoding.org/obitools3"
PLATFORMS = "posix"
LICENSE = "CeCILL-V2"
DESCRIPTION = "Tools and library for DNA metabarcoding",
DESCRIPTION = "A package for the management of analyses and data in DNA metabarcoding."
PYTHONMIN = '3.5'
SRC = 'python'
@ -145,17 +156,24 @@ classifiers=['Development Status :: 4 - Beta',
'Topic :: Utilities',
]
with open("README.md", "r") as fh:
long_description = fh.read()
setup(name=PACKAGE,
description=DESCRIPTION,
long_description=long_description,
long_description_content_type="text/markdown",
classifiers=classifiers,
version=VERSION,
author=AUTHOR,
author_email=EMAIL,
platforms=PLATFORMS,
license=LICENSE,
url=URL,
ext_modules=xx,
distclass=Distribution,
cmdclass={'build': build,
'bdist_egg': bdist_egg,
'build_clib': build_clib},
cobitools3=get_python_lib(),
packages = findPackage('python'),

View File

@ -157,7 +157,7 @@ int build_reference_db(const char* dms_name,
ecotx_t* lca_2 = NULL;
ecotx_t* lca = NULL;
index_t idx1, idx2;
index_t i, j, k;
index_t i, j, k, count;
int32_t taxid_array_length;
int32_t score_array_length;
int32_t taxid_array_writable_length;
@ -185,6 +185,7 @@ int build_reference_db(const char* dms_name,
matrix_view_name = strcpy(matrix_view_name, o_view_name);
strcat(matrix_view_name, "_matrix");
fprintf(stderr, "Aligning queries with reference database...\n");
if (obi_lcs_align_one_column(dms_name,
refs_view_name,
"",
@ -320,13 +321,19 @@ int build_reference_db(const char* dms_name,
return -1;
}
count = (matrix_with_lca_view->infos)->line_count;
fprintf(stderr, "Computing LCAs...\n");
// Compute all the LCAs
// For each pair
for (i=0; i<(matrix_with_lca_view->infos)->line_count; i++)
for (i=0; i<count; i++)
{
if (! keep_running)
return -1;
if (i%1000 == 0)
fprintf(stderr,"\rDone : %f %% ", (i / (float) count)*100);
// Read all taxids associated with the first sequence and compute their LCA
// Read line index
idx1 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx1_column, i, 0);
@ -363,6 +370,7 @@ int build_reference_db(const char* dms_name,
return -1;
}
}
fprintf(stderr,"\rDone : 100 %% \n");
// Clone refs view, add 2 arrays columns for lca and score, compute and write them
@ -442,13 +450,18 @@ int build_reference_db(const char* dms_name,
return -1;
}
fprintf(stderr, "Building LCA arrays...\n");
// For each sequence, look for all its alignments in the matrix, and for each different LCA taxid/score, order them and write them
// Going through matrix once, filling refs arrays on the go for efficiency
for (i=0; i<(matrix_with_lca_view->infos)->line_count; i++)
for (i=0; i<count; i++)
{
if (! keep_running)
return -1;
if (i%1000 == 0)
fprintf(stderr,"\rDone : %f %% ", (i / (float) count)*100);
// Read ref line indexes
idx1 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx1_column, i, 0);
idx2 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx2_column, i, 0);
@ -464,6 +477,8 @@ int build_reference_db(const char* dms_name,
// Read alignment score
score = obi_get_float_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_score_column, i, 0);
//fprintf(stderr, "\n\ntaxid_lca=%d, score=%f, idx1=%d, idx2=%d", taxid_lca, score, idx1, idx2);
///////////////// Compute for first sequence \\\\\\\\\\\\\\\\\\\\\\\ (TODO function)
// Read arrays
@ -480,9 +495,11 @@ int build_reference_db(const char* dms_name,
// return -1;
// }
//fprintf(stderr, "\n1st sequence");
// If empty, add values
if (taxid_array_length == 0)
{
//fprintf(stderr, "\nEmpty, add value");
if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx1, &taxid_lca, (uint8_t) (obi_sizeof(OBI_INT) * 8), 1) < 0)
{
obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@ -496,6 +513,8 @@ int build_reference_db(const char* dms_name,
}
else
{
//fprintf(stderr, "\nNot empty");
j = 0;
modified = false;
while (j < taxid_array_length)
@ -509,6 +528,9 @@ int build_reference_db(const char* dms_name,
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
//fprintf(stderr, "\nSame LCA, replace %d and %f with %d and %f", lca_taxid_array_writable[j],
// score_array_writable[j], taxid_lca, score);
// Better score for the same LCA, replace this LCA/score pair
lca_taxid_array_writable[j] = taxid_lca;
score_array_writable[j] = score;
@ -535,6 +557,8 @@ int build_reference_db(const char* dms_name,
{
if (score > score_array[j])
{
//fprintf(stderr, "\nInsert new");
memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
@ -579,10 +603,15 @@ int build_reference_db(const char* dms_name,
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
//fprintf(stderr, "\nAppend at the end");
// Append LCA
lca_taxid_array_writable[taxid_array_writable_length] = taxid_lca;
score_array_writable[score_array_writable_length] = score;
taxid_array_writable_length++;
score_array_writable_length++;
// Remove the previous (children) LCAs from the array if their score is equal or lower
while ((j>0) && (score_array_writable[j-1] <= score))
{
@ -603,6 +632,13 @@ int build_reference_db(const char* dms_name,
// Write new arrays
if (modified)
{
// fprintf(stderr, "\n\nnew array:");
// for (k=0;k<taxid_array_writable_length;k++)
// {
// lca = obi_taxo_get_taxon_with_taxid(tax, lca_taxid_array_writable[k]);
// fprintf(stderr, "\nLCA=%d, %s, score=%f", lca_taxid_array_writable[k], lca->name, score_array_writable[k]);
// }
if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx1, lca_taxid_array_writable, (uint8_t) (obi_sizeof(OBI_INT) * 8), taxid_array_writable_length) < 0)
{
obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@ -632,9 +668,13 @@ int build_reference_db(const char* dms_name,
// return -1;
// }
//fprintf(stderr, "\n2nd sequence");
// If empty, add values
if (taxid_array_length == 0)
{
//fprintf(stderr, "\nEmpty, add value");
if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx2, &taxid_lca, (uint8_t) (obi_sizeof(OBI_INT) * 8), 1) < 0)
{
obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@ -648,6 +688,8 @@ int build_reference_db(const char* dms_name,
}
else
{
//fprintf(stderr, "\nNot empty");
j = 0;
modified = false;
while (j < taxid_array_length)
@ -661,6 +703,9 @@ int build_reference_db(const char* dms_name,
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
//fprintf(stderr, "\nSame LCA, replace %d and %f with %d and %f", lca_taxid_array_writable[j],
// score_array_writable[j], taxid_lca, score);
// Better score for the same LCA, replace this LCA/score pair
lca_taxid_array_writable[j] = taxid_lca;
score_array_writable[j] = score;
@ -687,6 +732,8 @@ int build_reference_db(const char* dms_name,
{
if (score > score_array[j])
{
//fprintf(stderr, "\nInsert new");
memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
@ -727,6 +774,8 @@ int build_reference_db(const char* dms_name,
if (j == taxid_array_length) // same or parent LCA not found, need to be appended at the end
{
//fprintf(stderr, "\nAppend at the end");
memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
modified = true;
@ -735,6 +784,9 @@ int build_reference_db(const char* dms_name,
lca_taxid_array_writable[taxid_array_writable_length] = taxid_lca;
score_array_writable[score_array_writable_length] = score;
taxid_array_writable_length++;
score_array_writable_length++;
// Remove the previous (children) LCAs from the array if their score is equal or lower
while ((j>0) && (score_array_writable[j-1] <= score))
{
@ -769,11 +821,17 @@ int build_reference_db(const char* dms_name,
}
}
}
fprintf(stderr,"\rDone : 100 %% \n");
fprintf(stderr, "Writing results...\n");
count = (o_view->infos)->line_count;
// Fill empty LCA informations (because filling from potentially sparse alignment matrix) with the sequence taxid
score=1.0; // technically getting LCA of identical sequences
for (i=0; i<(o_view->infos)->line_count; i++)
for (i=0; i<count; i++)
{
if (i%1000 == 0)
fprintf(stderr,"\rDone : %f %% ", (i / (float) count)*100);
obi_get_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, i, &taxid_array_length);
if (taxid_array_length == 0) // no LCA set
{
@ -799,6 +857,7 @@ 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);
@ -858,7 +917,6 @@ int build_reference_db(const char* dms_name,
free(matrix_view_name);
free(matrix_with_lca_view_name);
fprintf(stderr,"\rDone : 100 %% \n");
return 0;
}

View File

@ -36,10 +36,12 @@ bool only_ATGC(const char* seq)
{
if (!((*c == 'A') || \
(*c == 'T') || \
(*c == 'U') || \
(*c == 'G') || \
(*c == 'C') || \
(*c == 'a') || \
(*c == 't') || \
(*c == 'u') || \
(*c == 'g') || \
(*c == 'c')))
{
@ -182,6 +184,8 @@ byte_t* encode_seq_on_2_bits(const char* seq, int32_t length)
break;
case 't':
case 'T':
case 'u':
case 'U':
seq_b[i/4] |= NUC_T_2b;
break;
default:
@ -288,6 +292,8 @@ byte_t* encode_seq_on_4_bits(const char* seq, int32_t length)
break;
case 't':
case 'T':
case 'u': // discussable
case 'U':
seq_b[i/2] |= NUC_T_4b;
break;
case 'r':

View File

@ -64,7 +64,7 @@ enum
/**
* @brief Checks if there are only 'atgcATGC' characters in a
* @brief Checks if there are only 'atgcuATGCU' characters in a
* character string.
*
* @param seq The sequence to check.
@ -129,12 +129,13 @@ byte_t get_nucleotide_from_encoded_seq(byte_t* seq, int32_t idx, uint8_t encodin
/**
* @brief Encodes a DNA sequence with each nucleotide coded on 2 bits.
*
* A or a : 00
* C or c : 01
* T or t : 10
* G or g : 11
* A or a : 00
* C or c : 01
* T or t or U or u : 10
* G or g : 11
*
* @warning The DNA sequence must contain only 'atgcATGC' characters.
* @warning The DNA sequence must contain only 'atgcuATGCU' characters.
* @warning Uracil ('U') bases are encoded as Thymine ('T') bases.
*
* @param seq The sequence to encode.
* @param length The length of the sequence to encode.
@ -169,23 +170,24 @@ char* decode_seq_on_2_bits(byte_t* seq_b, int32_t length_seq);
/**
* @brief Encodes a DNA sequence with each nucleotide coded on 4 bits.
*
* A or a : 0001
* C or c : 0010
* G or g : 0011
* T or t : 0100
* R or r : 0101
* Y or y : 0110
* S or s : 0111
* W or w : 1000
* K or k : 1001
* M or m : 1010
* B or b : 1011
* D or d : 1100
* H or h : 1101
* V or v : 1110
* N or n : 1111
* A or a : 0001
* C or c : 0010
* G or g : 0011
* T or t or U or u : 0100
* R or r : 0101
* Y or y : 0110
* S or s : 0111
* W or w : 1000
* K or k : 1001
* M or m : 1010
* B or b : 1011
* D or d : 1100
* H or h : 1101
* V or v : 1110
* N or n : 1111
*
* @warning The DNA sequence must contain only IUPAC characters.
* @warning Uracil ('U') bases are encoded as Thymine ('T') bases.
*
* @param seq The sequence to encode.
* @param length The length of the sequence to encode.

View File

@ -413,7 +413,13 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
return NULL;
}
score = max_common_kmers + kmer_size - 1; // aka the number of nucleotides in the longest stretch of kmers perfectly matching
if (max_common_kmers > 0)
score = max_common_kmers + kmer_size - 1; // aka an approximation of the number of nucleotides matching in the overlap of the alignment.
// It's an approximation because one mismatch produces kmer_size kmer mismatches if in the middle of the overlap,
// and less for mismatches located towards the ends of the overlap. The case where there are the most mismatches is assumed,
// meaning that the score will be often underestimated and never overestimated.
else
score = 0;
abs_shift = abs(best_shift);
// Save result in Obi_ali structure
@ -423,10 +429,15 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
ali->shift = abs_shift;
ali->consensus_seq = NULL;
ali->consensus_qual = NULL;
if (((best_shift <= 0) && (!switched_seqs)) || ((best_shift > 0) && switched_seqs))
strcpy(ali->direction, "left");
if (score == 0)
ali->direction[0] = '\0';
else
strcpy(ali->direction, "right");
{
if (((best_shift <= 0) && (!switched_seqs)) || ((best_shift > 0) && switched_seqs))
strcpy(ali->direction, "left");
else
strcpy(ali->direction, "right");
}
// Build the consensus sequence if asked
if (build_consensus)

View File

@ -27,7 +27,11 @@
* @brief Alignment structure, with informations about the similarity and to rebuild the alignment.
*/
typedef struct Obi_ali {
int score; /**< Alignment score, corresponding to the number of nucleotides in the longest stretch of kmers perfectly matching.
int score; /**< Alignment score, corresponding to an approximation of the number of
* nucleotides matching in the overlap of the alignment.
* It's an approximation because one mismatch produces kmer_size kmer mismatches if in the middle of the overlap,
* and less for mismatches located towards the ends of the overlap. The case where there are the most mismatches is assumed,
* meaning that the score will be often underestimated and never overestimated.
*/
int consensus_length; /**< Length of the final consensus sequence.
*/

View File

@ -246,7 +246,16 @@ int obi_clean(const char* dms_name,
// Open the sample column if there is one
if ((strcmp(sample_column_name, "") == 0) || (sample_column_name == NULL))
sample_column = NULL;
{
fprintf(stderr, "Info: No sample information provided, assuming one sample.\n");
sample_column = obi_view_get_column(i_view, COUNT_COLUMN);
if (sample_column == NULL)
{
obidebug(1, "\nError getting the COUNT column");
return -1;
}
sample_count = 1;
}
else
{
sample_column = obi_view_get_column(i_view, sample_column_name);
@ -255,6 +264,13 @@ int obi_clean(const char* dms_name,
obidebug(1, "\nError getting the sample column");
return -1;
}
sample_count = (sample_column->header)->nb_elements_per_line;
// Check that the sample column is a merged column with all sample informations
if (sample_count == 1)
{
obidebug(1, "\n\nError: If a sample column is provided, it must contain 'merged' sample counts as built by obi uniq with the -m option\n");
return -1;
}
}
// Create the output view, or a temporary one if heads only
@ -279,8 +295,6 @@ int obi_clean(const char* dms_name,
return -1;
}
sample_count = (sample_column->header)->nb_elements_per_line;
// Create the output columns
if (create_output_columns(o_view, sample_column, sample_count) < 0)
{
@ -409,8 +423,7 @@ int obi_clean(const char* dms_name,
stop = true;
}
#pragma omp parallel default(none) \
shared(thread_count, seq_count, blob_array, complete_sample_count_array, alignment_result_array, \
#pragma omp parallel shared(thread_count, seq_count, blob_array, complete_sample_count_array, alignment_result_array, \
stop, blob1, i, obi_errno, keep_running, stderr, max_ratio, iseq_column, i_view, \
similarity_mode, reference, normalize, threshold, ktable, status_column, o_view, sample_count)
{
@ -550,7 +563,7 @@ int obi_clean(const char* dms_name,
if (heads_only)
{
line_selection = malloc((o_view->infos)->line_count * sizeof(index_t));
line_selection = malloc((((o_view->infos)->line_count) + 1) * sizeof(index_t));
if (line_selection == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);

View File

@ -52,7 +52,8 @@
*
* @param dms A pointer on an OBIDMS.
* @param i_view_name The name of the input view.
* @param sample_column_name The name of the OBI_STR column in the input view where the sample information is kept.
* @param sample_column_name The name of the column in the input view where the sample information is kept.
* Must be merged informations as built by the obi uniq tool (checked by the function).
* NULL or "" (empty string) if there is no sample information.
* @param o_view_name The name of the output view where the results should be written (should not already exist).
* @param o_view_comments The comments that should be associated with the output view.

70
src/obi_ecopcr.c Executable file → Normal file
View File

@ -77,7 +77,8 @@ static int create_output_columns(Obiview_p o_view, bool kingdom_mode);
* @param err2 The number of errors in the second primer.
* @param strand The DNA strand direction of the amplicon (R(everse) or D(irect)).
* @param kingdom_mode Whether the kingdom or the superkingdom informations should be printed to the output.
* @param keep_nucleotides Number of nucleotides kept on each side of the amplicon.
* @param keep_nucleotides Number of nucleotides kept on each side of the amplicon (not including the primers if they are kept).
* @param keep_primers Whether to keep the primers.
* @param i_id_column A pointer on the input sequence identifier column.
* @param o_id_column A pointer on the output sequence identifier column.
* @param o_ori_seq_len_column A pointer on the original sequence length column.
@ -104,7 +105,8 @@ static int create_output_columns(Obiview_p o_view, bool kingdom_mode);
* @param o_temp1_column A pointer on the output column for the temperature for the first primer.
* @param o_temp2_column A pointer on the output column for the temperature for the second primer.
*
* @retval 0 if the operation was successfully completed.
* @retval 0 if the sequence was skipped (taxid not found, warning printed).
* @retval 1 if the sequence was successfully printed to the output.
* @retval -1 if an error occurred.
*
* @since July 2018
@ -124,6 +126,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
int32_t err1, int32_t err2,
char strand, bool kingdom_mode,
int keep_nucleotides,
bool keep_primers,
OBIDMS_column_p i_id_column, OBIDMS_column_p o_id_column, OBIDMS_column_p o_ori_seq_len_column,
OBIDMS_column_p o_amplicon_column, OBIDMS_column_p o_amplicon_length_column,
OBIDMS_column_p o_taxid_column, OBIDMS_column_p o_rank_column, OBIDMS_column_p o_name_column,
@ -328,6 +331,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
int32_t err1, int32_t err2,
char strand, bool kingdom_mode,
int keep_nucleotides,
bool keep_primers,
OBIDMS_column_p i_id_column, OBIDMS_column_p o_id_column, OBIDMS_column_p o_ori_seq_len_column,
OBIDMS_column_p o_amplicon_column, OBIDMS_column_p o_amplicon_length_column,
OBIDMS_column_p o_taxid_column, OBIDMS_column_p o_rank_column, OBIDMS_column_p o_name_column,
@ -363,6 +367,17 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
// TODO add check for primer longer than MAX_PAT_LEN (32)
// Get sequence id
seq_id = obi_get_str_with_elt_idx_and_col_p_in_view(i_view, i_id_column, i_idx, 0);
// Get the taxon structure
main_taxon = obi_taxo_get_taxon_with_taxid(taxonomy, taxid);
if (main_taxon == NULL)
{
obidebug(1, "\nWarning: error reading the taxonomic information of a sequence. Seq id: %s, taxid: %d. Probably deprecated taxid. Skipping this sequence.", seq_id, taxid);
return 0;
}
ldelta = (pos1 <= keep_nucleotides)?pos1:keep_nucleotides;
rdelta = ((pos2+keep_nucleotides)>=seq_len)?seq_len-pos2:keep_nucleotides;
@ -382,7 +397,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
oligo2[o1->patlen] = 0;
error2 = err1;
if (keep_nucleotides == 0)
if (!keep_primers)
amplicon+=o2->patlen;
else
{
@ -401,7 +416,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
oligo2[o2->patlen] = 0;
error2 = err2;
if (keep_nucleotides==0)
if (!keep_primers)
amplicon+=o1->patlen;
else
{
@ -411,16 +426,11 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
}
ecoComplementSequence(oligo2);
if (keep_nucleotides == 0)
if (!keep_primers)
amplicon[amplicon_len]=0;
else
{
amplicon_len = ldelta+rdelta+amplicon_len;
for (i=0; i<ldelta; i++)
amplicon[i]|=32;
for (i=1; i<=rdelta; i++)
amplicon[amplicon_len-i]|=32;
amplicon[amplicon_len] = 0;
}
@ -433,16 +443,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
if (isnan(tm2))
tm2 = OBIFloat_NA;
// Get the taxon structure
main_taxon = obi_taxo_get_taxon_with_taxid(taxonomy, taxid);
if (main_taxon == NULL)
{
obidebug(1, "\nError reading the taxonomic information of a sequence");
return -1;
}
// Write sequence id
seq_id = obi_get_str_with_elt_idx_and_col_p_in_view(i_view, i_id_column, i_idx, 0);
if (obi_set_str_with_elt_idx_and_col_p_in_view(o_view, o_id_column, o_idx, 0, seq_id) < 0)
{
obidebug(1, "\nError writing the sequence id");
@ -631,7 +632,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
return -1;
}
return 0;
return 1;
}
@ -659,6 +660,7 @@ int obi_ecopcr(const char* i_dms_name,
double salt,
int saltmethod,
int keep_nucleotides,
bool keep_primers,
bool kingdom_mode)
{
@ -699,6 +701,7 @@ int obi_ecopcr(const char* i_dms_name,
obiint_t taxid;
char* sequence;
int printed;
SeqPtr apatseq=NULL;
int32_t o1Hits;
@ -717,6 +720,9 @@ int obi_ecopcr(const char* i_dms_name,
signal(SIGINT, sig_handler);
if (keep_nucleotides > 0)
keep_primers = true;
if (circular)
{
circular = strlen(primer1);
@ -1055,14 +1061,14 @@ int obi_ecopcr(const char* i_dms_name,
length = 0;
if (posj > posi)
length = posj - posi - o1->patlen - o2->patlen;
if (posj < posi)
else if (circular > 0)
length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
if ((length>0) && // For when primers touch or overlap
(!min_len || (length >= min_len)) &&
(!max_len || (length <= max_len)))
{
// Print the found amplicon
if (print_seq(i_view, o_view,
printed = print_seq(i_view, o_view,
i_idx, o_idx,
taxonomy,
sequence,
@ -1076,6 +1082,7 @@ int obi_ecopcr(const char* i_dms_name,
erri, errj,
'D', kingdom_mode,
keep_nucleotides,
keep_primers,
i_id_column, o_id_column, o_ori_seq_len_column,
o_amplicon_column, o_amplicon_length_column,
o_taxid_column, o_rank_column, o_name_column,
@ -1087,12 +1094,14 @@ int obi_ecopcr(const char* i_dms_name,
o_strand_column,
o_primer1_column, o_primer2_column,
o_error1_column, o_error2_column,
o_temp1_column, o_temp2_column) < 0)
o_temp1_column, o_temp2_column);
if (printed < 0)
{
obidebug(1, "\nError writing the ecopcr result");
return -1;
}
o_idx++;
else if (printed > 0)
o_idx++;
}
}
}
@ -1142,14 +1151,14 @@ int obi_ecopcr(const char* i_dms_name,
length = 0;
if (posj > posi)
length = posj - posi + 1 - o2->patlen - o1->patlen; /* - o1->patlen : deleted by <EC> (prior to the OBITools3) */
if (posj < posi)
else if (circular > 0)
length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
if ((length>0) && // For when primers touch or overlap
(!min_len || (length >= min_len)) &&
(!max_len || (length <= max_len)))
{
// Print the found amplicon
if (print_seq(i_view, o_view,
printed = print_seq(i_view, o_view,
i_idx, o_idx,
taxonomy,
sequence,
@ -1163,6 +1172,7 @@ int obi_ecopcr(const char* i_dms_name,
erri, errj,
'R', kingdom_mode,
keep_nucleotides,
keep_primers,
i_id_column, o_id_column, o_ori_seq_len_column,
o_amplicon_column, o_amplicon_length_column,
o_taxid_column, o_rank_column, o_name_column,
@ -1174,12 +1184,14 @@ int obi_ecopcr(const char* i_dms_name,
o_strand_column,
o_primer1_column, o_primer2_column,
o_error1_column, o_error2_column,
o_temp1_column, o_temp2_column) < 0)
o_temp1_column, o_temp2_column);
if (printed < 0)
{
obidebug(1, "\nError writing the ecopcr result");
return -1;
}
o_idx++;
else if (printed > 0)
o_idx++;
}
}
}
@ -1220,7 +1232,7 @@ int obi_ecopcr(const char* i_dms_name,
return -1;
}
fprintf(stderr,"\rDone : 100 %% ");
fprintf(stderr,"\rDone : 100 %% \n");
return 0;
return 0;

View File

@ -81,8 +81,8 @@
* @param o_dms_name The path to the output DMS.
* @param o_view_name The name of the output view.
* @param o_view_comments The comments to associate with the output view.
* @param primer1 The first primer.
* @param primer2 The second primer.
* @param primer1 The first primer, length must be less than or equal to 32 (because of apat lib limitation).
* @param primer2 The second primer, length must be less than or equal to 32 (because of apat lib limitation).
* @param error_max The maximum number of errors allowed per primer for amplification.
* @param min_len The minimum length of an amplicon.
* @param max_len The maximum length of an amplicon.
@ -93,8 +93,8 @@
* @param salt_concentration The salt concentration used for estimating the Tm.
* @param salt_correction_method The method used for estimating the Tm (melting temperature) between the primers and their corresponding
* target sequences. SANTALUCIA: 1, or OWCZARZY: 2.
* @param keep_nucleotides The number of nucleotides to keep on each side of the in silico amplified sequences
* (already including the amplified DNA fragment plus the two target sequences of the primers).
* @param keep_nucleotides The number of nucleotides to keep on each side of the in silico amplified sequences, not including primers (primers automatically entirely kept if > 0).
* @param keep_primers Whether primers are kept attached to the output sequences.
* @param kingdom_mode Whether the kingdom or the superkingdom informations should be printed to the output.
*
* @returns A value indicating the success of the operation.
@ -121,6 +121,7 @@ int obi_ecopcr(const char* i_dms_name,
double salt_concentration,
int salt_correction_method,
int keep_nucleotides,
bool keep_primers,
bool kingdom_mode);
#endif /* OBI_ECOPCR_H_ */

View File

@ -71,9 +71,12 @@ static int create_output_columns(Obiview_p o_view);
* @param name The assigned scientific name.
* @param assigned_status_column A pointer on the column where the assigned status should be written.
* @param assigned The assigned status (whether the sequence was assigned to a taxon or not).
* @param best_match_column A pointer on the column where the list of ids of the best matches should be written.
* @param best_match_ids_column A pointer on the column where the list of ids of the best matches should be written.
* @param best_match_ids The list of ids of the best matches as an array of the concatenated ids separated by '\0'.
* @param best_match_ids_length The total length of the array of ids of best matches.
* @param best_match_taxids_column A pointer on the column where the list of taxids of the best matches should be written.
* @param best_match_taxids The list of taxids of the best matches as an array of the taxids.
* @param best_match_taxids_length The length of the array of taxids of best matches.
* @param score_column A pointer on the column where the score should be written.
* @param score The similarity score of the sequence with its best match(es).
*
@ -87,7 +90,8 @@ int print_assignment_result(Obiview_p output_view, index_t line,
OBIDMS_column_p assigned_taxid_column, int32_t taxid,
OBIDMS_column_p assigned_name_column, const char* name,
OBIDMS_column_p assigned_status_column, bool assigned,
OBIDMS_column_p best_match_column, const char* best_match_ids, int best_match_ids_length,
OBIDMS_column_p best_match_ids_column, const char* best_match_ids, int best_match_ids_length,
OBIDMS_column_p best_match_taxids_column, const int32_t* best_match_taxids, int best_match_taxids_length,
OBIDMS_column_p score_column, double score);
@ -100,37 +104,44 @@ int print_assignment_result(Obiview_p output_view, index_t line,
static int create_output_columns(Obiview_p o_view)
{
// Score column
if (obi_view_add_column(o_view, ECOTAG_SCORE_COLUMN_NAME, -1, NULL, OBI_FLOAT, 0, 1, NULL, false, false, false, NULL, NULL, -1, ECOTAG_SCORE_COLUMN_NAME, true) < 0)
if (obi_view_add_column(o_view, ECOTAG_SCORE_COLUMN_NAME, -1, NULL, OBI_FLOAT, 0, 1, NULL, false, false, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the score in ecotag");
return -1;
}
// Assigned taxid column
if (obi_view_add_column(o_view, ECOTAG_TAXID_COLUMN_NAME, -1, NULL, OBI_INT, 0, 1, NULL, false, false, false, NULL, NULL, -1, ECOTAG_TAXID_COLUMN_NAME, true) < 0)
if (obi_view_add_column(o_view, ECOTAG_TAXID_COLUMN_NAME, -1, NULL, OBI_INT, 0, 1, NULL, false, false, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the assigned taxid in ecotag");
return -1;
}
// Assigned scientific name column
if (obi_view_add_column(o_view, ECOTAG_NAME_COLUMN_NAME, -1, NULL, OBI_STR, 0, 1, NULL, false, false, false, NULL, NULL, -1, ECOTAG_NAME_COLUMN_NAME, true) < 0)
if (obi_view_add_column(o_view, ECOTAG_NAME_COLUMN_NAME, -1, NULL, OBI_STR, 0, 1, NULL, false, false, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the assigned scientific name in ecotag");
return -1;
}
// Assignement status column
if (obi_view_add_column(o_view, ECOTAG_STATUS_COLUMN_NAME, -1, NULL, OBI_BOOL, 0, 1, NULL, false, false, false, NULL, NULL, -1, ECOTAG_STATUS_COLUMN_NAME, true) < 0)
if (obi_view_add_column(o_view, ECOTAG_STATUS_COLUMN_NAME, -1, NULL, OBI_BOOL, 0, 1, NULL, false, false, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the assignment status in ecotag");
return -1;
}
// Column for array of best match ids
if (obi_view_add_column(o_view, ECOTAG_BEST_MATCH_IDS_COLUMN_NAME, -1, NULL, OBI_STR, 0, 1, NULL, false, true, false, NULL, NULL, -1, ECOTAG_BEST_MATCH_IDS_COLUMN_NAME, true) < 0)
if (obi_view_add_column(o_view, ECOTAG_BEST_MATCH_IDS_COLUMN_NAME, -1, NULL, OBI_STR, 0, 1, NULL, false, true, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the array of ids of the best match in ecotag");
obidebug(1, "\nError creating the column for the array of ids of best matches in ecotag");
return -1;
}
// Column for array of best match taxids
if (obi_view_add_column(o_view, ECOTAG_BEST_MATCH_TAXIDS_COLUMN_NAME, -1, NULL, OBI_INT, 0, 1, NULL, false, true, false, NULL, NULL, -1, "{}", true) < 0)
{
obidebug(1, "\nError creating the column for the array of taxids of best matches in ecotag");
return -1;
}
@ -142,7 +153,8 @@ int print_assignment_result(Obiview_p output_view, index_t line,
OBIDMS_column_p assigned_taxid_column, int32_t taxid,
OBIDMS_column_p assigned_name_column, const char* name,
OBIDMS_column_p assigned_status_column, bool assigned,
OBIDMS_column_p best_match_column, const char* best_match_ids, int best_match_ids_length,
OBIDMS_column_p best_match_ids_column, const char* best_match_ids, int best_match_ids_length,
OBIDMS_column_p best_match_taxids_column, const int32_t* best_match_taxids, int best_match_taxids_length,
OBIDMS_column_p score_column, double score)
{
// Write the assigned taxid
@ -167,9 +179,16 @@ int print_assignment_result(Obiview_p output_view, index_t line,
}
// Write the best match ids
if (obi_set_array_with_col_p_in_view(output_view, best_match_column, line, best_match_ids, (uint8_t)(sizeof(char)*8), best_match_ids_length) < 0)
if (obi_set_array_with_col_p_in_view(output_view, best_match_ids_column, line, best_match_ids, (uint8_t)(sizeof(char)*8), best_match_ids_length) < 0)
{
obidebug(1, "\nError writing a assignment status in a column when writing ecotag results");
obidebug(1, "\nError writing the array of best match ids in a column when writing ecotag results");
return -1;
}
// Write the best match taxids
if (obi_set_array_with_col_p_in_view(output_view, best_match_taxids_column, line, best_match_taxids, (uint8_t)(sizeof(OBI_INT)*8), best_match_taxids_length) < 0)
{
obidebug(1, "\nError writing the array of best match taxids in a column when writing ecotag results");
return -1;
}
@ -199,7 +218,8 @@ int obi_ecotag(const char* dms_name,
const char* taxonomy_name,
const char* output_view_name,
const char* output_view_comments,
double ecotag_threshold) // TODO different threshold for the similarity sphere around ref seqs
double ecotag_threshold,
double bubble_threshold)
{
// For each sequence
@ -220,6 +240,7 @@ int obi_ecotag(const char* dms_name,
index_t query_seq_idx, ref_seq_idx;
double score, best_score;
double threshold;
double lca_threshold;
int lcs_length;
int ali_length;
Kmer_table_p ktable;
@ -235,6 +256,8 @@ int obi_ecotag(const char* dms_name,
char* best_match_ids;
char* best_match_ids_to_store;
int32_t best_match_ids_length;
int32_t* best_match_taxids;
int32_t* best_match_taxids_to_store;
int best_match_count;
int buffer_size;
int best_match_ids_buffer_size;
@ -263,7 +286,8 @@ int obi_ecotag(const char* dms_name,
OBIDMS_column_p assigned_taxid_column = NULL;
OBIDMS_column_p assigned_name_column = NULL;
OBIDMS_column_p assigned_status_column = NULL;
OBIDMS_column_p best_match_column = NULL;
OBIDMS_column_p best_match_ids_column = NULL;
OBIDMS_column_p best_match_taxids_column = NULL;
OBIDMS_column_p lca_taxid_a_column = NULL;
OBIDMS_column_p score_a_column = NULL;
OBIDMS_column_p ref_taxid_column = NULL;
@ -367,10 +391,10 @@ int obi_ecotag(const char* dms_name,
return -1;
}
free(db_threshold_str);
if (ecotag_threshold < db_threshold)
if (bubble_threshold < db_threshold)
{
fprintf(stderr, "\nError: The threshold demanded (%f) is lower than the threshold used to build the reference database (%f).\n\n",
ecotag_threshold, db_threshold);
bubble_threshold, db_threshold);
return -1;
}
@ -396,7 +420,8 @@ int obi_ecotag(const char* dms_name,
assigned_taxid_column = obi_view_get_column(output_view, ECOTAG_TAXID_COLUMN_NAME);
assigned_name_column = obi_view_get_column(output_view, ECOTAG_NAME_COLUMN_NAME);
assigned_status_column = obi_view_get_column(output_view, ECOTAG_STATUS_COLUMN_NAME);
best_match_column = obi_view_get_column(output_view, ECOTAG_BEST_MATCH_IDS_COLUMN_NAME);
best_match_ids_column = obi_view_get_column(output_view, ECOTAG_BEST_MATCH_IDS_COLUMN_NAME);
best_match_taxids_column = obi_view_get_column(output_view, ECOTAG_BEST_MATCH_TAXIDS_COLUMN_NAME);
score_column = obi_view_get_column(output_view, ECOTAG_SCORE_COLUMN_NAME);
// Open the used reference columns
@ -453,9 +478,17 @@ int obi_ecotag(const char* dms_name,
return -1;
}
best_match_taxids = (int32_t*) malloc(buffer_size* sizeof(int32_t));
if (best_match_taxids == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for the best match taxid array in ecotag");
return -1;
}
for (i=0; i < query_count; i++)
{
if (i%100 == 0)
if (i%1000 == 0)
fprintf(stderr,"\rDone : %f %% ", (i / (float) query_count)*100);
best_match_count = 0;
@ -514,7 +547,7 @@ int obi_ecotag(const char* dms_name,
// Store in best match array
// Grow match array if needed
// Grow match and taxid array if needed
if (best_match_count == buffer_size)
{
buffer_size = buffer_size*2;
@ -525,6 +558,13 @@ int obi_ecotag(const char* dms_name,
obidebug(1, "\nError reallocating match array when assigning");
return -1;
}
best_match_taxids = (int32_t*) realloc(best_match_taxids, buffer_size*sizeof(int32_t));
if (best_match_taxids == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError reallocating match taxids array when assigning");
return -1;
}
}
id = obi_get_str_with_elt_idx_and_col_p_in_view(ref_view, ref_id_column, j, 0);
@ -545,6 +585,7 @@ 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++;
strcpy(best_match_ids+best_match_ids_length, id);
best_match_ids_length = best_match_ids_length + id_len + 1;
@ -558,11 +599,16 @@ int obi_ecotag(const char* dms_name,
{
best_match_idx = best_match_array[j];
// Find the LCA for the chosen threshold
// Find the LCA for the highest threshold between best_score and the chosen bubble threshold
score_array = obi_get_array_with_col_p_in_view(ref_view, score_a_column, best_match_idx, &lca_array_length);
if (bubble_threshold < best_score)
lca_threshold = best_score;
else
lca_threshold = bubble_threshold;
k = 0;
while ((k < lca_array_length) && (score_array[k] >= ecotag_threshold))
while ((k < lca_array_length) && (score_array[k] >= lca_threshold))
k++;
if (k>0)
@ -570,12 +616,12 @@ int obi_ecotag(const char* dms_name,
lca_array = obi_get_array_with_col_p_in_view(ref_view, lca_taxid_a_column, best_match_idx, &lca_array_length);
if (j>0)
{
lca = obi_taxo_get_taxon_with_taxid(taxonomy, lca_taxid);
if (lca == NULL)
{
obidebug(1, "\nError getting a taxon from a taxid when doing taxonomic assignment");
return -1;
}
// lca = obi_taxo_get_taxon_with_taxid(taxonomy, lca_taxid);
// if (lca == NULL)
// {
// obidebug(1, "\nError getting a taxon from a taxid when doing taxonomic assignment");
// return -1;
// }
lca_in_array = obi_taxo_get_taxon_with_taxid(taxonomy, lca_array[k-1]);
if (lca_in_array == NULL)
{
@ -629,6 +675,7 @@ int obi_ecotag(const char* dms_name,
else
lca_name = lca->name;
best_match_ids_to_store = best_match_ids;
best_match_taxids_to_store = best_match_taxids;
}
else
{
@ -636,6 +683,7 @@ int obi_ecotag(const char* dms_name,
lca_name = OBIStr_NA;
lca_taxid = OBIInt_NA;
best_match_ids_to_store = OBITuple_NA;
best_match_taxids_to_store = OBITuple_NA;
score = OBIFloat_NA;
}
@ -644,7 +692,8 @@ int obi_ecotag(const char* dms_name,
assigned_taxid_column, lca_taxid,
assigned_name_column, lca_name,
assigned_status_column, assigned,
best_match_column, best_match_ids_to_store, best_match_ids_length,
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,
score_column, best_score
) < 0)
return -1;
@ -652,6 +701,7 @@ int obi_ecotag(const char* dms_name,
free(best_match_array);
free(best_match_ids);
free(best_match_taxids);
obi_close_taxonomy(taxonomy);
obi_save_and_close_view(query_view);

View File

@ -23,7 +23,8 @@
#define ECOTAG_TAXID_COLUMN_NAME "TAXID"
#define ECOTAG_NAME_COLUMN_NAME "SCIENTIFIC_NAME"
#define ECOTAG_STATUS_COLUMN_NAME "ID_STATUS"
#define ECOTAG_BEST_MATCH_IDS_COLUMN_NAME "BEST_MATCH"
#define ECOTAG_BEST_MATCH_IDS_COLUMN_NAME "BEST_MATCH_IDS"
#define ECOTAG_BEST_MATCH_TAXIDS_COLUMN_NAME "BEST_MATCH_TAXIDS"
#define ECOTAG_SCORE_COLUMN_NAME "BEST_IDENTITY"
@ -41,12 +42,14 @@
* @param output_view_name The name to give to the output view.
* @param output_view_comments The comments to associate to the output view.
* @param ecotag_threshold The threshold at which to assign.
* @param bubble_threshold The threshold at which to look for an LCA (i.e. minimum identity considered for the assignment circle);
* the threshold actually used will be the highest between this value and the best assignment score found.
*
* The algorithm works like this:
* For each query sequence:
* Align with reference database
* Keep the indices of all the best matches
* For each kept index, get the LCA at that threshold as stored in the reference database, then the LCA of those LCAs
* For each kept index, get the LCA at the highest threshold between bubble_threshold and the best assignment score found (as stored in the reference database), then the LCA of those LCAs
* Write result (max score, threshold, taxid and scientific name of the LCA assigned, list of the ids of the best matches)
*
* @returns A value indicating the success of the operation.
@ -64,7 +67,8 @@ int obi_ecotag(const char* dms_name,
const char* taxonomy_name,
const char* output_view_name,
const char* output_view_comments,
double ecotag_threshold);
double ecotag_threshold,
double bubble_threshold);
#endif /* OBI_ECOTAG_H_ */

View File

@ -648,7 +648,7 @@ int truncate_avl_data_to_size_used(OBIDMS_avl_data_p avl_data) // TODO is it nec
new_data_size = ((index_t) multiple) * getpagesize();
// Check that it is actually greater than the current size of the file, otherwise no need to truncate
if ((avl_data->header)->data_size_max == new_data_size)
if ((avl_data->header)->data_size_max >= new_data_size)
return 0;
// Get the file descriptor
@ -667,7 +667,7 @@ int truncate_avl_data_to_size_used(OBIDMS_avl_data_p avl_data) // TODO is it nec
if (ftruncate(file_descriptor, file_size) < 0)
{
obi_set_errno(OBI_AVL_ERROR);
obidebug(1, "\nError truncating an AVL data file");
obidebug(1, "\nError truncating an AVL data file, old data size = %lld, new data size = %lld", (avl_data->header)->data_size_max, new_data_size);
return -1;
}
@ -2259,7 +2259,13 @@ index_t obi_avl_add(OBIDMS_avl_p avl, Obi_blob_p value)
parent = next;
// Compare the crc of the value with the crc of the current node
comp = (current_node->crc64) - crc;
//comp = (current_node->crc64) - crc;
if ((current_node->crc64) == crc)
comp = 0;
else if ((current_node->crc64) > crc)
comp = 1;
else
comp = -1;
if (comp == 0)
{ // check if really same value
@ -2354,7 +2360,13 @@ index_t obi_avl_find(OBIDMS_avl_p avl, Obi_blob_p value)
current_node = (avl->tree)+next;
// Compare the crc of the value with the crc of the current node
comp = (current_node->crc64) - crc;
//comp = (current_node->crc64) - crc;
if ((current_node->crc64) == crc)
comp = 0;
else if ((current_node->crc64) > crc)
comp = 1;
else
comp = -1;
if (comp == 0)
{ // Check if really same value

View File

@ -696,6 +696,12 @@ int obi_clean_dms(const char* dms_path)
// return -1;
// }
if (obi_close_dms(dms, true) < 0)
{
obidebug(1, "\nError closing a DMS after cleaning");
return -1;
}
return 0;
}
@ -1403,6 +1409,107 @@ DIR* opendir_in_dms(OBIDMS_p dms, const char* path_name)
}
char* obi_dms_formatted_infos(OBIDMS_p dms, bool detailed)
{
char* dms_infos = NULL;
char* view_infos = NULL;
char* view_name = NULL;
char* tax_name = NULL;
char* all_tax_dir_path = NULL;
int i;
struct dirent* dp;
Obiview_p view;
// DMS name
dms_infos = (char*) malloc((strlen("# DMS name: ")+strlen(dms->dms_name)+strlen("\n# Views:\n")+1) * sizeof(char));
if (dms_infos == NULL)
{
obidebug(1, "\nError allocating memory for DMS formatted infos");
return NULL;
}
strcpy(dms_infos, "# DMS name: ");
strcat(dms_infos, dms->dms_name);
strcat(dms_infos, "\n# Views:\n");
// Go through views and get their infos
rewinddir(dms->view_directory);
while ((dp = readdir(dms->view_directory)) != NULL)
{
if ((dp->d_name)[0] == '.')
continue;
i=0;
while ((dp->d_name)[i] != '.')
i++;
view_name = (char*) malloc((i+1) * sizeof(char));
if (view_name == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for a view name when getting formatted DMS infos: file %s", dp->d_name);
return NULL;
}
strncpy(view_name, dp->d_name, i);
view_name[i] = '\0';
view = obi_open_view(dms, view_name);
if (view == NULL)
{
obidebug(1, "\nError opening a view to get DMS formatted infos");
return NULL;
}
if (detailed)
view_infos = obi_view_formatted_infos(view, detailed);
else
view_infos = obi_view_formatted_infos_one_line(view);
if (view_infos == NULL)
{
obidebug(1, "\nError getting a view infos to get DMS formatted infos");
return NULL;
}
dms_infos = realloc(dms_infos, (strlen(dms_infos)+strlen(view_infos)+1) * sizeof(char));
if (dms_infos == NULL)
{
obidebug(1, "\nError reallocating memory for DMS formatted infos");
return NULL;
}
strcat(dms_infos, view_infos);
if (obi_save_and_close_view(view) < 0)
{
obidebug(1, "\nError closing view while getting DMS formatted infos");
return NULL;
}
if (detailed)
{
dms_infos = realloc(dms_infos, (strlen(dms_infos)+2) * sizeof(char));
strcat(dms_infos, "\n");
}
}
// Add taxonomies
dms_infos = realloc(dms_infos, (strlen(dms_infos)+strlen("\n# Taxonomies:\n")+1) * sizeof(char));
if (dms_infos == NULL)
{
obidebug(1, "\nError reallocating memory for DMS formatted infos");
return NULL;
}
strcat(dms_infos, "# Taxonomies:\n");
rewinddir(dms->tax_directory);
while ((dp = readdir(dms->tax_directory)) != NULL)
{
if ((dp->d_name)[0] == '.')
continue;
tax_name = dp->d_name;
dms_infos = realloc(dms_infos, (strlen(dms_infos)+strlen(" # ")+strlen(view_infos)+1) * sizeof(char));
if (dms_infos == NULL)
{
obidebug(1, "\nError reallocating memory for DMS formatted infos");
return NULL;
}
strcat(dms_infos, " # ");
strcat(dms_infos, tax_name);
}
return dms_infos;
}
// TODO move somewhere else maybe
// TODO discuss arguments
obiversion_t obi_import_column(const char* dms_path_1, const char* dms_path_2, const char* column_name, obiversion_t version_number)
@ -1490,7 +1597,7 @@ obiversion_t obi_import_column(const char* dms_path_1, const char* dms_path_2, c
memcpy(column_2->data, column_1->data, header_1->data_size);
// Copy the AVL files if there are some (overwriting the automatically created files)
if ((header_1->returned_data_type == OBI_STR) || (header_1->returned_data_type == OBI_SEQ) || (header_1->returned_data_type == OBI_QUAL))
if ((header_1->tuples) || ((header_1->returned_data_type == OBI_STR) || (header_1->returned_data_type == OBI_SEQ) || (header_1->returned_data_type == OBI_QUAL)))
{
avl_name_1 = (char*) malloc((strlen(header_1->indexer_name) + 1) * sizeof(char));
if (avl_name_1 == NULL)
@ -1653,6 +1760,12 @@ int obi_import_view(const char* dms_path_1, const char* dms_path_2, const char*
else // Non-typed view
view_2 = obi_new_view(dms_2, view_name_2, NULL, NULL, (view_1->infos)->comments);
if (view_2 == NULL)
{
obidebug(1, "\nError creating the new view to import a view in a DMS");
return -1;
}
// Import line count
view_2->infos->line_count = view_1->infos->line_count;

View File

@ -459,6 +459,23 @@ char* obi_dms_get_full_path(OBIDMS_p dms, const char* path_name);
DIR* opendir_in_dms(OBIDMS_p dms, const char* path_name);
/**
* @brief Returns the informations of a DMS with a human readable format (dms name, taxonomies and view infos).
*
* @warning The returned pointer has to be freed by the caller.
*
* @param column A pointer on a DMS.
* @param detailed Whether the informations should contain detailed view infos.
*
* @returns A pointer on a character array where the formatted DMS informations are stored.
* @retval NULL if an error occurred.
*
* @since September 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_dms_formatted_infos(OBIDMS_p dms, bool detailed);
/**
* @brief Imports a column, copying it from a DMS to another DMS, and returns the version of the column in the destination DMS.
*

View File

@ -2376,9 +2376,10 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
// and the deleted taxids with no current reference. An element of the list is composed of the taxid, and the index
// of the taxon in the taxa structure, or -1 for deleted taxids.
// Creating the merged list requires to merge the 3 ordered lists into one.
while (((nT < (tax->taxa)->count) && ((tax->taxa)->taxon[nT].taxid < old_taxid)) || ((nD >= 0) && (delnodes[nD] < old_taxid)))
while (((nT < (tax->taxa)->count) && ((tax->taxa)->taxon[nT].taxid < old_taxid)) ||
((nD >= 0) && (delnodes[nD] < old_taxid)))
{
if ((tax->taxa)->taxon[nT].taxid < delnodes[nD])
if ((nT < (tax->taxa)->count) && (tax->taxa)->taxon[nT].taxid < delnodes[nD])
{ // Add element from taxa list
// Enlarge structure if needed
if (n == buffer_size)
@ -2401,7 +2402,7 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
nT++;
n++;
}
else if (delnodes[nD] < (tax->taxa)->taxon[nT].taxid)
else
{ // Add element from deleted taxids list
// Enlarge structure if needed
if (n == buffer_size)
@ -3036,12 +3037,12 @@ OBIDMS_taxonomy_p obi_read_taxonomy(OBIDMS_p dms, const char* taxonomy_name, boo
strcpy(tax->tax_name, taxonomy_name);
buffer_size = 2048;
taxonomy_path = get_taxonomy_path(dms, taxonomy_name);
if (taxonomy_path == NULL)
return NULL;
buffer_size = strlen(taxonomy_path) + strlen(taxonomy_name) + 6;
// Read ranks
ranks_file_name = (char*) malloc(buffer_size*sizeof(char));
if (ranks_file_name == NULL)
@ -3648,6 +3649,18 @@ ecotx_t* obi_taxo_get_taxon_with_taxid(OBIDMS_taxonomy_p taxonomy, int32_t taxid
}
char* obi_taxo_get_name_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx)
{
return (((taxonomy->names)->names)[idx]).name;
}
ecotx_t* obi_taxo_get_taxon_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx)
{
return (((taxonomy->names)->names)[idx]).taxon;
}
int obi_taxo_is_taxon_under_taxid(ecotx_t* taxon, int32_t other_taxid) // TODO discuss that this doesn't work with deprecated taxids
{
ecotx_t* next_parent;

View File

@ -447,8 +447,51 @@ ecotx_t* obi_taxo_get_superkingdom(ecotx_t* taxon, OBIDMS_taxonomy_p taxonomy);
const char* obi_taxo_rank_index_to_label(int32_t rank_idx, ecorankidx_t* ranks);
// TODO
/**
* @brief Function checking whether a taxid is included in a subset of the taxonomy.
*
* @param taxonomy A pointer on the taxonomy structure.
* @param restrict_to_taxids An array of taxids. The researched taxid must be under at least one of those array taxids.
* @param count Number of taxids in restrict_to_taxids.
* @param taxid The taxid to check.
*
* @returns A value indicating whether the taxid is included in the chosen subset of the taxonomy.
* @retval 0 if the taxid is not included in the subset of the taxonomy.
* @retval 1 if the taxid is included in the subset of the taxonomy.
*
* @since October 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
int obi_taxo_is_taxid_included(OBIDMS_taxonomy_p taxonomy,
int32_t* restrict_to_taxids,
int32_t count,
int32_t taxid);
/**
* @brief Function returning the name of a taxon from its index in the taxonomy name index (econameidx_t).
*
* @param taxonomy A pointer on the taxonomy structure.
* @param idx The index at which the name is in the taxonomy name index (econameidx_t).
*
* @returns The taxon name.
*
* @since October 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_taxo_get_name_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx);
/**
* @brief Function returning a taxon structure from its index in the taxonomy name index (econameidx_t).
*
* @param taxonomy A pointer on the taxonomy structure.
* @param idx The index at which the taxon is in the taxonomy name index (econameidx_t).
*
* @returns The taxon structure.
*
* @since October 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
ecotx_t* obi_taxo_get_taxon_from_name_idx(OBIDMS_taxonomy_p taxonomy, int32_t idx);

View File

@ -1312,19 +1312,10 @@ OBIDMS_column_p obi_create_column(OBIDMS_p dms,
return NULL;
}
// Store the associated column reference if needed // TODO discuss cases
if (data_type == OBI_QUAL)
// Store the associated column reference if needed
if ((associated_column_name != NULL) && (*associated_column_name != '\0'))
{
if ((associated_column_name == NULL) || (*associated_column_name == '\0'))
{
obidebug(1, "\nError: The name of the associated column when creating a new column is NULL");
munmap(new_column->header, header_size);
close(column_file_descriptor);
free(new_column);
return NULL;
}
strcpy((header->associated_column).column_name, associated_column_name);
if (associated_column_version == -1)
{
obidebug(1, "\nError: The version of the associated column when creating a new column is not defined");
@ -1336,6 +1327,7 @@ OBIDMS_column_p obi_create_column(OBIDMS_p dms,
(header->associated_column).version = associated_column_version;
}
// If the data type is OBI_STR, OBI_SEQ or OBI_QUAL, the associated obi_indexer is opened or created
if ((returned_data_type == OBI_STR) || (returned_data_type == OBI_SEQ) || (returned_data_type == OBI_QUAL) || tuples)
{
@ -1350,6 +1342,8 @@ OBIDMS_column_p obi_create_column(OBIDMS_p dms,
}
strncpy(header->indexer_name, final_indexer_name, INDEXER_MAX_NAME);
}
else
new_column->indexer = NULL;
// Fill the data with NA values
obi_ini_to_NA_values(new_column, 0, nb_lines);
@ -1558,6 +1552,8 @@ OBIDMS_column_p obi_open_column(OBIDMS_p dms,
return NULL;
}
}
else
column->indexer = NULL;
if (close(column_file_descriptor) < 0)
{
@ -1693,8 +1689,8 @@ int obi_close_column(OBIDMS_column_p column)
if (obi_dms_unlist_column(column->dms, column) < 0)
ret_val = -1;
// If the data type is OBI_STR, OBI_SEQ or OBI_QUAL, the associated indexer is closed
if (((column->header)->returned_data_type == OBI_STR) || ((column->header)->returned_data_type == OBI_SEQ) || ((column->header)->returned_data_type == OBI_QUAL))
// If it's a tuple column or the data type is OBI_STR, OBI_SEQ or OBI_QUAL, the associated indexer is closed
if ((column->indexer) != NULL)
if (obi_close_indexer(column->indexer) < 0)
ret_val = -1;
@ -1729,16 +1725,32 @@ int obi_close_column(OBIDMS_column_p column)
int obi_clone_column_indexer(OBIDMS_column_p column)
{
char* new_indexer_name;
int i;
new_indexer_name = obi_build_indexer_name((column->header)->name, (column->header)->version);
if (new_indexer_name == NULL)
return -1;
column->indexer = obi_clone_indexer(column->indexer, new_indexer_name); // TODO Need to lock this somehow?
if (column->indexer == NULL)
i=0;
while (true) // find avl name not already used
{
obidebug(1, "\nError cloning a column's indexer to make it writable");
return -1;
new_indexer_name = obi_build_indexer_name((column->header)->name, ((column->header)->version)+i);
if (new_indexer_name == NULL)
return -1;
column->indexer = obi_clone_indexer(column->indexer, new_indexer_name); // TODO Need to lock this somehow?
if (column->indexer == NULL)
{
if (errno == EEXIST)
{
free(new_indexer_name);
i++;
}
else
{
free(new_indexer_name);
obidebug(1, "\nError cloning a column's indexer to make it writable");
return -1;
}
}
else
break;
}
strcpy((column->header)->indexer_name, new_indexer_name);
@ -1973,7 +1985,11 @@ int obi_enlarge_column(OBIDMS_column_p column)
// Calculate the new file size
old_line_count = (column->header)->line_count;
new_line_count = old_line_count * COLUMN_GROWTH_FACTOR;
new_line_count = ceil((double) old_line_count * (double) COLUMN_GROWTH_FACTOR);
if (new_line_count > old_line_count+100000)
new_line_count = old_line_count+100000;
else if (new_line_count < old_line_count+1000)
new_line_count = old_line_count+1000;
if (new_line_count > MAXIMUM_LINE_COUNT)
{
@ -2381,6 +2397,118 @@ char* obi_get_elements_names(OBIDMS_column_p column)
}
char* obi_get_formatted_elements_names(OBIDMS_column_p column)
{
char* elements_names;
int i, j;
int elt_idx;
int len;
elements_names = (char*) malloc(((column->header)->elements_names_length + (column->header)->nb_elements_per_line) * sizeof(char));
if (elements_names == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for elements names");
return NULL;
}
j = 0;
for (i=0; i < (column->header)->nb_elements_per_line; i++)
{
elt_idx = ((column->header)->elements_names_idx)[i];
len = strlen(((column->header)->elements_names)+elt_idx);
memcpy(elements_names+j, ((column->header)->elements_names)+elt_idx, len*sizeof(char));
j = j + len;
elements_names[j] = ';';
j++;
elements_names[j] = ' ';
j++;
}
elements_names[j - 1] = '\0';
return elements_names;
}
char* obi_column_formatted_infos(OBIDMS_column_p column, bool detailed)
{
char* column_infos = NULL;
char* elt_names = NULL;
char* data_type_str = NULL;
char* comments = NULL;
// Get element names informations
elt_names = obi_get_formatted_elements_names(column);
if (elt_names == NULL)
{
obidebug(1, "\nError getting formatted elements names for formatted columns infos");
return NULL;
}
// Get data type informations
data_type_str = name_data_type((column->header)->returned_data_type);
if (data_type_str == NULL)
{
obidebug(1, "\nError getting formatted data type for formatted columns infos");
return NULL;
}
// Get commments if detailed informations required
if (detailed)
comments = (column->header)->comments;
// Build the string of formatted infos, allocating memory as needed
// Data type
column_infos = (char*) malloc((strlen("data type: ")+strlen(data_type_str)+1) * sizeof(char));
if (column_infos == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for formatted column infos");
return NULL;
}
strcpy(column_infos, "data type: ");
strcat(column_infos, data_type_str);
// Element names if more than 1
if ((column->header)->nb_elements_per_line > 1)
{
column_infos = realloc(column_infos, (strlen(column_infos)+strlen(", elements: ")+strlen(elt_names)+1) * sizeof(char));
if (column_infos == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for formatted column infos");
return NULL;
}
strcat(column_infos, ", elements: ");
strcat(column_infos, elt_names);
}
if (detailed && (strlen(comments)>2)) // Add all comments if required and not empty
{
column_infos = realloc(column_infos, (strlen(column_infos)+strlen("\nComments:\n")+strlen(comments)+1) * sizeof(char));
if (column_infos == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for formatted column infos");
return NULL;
}
strcat(column_infos, "\nComments:\n");
strcat(column_infos, comments);
}
// "data type: OBI_TYPE, element names: [formatted element names](, all comments)"
free(elt_names);
free(data_type_str);
return column_infos;
}
int obi_column_prepare_to_set_value(OBIDMS_column_p column, index_t line_nb, index_t elt_idx)
{
@ -2424,7 +2552,6 @@ int obi_column_prepare_to_set_value(OBIDMS_column_p column, index_t line_nb, ind
}
int obi_column_prepare_to_get_value(OBIDMS_column_p column, index_t line_nb)
{
if ((line_nb+1) > ((column->header)->line_count))

View File

@ -505,6 +505,39 @@ index_t obi_column_get_element_index_from_name(OBIDMS_column_p column, const cha
char* obi_get_elements_names(OBIDMS_column_p column);
/**
* @brief Recovers the elements names of the lines of a column with a human readable format ("0; 1; 2; ...; n\0").
*
* @warning The returned pointer has to be freed by the caller.
*
* @param column A pointer on an OBIDMS column.
*
* @returns A pointer on a character array where the elements names are stored.
* @retval NULL if an error occurred.
*
* @since September 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_get_formatted_elements_names(OBIDMS_column_p column);
/**
* @brief Returns the informations of a column with a human readable format (data type, element names, comments).
*
* @warning The returned pointer has to be freed by the caller.
*
* @param column A pointer on an OBIDMS column.
* @param detailed Whether the informations should contain column comments or just data type and element names.
*
* @returns A pointer on a character array where the formatted column informations are stored.
* @retval NULL if an error occurred.
*
* @since September 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_column_formatted_infos(OBIDMS_column_p column, bool detailed);
/**
* @brief Prepares a column to set a value.
*

View File

@ -34,8 +34,8 @@
* @brief enum for the boolean OBIType.
*/
typedef enum OBIBool {
FALSE = 0,
TRUE = 1,
OBIFalse = 0,
OBITrue = 1,
OBIBool_NA = 2
} obibool_t, *obibool_p; /**< a boolean true/false value */ // TODO check name convention?

View File

@ -17,6 +17,7 @@
#include <sys/mman.h>
#include <inttypes.h>
#include <math.h>
#include <time.h>
//#include <ctype.h>
#include "obiview.h"
@ -254,11 +255,15 @@ static int update_lines(Obiview_p view, index_t line_count);
/**
* @brief Internal function to clone a column in the context of a view.
*
* Used to edit a closed column.
*
* Clones with the right line selection and replaces the cloned columns with the new ones in the view.
* If there is a line selection, all columns have to be cloned, otherwise only the column of interest is cloned.
*
* @param view A pointer on the view.
* @param column_name The name of the column in the view that should be cloned.
* @param clone_associated Whether to clone the associated column
* (should always be true except when calling from the function itself to avoid infinite recursion).
*
* @returns A pointer on the new column.
* @retval NULL if an error occurred.
@ -266,7 +271,7 @@ static int update_lines(Obiview_p view, index_t line_count);
* @since February 2016
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
static OBIDMS_column_p clone_column_in_view(Obiview_p view, const char* column_name);
static OBIDMS_column_p clone_column_in_view(Obiview_p view, const char* column_name, bool clone_associated);
/**
@ -845,7 +850,7 @@ static int update_lines(Obiview_p view, index_t line_count)
// Clone the column first if needed
if (!(column->writable))
{
column = clone_column_in_view(view, (((view->infos)->column_references)[i]).alias);
column = clone_column_in_view(view, (((view->infos)->column_references)[i]).alias, true);
if (column == NULL)
{
obidebug(1, "\nError cloning a column in a view when updating its line count");
@ -870,12 +875,14 @@ static int update_lines(Obiview_p view, index_t line_count)
}
static OBIDMS_column_p clone_column_in_view(Obiview_p view, const char* column_name)
static OBIDMS_column_p clone_column_in_view(Obiview_p view, const char* column_name, bool clone_associated)
{
int i;
int i, j;
OBIDMS_column_p column = NULL;
OBIDMS_column_p new_column = NULL;
OBIDMS_column_p column_buffer;
OBIDMS_column_p associated_cloned_column = NULL;
char* associated_column_alias = NULL;
// Check that the view is not read-only
if (view->read_only)
@ -916,11 +923,62 @@ static OBIDMS_column_p clone_column_in_view(Obiview_p view, const char* column_n
return NULL;
}
// Look for associated column to clone and reassociate
if ((column_buffer->header->associated_column).column_name[0] != '\0')
{
// Get the associated column alias
j=0;
while (((strcmp((((view->infos)->column_references)[j]).column_refs.column_name, (column_buffer->header->associated_column).column_name)) ||
((((view->infos)->column_references)[j]).column_refs.version != (column_buffer->header->associated_column).version)) &&
j<(view->infos)->column_count) // TODO function for that
j++;
if (j == (view->infos)->column_count) // not found
{
obi_set_errno(OBIVIEW_ERROR);
obidebug(1, "\nCould not find associated column when cloning a column for editing");
return NULL;
}
// No line selection: only this column is cloned, clone and reassociate the associated column
if ((view->line_selection == NULL) && clone_associated)
{
associated_column_alias = (((view->infos)->column_references)[j]).alias;
// Clone the associated column
associated_cloned_column = clone_column_in_view(view, associated_column_alias, false);
// Reassociate both ways
strcpy((associated_cloned_column->header->associated_column).column_name, column->header->name);
(associated_cloned_column->header->associated_column).version = column->header->version;
strcpy((column->header->associated_column).column_name, associated_cloned_column->header->name);
(column->header->associated_column).version = associated_cloned_column->header->version;
}
else
{
// Line selection: all columns are cloned, check if associated column has been cloned previously (it precedes this one in the list) to reassociate
if (j < i)
{
// Get pointer to associated column
associated_cloned_column = *((OBIDMS_column_p*)ll_get(view->columns, j));
if (associated_cloned_column == NULL)
{
obi_set_errno(OBIVIEW_ERROR);
obidebug(1, "\nError getting a column to clone from the linked list of column pointers of a view");
return NULL;
}
// Reassociate both ways
strcpy((associated_cloned_column->header->associated_column).column_name, column->header->name);
(associated_cloned_column->header->associated_column).version = column->header->version;
strcpy((column->header->associated_column).column_name, associated_cloned_column->header->name);
(column->header->associated_column).version = associated_cloned_column->header->version;
}
}
}
// Close old cloned column
obi_close_column(column_buffer);
if (!strcmp((((view->infos)->column_references)[i]).alias, column_name))
// Found the column to return
// Get the column to return
new_column = column;
}
}
@ -1037,8 +1095,9 @@ static int finish_view(Obiview_p view)
return -1;
}
// Add count column if it's a NUC_SEQ_VIEW with no count column // TODO discuss
if ((!strcmp((view->infos)->view_type, VIEW_TYPE_NUC_SEQS)) && (!obi_view_column_exists(view, COUNT_COLUMN)))
// Add count column if it's a NUC_SEQ_VIEW with no count column (and there's no MERGED_sample column) // TODO discuss
if ((!strcmp((view->infos)->view_type, VIEW_TYPE_NUC_SEQS)) && (!obi_view_column_exists(view, COUNT_COLUMN))
&& (!obi_view_column_exists(view, "MERGED_sample"))) // TODO should eventually compute from merged samples?
{
if (obi_create_auto_count_column(view) < 0)
{
@ -1127,6 +1186,7 @@ static int close_view(Obiview_p view)
obidebug(1, "\nError getting a column to close from the linked list of column pointers of a view");
return -1;
}
if (obi_close_column(column) < 0)
{
obidebug(1, "\nError closing a column while closing a view");
@ -1192,7 +1252,7 @@ static int prepare_to_set_value_in_column(Obiview_p view, OBIDMS_column_p* colum
return -1;
}
(*column_pp) = clone_column_in_view(view, column_name);
(*column_pp) = clone_column_in_view(view, column_name, true);
if ((*column_pp) == NULL)
{
obidebug(1, "\nError trying to clone a column to modify it");
@ -1407,7 +1467,7 @@ static char* view_check_qual_match_seqs(Obiview_p view)
// Test that the lengths of the quality and the sequence are equal
if ((size_t)qual_len != strlen(seq))
{
obidebug(1, "\nError checking the predicate for view %s: The sequences and sequence quality arrays match.", (view->infos)->name);
obidebug(1, "\nError checking the predicate for view %s: The sequences and sequence quality arrays match (index %lld, seq=%s, quality length = %d).", (view->infos)->name, j, seq, qual_len);
return NULL;
}
}
@ -1843,6 +1903,7 @@ Obiview_p obi_new_view_nuc_seqs(OBIDMS_p dms, const char* view_name, Obiview_p v
{
Obiview_p view;
OBIDMS_column_p associated_nuc_column;
OBIDMS_column_p associated_qual_column;
int nb_predicates;
if (view_to_clone != NULL)
@ -1895,6 +1956,10 @@ Obiview_p obi_new_view_nuc_seqs(OBIDMS_p dms, const char* view_name, Obiview_p v
obidebug(1, "Error adding an obligatory column in a nucleotide sequences view");
return NULL;
}
// Associating both ways: associating nuc sequences column to quality column
associated_qual_column = obi_view_get_column(view, QUALITY_COLUMN);
strcpy((associated_nuc_column->header->associated_column).column_name, associated_qual_column->header->name);
(associated_nuc_column->header->associated_column).version = associated_qual_column->header->version;
}
}
@ -1921,7 +1986,7 @@ Obiview_p obi_new_view_nuc_seqs(OBIDMS_p dms, const char* view_name, Obiview_p v
(view->predicate_functions)[(view->nb_predicates)] = view_has_nuc_sequence_column;
(view->predicate_functions)[(view->nb_predicates) + 1] = view_has_id_column;
(view->predicate_functions)[(view->nb_predicates) + 2] = view_has_definition_column;
// if (quality_column) # TODO discuss. Commented bc for example with obi annotate, clone view so clone predicate, then modify seq, so quality is deleted, and predicate boom
// if (quality_column) # TODO fix by triggering predicate deleting if quality deleting. Commented bc for example with obi annotate, clone view so clone predicate, then modify seq, so quality is deleted, and predicate boom
// (view->predicate_functions)[(view->nb_predicates) + 3] = view_has_quality_column;
view->nb_predicates = nb_predicates;
@ -2211,7 +2276,7 @@ Obiview_p obi_open_view(OBIDMS_p dms, const char* view_name)
// TODO return a pointer on the column?
int obi_view_add_column(Obiview_p view,
char* column_name,
char* column_name,
obiversion_t version_number,
const char* alias,
OBIType_t data_type,
@ -2380,11 +2445,12 @@ int obi_view_add_column(Obiview_p view,
}
int obi_view_delete_column(Obiview_p view, const char* column_name)
int obi_view_delete_column(Obiview_p view, const char* column_name, bool delete_file)
{
int i;
bool found;
OBIDMS_column_p column;
char* col_to_delete_path;
// Check that the view is not read-only
if (view->read_only)
@ -2406,8 +2472,31 @@ int obi_view_delete_column(Obiview_p view, const char* column_name)
obidebug(1, "\nError getting a column from the linked list of column pointers of a view when deleting a column from a view");
return -1;
}
// Keep column path if need to delete the file
if (delete_file)
{
col_to_delete_path = obi_column_full_path(view->dms, column->header->name, column->header->version);
if (col_to_delete_path == NULL)
{
obidebug(1, "\nError getting a column file path when deleting a column");
return -1;
}
}
obi_close_column(column);
// Delete file if needed
if (delete_file)
{
if (remove(col_to_delete_path) < 0)
{
obi_set_errno(OBICOL_UNKNOWN_ERROR);
obidebug(1, "\nError deleting a column file when deleting unfinished columns: file %s", col_to_delete_path);
return -1;
}
free(col_to_delete_path);
}
view->columns = ll_delete(view->columns, i);
// TODO how do we check for error? NULL can be empty list
found = true;
@ -2516,6 +2605,144 @@ int obi_view_create_column_alias(Obiview_p view, const char* current_name, const
}
char* obi_view_formatted_infos(Obiview_p view, bool detailed)
{
int i;
char* view_infos = NULL;
char* view_name = NULL;
time_t creation_date;
char* creation_date_str = NULL;
index_t line_count;
char line_count_str[256];
OBIDMS_column_p column;
char* column_alias = NULL;
char* column_infos = NULL;
char* comments = NULL;
// View name
view_name = (view->infos)->name;
view_infos = (char*) malloc((strlen("# View name:\n")+strlen(view_name)+1) * sizeof(char));
strcpy(view_infos, "# View name:\n");
strcat(view_infos, view_name);
// Date created
if (view->read_only) // Date not saved until view is finished writing
{
creation_date = (view->infos)->creation_date;
creation_date_str = ctime(&creation_date);
view_infos = realloc(view_infos, (strlen(view_infos)+strlen("\n# Date created:\n")+strlen(creation_date_str)+1) * sizeof(char));
strcat(view_infos, "\n# Date created:\n");
strcat(view_infos, creation_date_str);
}
// Line count
line_count = (view->infos)->line_count;
snprintf(line_count_str, sizeof line_count_str, "%lld", line_count);
view_infos = realloc(view_infos, (strlen(view_infos)+strlen("\n# Line count:\n")+strlen(line_count_str)+1) * sizeof(char));
strcat(view_infos, "# Line count:\n");
strcat(view_infos, line_count_str);
// Columns: go through each, print their alias then their infos
view_infos = realloc(view_infos, (strlen(view_infos)+strlen("\n# Columns:")+1) * sizeof(char));
strcat(view_infos, "\n# Columns:");
for (i=0; i<((view->infos)->column_count); i++)
{
column = *((OBIDMS_column_p*)ll_get(view->columns, i));
if (column == NULL)
{
obidebug(1, "\nError getting a column from the linked list of column pointers of a view to format view infos");
return NULL;
}
// Column alias
column_alias = (((view->infos)->column_references)[i]).alias;
view_infos = realloc(view_infos, (strlen(view_infos)+strlen("\n")+strlen(column_alias)+strlen(", ")+1) * sizeof(char));
strcat(view_infos, "\n");
strcat(view_infos, column_alias);
strcat(view_infos, ", ");
// Column infos
column_infos = obi_column_formatted_infos(column, detailed);
if (column_infos == NULL)
{
obidebug(1, "\nError getting column infos to format view infos");
return NULL;
}
view_infos = realloc(view_infos, (strlen(view_infos)+strlen(column_infos)+1) * sizeof(char));
strcat(view_infos, column_infos);
free(column_infos);
}
// Get commments if detailed informations required
if (detailed)
{
comments = (view->infos)->comments;
if (strlen(comments)>2) // Add all comments if not empty
{
view_infos = realloc(view_infos, (strlen(view_infos)+strlen("\n# Comments:\n")+strlen(comments)+1) * sizeof(char));
if (view_infos == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for formatted view infos");
return NULL;
}
strcat(view_infos, "\n# Comments:\n");
strcat(view_infos, comments);
}
}
view_infos = realloc(view_infos, (strlen(view_infos)+2) * sizeof(char));
strcat(view_infos, "\n");
return view_infos;
}
char* obi_view_formatted_infos_one_line(Obiview_p view)
{
int i;
char* view_infos = NULL;
char* view_name = NULL;
time_t creation_date;
char* creation_date_str = NULL;
index_t line_count;
char line_count_str[256];
// View name
view_name = (view->infos)->name;
view_infos = (char*) malloc((strlen(" # ")+strlen(view_name)+2) * sizeof(char));
strcpy(view_infos, " # ");
strcat(view_infos, view_name);
strcat(view_infos, ":");
// Date created
if (view->read_only) // Date not saved until view is finished writing
{
creation_date = (view->infos)->creation_date;
creation_date_str = ctime(&creation_date);
// Delete \n added by ctime
creation_date_str[strlen(creation_date_str)-1] = '\0';
view_infos = realloc(view_infos, (strlen(view_infos)+strlen(" Date created: ")+strlen(creation_date_str)+1) * sizeof(char));
strcat(view_infos, " Date created: ");
strcat(view_infos, creation_date_str);
}
// Line count
line_count = (view->infos)->line_count;
snprintf(line_count_str, sizeof line_count_str, "%lld", line_count);
view_infos = realloc(view_infos, (strlen(view_infos)+strlen(" ; Line count: ")+strlen(line_count_str)+1) * sizeof(char));
strcat(view_infos, " ; Line count: ");
strcat(view_infos, line_count_str);
view_infos = realloc(view_infos, (strlen(view_infos)+2) * sizeof(char));
strcat(view_infos, "\n");
return view_infos;
}
int obi_view_write_comments(Obiview_p view, const char* comments)
{
size_t new_size;
@ -3047,7 +3274,7 @@ int obi_create_auto_id_column(Obiview_p view, const char* prefix)
// Delete old ID column if it exists
if (obi_view_get_column(view, ID_COLUMN) != NULL)
{
if (obi_view_delete_column(view, ID_COLUMN) < 0)
if (obi_view_delete_column(view, ID_COLUMN, false) < 0)
{
obidebug(1, "Error deleting an ID column to replace it in a view");
return -1;

View File

@ -30,45 +30,56 @@
#include "obiblob.h"
#define OBIVIEW_NAME_MAX_LENGTH (249) /**< The maximum length of an OBIDMS view name, without the extension.
*/
#define VIEW_TYPE_MAX_LENGTH (1024) /**< The maximum length of the type name of a view.
*/
#define LINES_COLUMN_NAME "LINES" /**< The name of the column containing the line selections
* in all views.
*/
#define VIEW_TYPE_NUC_SEQS "NUC_SEQS_VIEW" /**< The type name of views based on nucleotide sequences
* and their metadata.
*/
#define NUC_SEQUENCE_COLUMN "NUC_SEQ" /**< The name of the column containing the nucleotide sequences
* in NUC_SEQS_VIEW views.
*/
#define ID_COLUMN "ID" /**< The name of the column containing the sequence identifiers
* in NUC_SEQS_VIEW views.
*/
#define DEFINITION_COLUMN "DEFINITION" /**< The name of the column containing the sequence definitions
* in NUC_SEQS_VIEW views.
*/
#define QUALITY_COLUMN "QUALITY" /**< The name of the column containing the sequence qualities
* in NUC_SEQS_VIEW views.
*/
#define COUNT_COLUMN "COUNT" /**< The name of the column containing the sequence counts
* in NUC_SEQS_VIEW views.
*/
#define TAXID_COLUMN "TAXID" /**< The name of the column containing the taxids. TODO subtype of INT column?
*/
#define MERGED_TAXID_COLUMN "MERGED_TAXID" /**< The name of the column containing the merged taxids information.
*/
#define MERGED_PREFIX "MERGED_" /**< The prefix to prepend to column names when merging informations during obi uniq.
*/
#define TAXID_DIST_COLUMN "TAXID_DIST" /**< The name of the column containing a dictionary of taxid:[list of ids] when merging informations during obi uniq.
*/
#define MERGED_COLUMN "MERGED" /**< The name of the column containing a list of ids when merging informations during obi uniq.
*/
#define ID_PREFIX "seq" /**< The default prefix of sequence identifiers in automatic ID columns.
*/
#define PREDICATE_KEY "predicates" /**< The key used in the json-formatted view comments to store predicates.
*/
#define OBIVIEW_NAME_MAX_LENGTH (249) /**< The maximum length of an OBIDMS view name, without the extension.
*/
#define VIEW_TYPE_MAX_LENGTH (1024) /**< The maximum length of the type name of a view.
*/
#define LINES_COLUMN_NAME "LINES" /**< The name of the column containing the line selections
* in all views.
*/
#define VIEW_TYPE_NUC_SEQS "NUC_SEQS_VIEW" /**< The type name of views based on nucleotide sequences
* and their metadata.
*/
#define NUC_SEQUENCE_COLUMN "NUC_SEQ" /**< The name of the column containing the nucleotide sequences
* in NUC_SEQS_VIEW views.
*/
#define ID_COLUMN "ID" /**< The name of the column containing the sequence identifiers
* in NUC_SEQS_VIEW views.
*/
#define DEFINITION_COLUMN "DEFINITION" /**< The name of the column containing the sequence definitions
* in NUC_SEQS_VIEW views.
*/
#define QUALITY_COLUMN "QUALITY" /**< The name of the column containing the sequence qualities
* in NUC_SEQS_VIEW views.
*/
#define REVERSE_QUALITY_COLUMN "REVERSE_QUALITY" /**< The name of the column containing the sequence qualities
* of the reverse read (generated by ngsfilter, used by alignpairedend).
*/
#define REVERSE_SEQUENCE_COLUMN "REVERSE_SEQUENCE" /**< The name of the column containing the sequence
* of the reverse read (generated by ngsfilter, used by alignpairedend).
*/
#define QUALITY_COLUMN "QUALITY" /**< The name of the column containing the sequence qualities
* in NUC_SEQS_VIEW views.
*/
#define COUNT_COLUMN "COUNT" /**< The name of the column containing the sequence counts
* in NUC_SEQS_VIEW views.
*/
#define SCIENTIFIC_NAME_COLUMN "SCIENTIFIC_NAME" /**< The name of the column containing the taxon scientific name.
*/
#define TAXID_COLUMN "TAXID" /**< The name of the column containing the taxids. TODO subtype of INT column?
*/
#define MERGED_TAXID_COLUMN "MERGED_TAXID" /**< The name of the column containing the merged taxids information.
*/
#define MERGED_PREFIX "MERGED_" /**< The prefix to prepend to column names when merging informations during obi uniq.
*/
#define TAXID_DIST_COLUMN "TAXID_DIST" /**< The name of the column containing a dictionary of taxid:[list of ids] when merging informations during obi uniq.
*/
#define MERGED_COLUMN "MERGED" /**< The name of the column containing a list of ids when merging informations during obi uniq.
*/
#define ID_PREFIX "seq" /**< The default prefix of sequence identifiers in automatic ID columns.
*/
#define PREDICATE_KEY "predicates" /**< The key used in the json-formatted view comments to store predicates.
*/
/**
@ -397,7 +408,7 @@ Obiview_p obi_open_view(OBIDMS_p dms, const char* view_name);
* @param associated_column_name The name of the associated column if there is one (otherwise NULL or ""), if the column is created.
* @param associated_column_version The version of the associated column if there is one (otherwise -1), if the column is created.
* @param comments Optional comments associated with the column if it is created (NULL or "" if no comments associated).
* @param create Whether the column should be created (create == true) or opened (create == false).
* @param create Whether the column should be created (create == true) or already exists (create == false).
*
* @returns A value indicating the success of the operation.
* @retval 0 if the operation was successfully completed.
@ -407,7 +418,7 @@ Obiview_p obi_open_view(OBIDMS_p dms, const char* view_name);
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
int obi_view_add_column(Obiview_p view,
char* column_name,
char* column_name,
obiversion_t version_number,
const char* alias,
OBIType_t data_type,
@ -431,6 +442,7 @@ int obi_view_add_column(Obiview_p view,
*
* @param view A pointer on the view.
* @param column_name The name of the column that should be deleted from the view.
* @param delete_file Whether the column file should be deleted. Use carefully re: dependencies.
*
* @returns A value indicating the success of the operation.
* @retval 0 if the operation was successfully completed.
@ -439,7 +451,7 @@ int obi_view_add_column(Obiview_p view,
* @since February 2016
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
int obi_view_delete_column(Obiview_p view, const char* column_name);
int obi_view_delete_column(Obiview_p view, const char* column_name, bool delete_file);
/**
@ -509,6 +521,39 @@ OBIDMS_column_p* obi_view_get_pointer_on_column_in_view(Obiview_p view, const ch
int obi_view_create_column_alias(Obiview_p view, const char* current_name, const char* alias);
/**
* @brief Returns the informations of a view with a human readable format (view name, date created, line count, column informations, comments).
*
* @warning The returned pointer has to be freed by the caller.
*
* @param column A pointer on a view.
* @param detailed Whether the informations should contain view comments.
*
* @returns A pointer on a character array where the formatted view informations are stored.
* @retval NULL if an error occurred.
*
* @since September 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_view_formatted_infos(Obiview_p view, bool detailed);
/**
* @brief Returns the informations of a view with a human readable format on one line (view name, date created, line count).
*
* @warning The returned pointer has to be freed by the caller.
*
* @param column A pointer on a view.
*
* @returns A pointer on a character array where the formatted view informations are stored.
* @retval NULL if an error occurred.
*
* @since September 2020
* @author Celine Mercier (celine.mercier@metabarcoding.org)
*/
char* obi_view_formatted_infos_one_line(Obiview_p view);
/**
* @brief Internal function writing new comments in a view file.
*

View File

@ -686,6 +686,9 @@ int calculateSizeToAllocate(int maxLen, int LCSmin)
size *= 3;
size += 16;
size += 10; // band-aid for memory bug I don't understand (triggered on specific db on ubuntu)
// bug might have to do with the way different systems behave when aligning the address in obi_get_memory_aligned_on_16
return(size*sizeof(int16_t));
}
@ -951,15 +954,15 @@ double generic_sse_banded_lcs_align(char* seq1, char* seq2, double threshold, bo
// Put the DNA sequences in the int arrays. Longest sequence must be first argument of sse_align function
if (l2 > l1)
{
putSeqInSeq(iseq1, seq2, l2, TRUE);
putSeqInSeq(iseq2, seq1, l1, FALSE);
putSeqInSeq(iseq1, seq2, l2, true);
putSeqInSeq(iseq2, seq1, l1, false);
// Compute alignment
id = sse_banded_lcs_align(iseq1, iseq2, l2, l1, normalize, reference, similarity_mode, address, LCSmin, lcs_length, ali_length);
}
else
{
putSeqInSeq(iseq1, seq1, l1, TRUE);
putSeqInSeq(iseq2, seq2, l2, FALSE);
putSeqInSeq(iseq1, seq1, l1, true);
putSeqInSeq(iseq2, seq2, l2, false);
// Compute alignment
id = sse_banded_lcs_align(iseq1, iseq2, l1, l2, normalize, reference, similarity_mode, address, LCSmin, lcs_length, ali_length);
}
@ -1054,15 +1057,15 @@ double obiblob_sse_banded_lcs_align(Obi_blob_p seq1, Obi_blob_p seq2, double thr
// Put the DNA sequences in the int arrays. Longest sequence must be first argument of sse_align function
if (l2 > l1)
{
putBlobInSeq(iseq1, seq2, l2, TRUE);
putBlobInSeq(iseq2, seq1, l1, FALSE);
putBlobInSeq(iseq1, seq2, l2, true);
putBlobInSeq(iseq2, seq1, l1, false);
// Compute alignment
id = sse_banded_lcs_align(iseq1, iseq2, l2, l1, normalize, reference, similarity_mode, address, LCSmin, lcs_length, ali_length);
}
else
{
putBlobInSeq(iseq1, seq1, l1, TRUE);
putBlobInSeq(iseq2, seq2, l2, FALSE);
putBlobInSeq(iseq1, seq1, l1, true);
putBlobInSeq(iseq2, seq2, l2, false);
// Compute alignment
id = sse_banded_lcs_align(iseq1, iseq2, l1, l2, normalize, reference, similarity_mode, address, LCSmin, lcs_length, ali_length);
}