Source code for tool.tb_parse_mapping
"""
.. See the NOTICE file distributed with this work for additional information
regarding copyright ownership.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
from __future__ import print_function
import sys
from basic_modules.tool import Tool
from utils import logger
try:
if hasattr(sys, '_run_from_cmdl') is True:
raise ImportError
from pycompss.api.parameter import FILE_IN, FILE_OUT, IN
from pycompss.api.task import task
# from pycompss.api.constraint import constraint
from pycompss.api.api import compss_wait_on
except ImportError:
logger.info("[Warning] Cannot import \"pycompss\" API packages.")
logger.info(" Using mock decorators.")
from utils.dummy_pycompss import FILE_IN, FILE_OUT, IN # pylint: disable=ungrouped-imports
from utils.dummy_pycompss import task # pylint: disable=ungrouped-imports
# from utils.dummy_pycompss import constraint # pylint: disable=ungrouped-imports
from utils.dummy_pycompss import compss_wait_on # pylint: disable=ungrouped-imports
from pytadbit.parsers.genome_parser import parse_fasta # pylint: disable=import-error
from pytadbit.parsers.map_parser import parse_map # pylint: disable=import-error
from pytadbit.mapping import get_intersection # pylint: disable=import-error
# ------------------------------------------------------------------------------
[docs]class tbParseMappingTool(Tool): # pylint: disable=invalid-name
"""
Tool for parsing the mapped reads and generating the list of paired ends
that have a match at both ends.
"""
def __init__(self):
"""
Init function
"""
logger.info("TADbit parse mapping")
Tool.__init__(self)
[docs] @task(
genome_seq=IN, enzyme_name=IN,
window1_1=FILE_IN, window1_2=FILE_IN, window1_3=FILE_IN, window1_4=FILE_IN,
window2_1=FILE_IN, window2_2=FILE_IN, window2_3=FILE_IN, window2_4=FILE_IN,
reads=FILE_OUT)
def tb_parse_mapping_iter( # pylint: disable=no-self-use,too-many-arguments,too-many-locals
self, genome_seq, enzyme_name,
window1_1, window1_2, window1_3, window1_4,
window2_1, window2_2, window2_3, window2_4,
reads, ncpus=1):
"""
Function to map the aligned reads and return the matching pairs
Parameters
----------
genome_seq : dict
Object containing the sequence of each of the chromosomes
enzyme_name : str
Name of the enzyme used to digest the genome
window1_1 : str
Location of the first window index file
window1_2 : str
Location of the second window index file
window1_3 : str
Location of the third window index file
window1_4 : str
Location of the fourth window index file
window2_1 : str
Location of the first window index file
window2_2 : str
Location of the second window index file
window2_3 : str
Location of the third window index file
window2_4 : str
Location of the fourth window index file
reads : str
Location of the reads thats that has a matching location at both
ends of the paired reads
Returns
-------
reads : str
Location of the intersection of mapped reads that have matching
reads in both pair end files
"""
reads1 = reads + '_reads_1.tsv'
reads2 = reads + '_reads_2.tsv'
reads_both = reads + '_reads_both.tsv'
wind1 = [window1_1]
wind2 = [window2_1]
if window1_2 and window2_2:
wind1 += [window1_2]
wind2 += [window2_2]
if window1_3 and window2_3:
wind1 += [window1_3]
wind2 += [window2_3]
if window1_4 and window2_4:
wind1 += [window1_4]
wind2 += [window2_4]
parse_map(
wind1,
wind2,
out_file1=reads1,
out_file2=reads2,
genome_seq=genome_seq,
re_name=enzyme_name,
verbose=True,
ncpus=ncpus
)
counts, _ = get_intersection(reads1, reads2, reads_both, verbose=True)
with open(reads, "wb") as f_out:
with open(reads_both, "rb") as f_in:
f_out.write(f_in.read())
return counts
[docs] @task(
genome_seq=IN, enzyme_name=IN,
window1_full=FILE_IN, window1_frag=FILE_IN,
window2_full=FILE_IN, window2_frag=FILE_IN,
reads=FILE_OUT)
def tb_parse_mapping_frag( # pylint: disable=no-self-use,too-many-arguments,too-many-locals
self, genome_seq, enzyme_name,
window1_full, window1_frag, window2_full, window2_frag,
reads, ncpus=1):
"""
Function to map the aligned reads and return the matching pairs
Parameters
----------
genome_seq : dict
Object containing the sequence of each of the chromosomes
enzyme_name : str
Name of the enzyme used to digest the genome
window1_full : str
Location of the first window index file
window1_frag : str
Location of the second window index file
window2_full : str
Location of the first window index file
window2_frag : str
Location of the second window index file
reads : str
Location of the reads thats that has a matching location at both
ends of the paired reads
Returns
-------
reads : str
Location of the intersection of mapped reads that have matching
reads in both pair end files
"""
logger.info("TB WINDOWS - full 1 {0}".format(window1_full))
logger.info("TB WINDOWS - frag 1 {0}".format(window1_frag))
logger.info("TB WINDOWS - full 2 {0}".format(window2_full))
logger.info("TB WINDOWS - frag 2 {0}".format(window2_frag))
# root_name = reads.split("/")
# reads1 = "/".join(root_name) + '/reads_1.tsv'
# reads2 = "/".join(root_name) + '/reads_2.tsv'
reads1 = reads + '_reads_1.tsv'
reads2 = reads + '_reads_2.tsv'
reads_both = reads + '_reads_both.tsv'
parse_map(
[window1_frag, window1_full],
[window2_frag, window2_full],
out_file1=reads1,
out_file2=reads2,
genome_seq=genome_seq,
re_name=enzyme_name,
verbose=True,
ncpus=ncpus
)
counts, _ = get_intersection(reads1, reads2, reads_both, verbose=True)
with open(reads, "wb") as f_out:
with open(reads_both, "rb") as f_in:
f_out.write(f_in.read())
return counts
[docs] def run(self, input_files, input_metadata, output_files): # pylint: disable=too-many-locals
"""
The main function to map the aligned reads and return the matching
pairs. Parsing of the mappings can be either iterative of fragment
based. If it is to be iteractive then the locations of 4 output file
windows for each end of the paired end window need to be provided. If
it is fragment based, then only 2 window locations need to be provided
along within an enzyme name.
Parameters
----------
input_files : list
genome_file : str
Location of the genome FASTA file
window1_1 : str
Location of the first window index file
window1_2 : str
Location of the second window index file
window1_3 : str
[OPTIONAL] Location of the third window index file
window1_4 : str
[OPTIONAL] Location of the fourth window index file
window2_1 : str
Location of the first window index file
window2_2 : str
Location of the second window index file
window2_3 : str
[OPTIONAL] Location of the third window index file
window2_4 : str
[OPTIONAL] Location of the fourth window index file
metadata : dict
windows : list
List of lists with the window sizes to be computed
enzyme_name : str
Restricture enzyme name
mapping : list
The mapping function used. The options are iter or frag.
Returns
-------
output_files : list
List of locations for the output files.
output_metadata : dict
Dict of matching metadata dict objects
Example
-------
Iterative:
.. code-block:: python
from tool import tb_parse_mapping
genome_file = 'genome.fasta'
root_name_1 = "/tmp/data/expt_source_1".split
root_name_2 = "/tmp/data/expt_source_2".split
windows = [[1,25], [1,50], [1,75], [1,100]]
windows1 = []
windows2 = []
for w in windows:
tail = "_full_" + w[0] + "-" + w[1] + ".map"
windows1.append('/'.join(root_name_1) + tail)
windows2.append('/'.join(root_name_2) + tail)
files = [genome_file] + windows1 + windows2
tpm = tb_parse_mapping.tb_parse_mapping()
metadata = {'enzyme_name' : 'MboI', 'mapping' : ['iter', 'iter'], 'expt_name' = 'test'}
tpm_files, tpm_meta = tpm.run(files, metadata)
Fragment based mapping:
.. code-block:: python
from tool import tb_parse_mapping
genome_file = 'genome.fasta'
root_name_1 = "/tmp/data/expt_source_1".split
root_name_2 = "/tmp/data/expt_source_2".split
windows = [[1,100]]
start = windows[0][0]
end = windows[0][1]
window1_1 = '/'.join(root_name_1) + "_full_" + start + "-" + end + ".map"
window1_2 = '/'.join(root_name_1) + "_frag_" + start + "-" + end + ".map"
window2_1 = '/'.join(root_name_2) + "_full_" + start + "-" + end + ".map"
window2_2 = '/'.join(root_name_2) + "_frag_" + start + "-" + end + ".map"
files = [
genome_file,
window1_1, window1_2,
window2_1, window2_2,
]
tpm = tb_parse_mapping.tb_parse_mapping()
metadata = {'enzyme_name' : 'MboI', 'mapping' : ['frag', 'frag'], 'expt_name' = 'test'}
tpm_files, tpm_meta = tpm.run(files, metadata)
"""
genome_file = input_files[0]
enzyme_name = input_metadata['enzyme_name']
mapping_list = input_metadata['mapping']
expt_name = input_metadata['expt_name']
filter_chrom = None
if 'chromosomes' in input_metadata and input_metadata['chromosomes'] != '':
filter_chrom = input_metadata['chromosomes'].split(',')
root_name = input_files[1].split("/")
reads = "/".join(root_name[0:-1]) + '/'
genome_seq = parse_fasta( # pylint: disable=unexpected-keyword-arg
genome_file, chr_filter=filter_chrom,
chr_regexp="^(chr)?[A-Za-z]?[0-9]{0,3}[XVI]{0,3}(?:ito)?[A-Z-a-z]?$",
save_cache=False, reload_cache=True)
if not genome_seq:
genome_seq = parse_fasta( # pylint: disable=unexpected-keyword-arg
genome_file, chr_filter=filter_chrom,
save_cache=False, reload_cache=True)
if not genome_seq:
logger.fatal("Reference genome FASTA file does not contain any valid chromosome")
raise ValueError('No valid chromosomes in FASTA.')
else:
logger.warn("No chromosomes headers found in fasta, found {0}.".format(
[k for k in genome_seq]))
logger.warn("Using them instead")
chromosome_meta = []
for k in genome_seq:
chromosome_meta.append([k, len(genome_seq[k])])
# input and output share most metadata
output_metadata = {
'chromosomes': chromosome_meta
}
if mapping_list[0] == mapping_list[1]:
if mapping_list[0] == 'iter':
window1_1 = input_files[1]
window1_2 = input_files[2]
window1_3 = input_files[3]
window1_4 = input_files[4]
window2_1 = input_files[5]
window2_2 = input_files[6]
window2_3 = input_files[7]
window2_4 = input_files[8]
read_iter = reads + expt_name + '_iter.tsv'
results = self.tb_parse_mapping_iter(
genome_seq, enzyme_name,
window1_1, window1_2, window1_3, window1_4,
window2_1, window2_2, window2_3, window2_4,
read_iter, ncpus=(
1 if 'ncpus' not in input_metadata else input_metadata['ncpus']))
results = compss_wait_on(results)
if results == 0:
output_metadata = {
'error': 'No interactions found, \
please verify input data and chromosome filtering'
}
return ([], output_metadata)
return ([read_iter], output_metadata)
if mapping_list[0] == 'frag':
window1_full = input_files[1]
window1_frag = input_files[2]
window2_full = input_files[3]
window2_frag = input_files[4]
read_frag = reads + expt_name + '_frag.tsv'
results = self.tb_parse_mapping_frag(
genome_seq, enzyme_name,
window1_full, window1_frag,
window2_full, window2_frag,
read_frag, ncpus=(
1 if 'ncpus' not in input_metadata else input_metadata['ncpus']))
results = compss_wait_on(results)
if results == 0:
output_metadata = {
'error': 'No interactions found, \
please verify input data and chromosome filtering'
}
return ([], output_metadata)
return ([read_frag], output_metadata)
reads = None
return ([reads], output_metadata)
return ([], [])
# ------------------------------------------------------------------------------