"""
09 Jun 2015
"""
from __future__ import print_function
from builtins import next
import os
import re
import copy
import multiprocessing as mu
from warnings import warn
from tempfile import gettempdir, mkstemp
from subprocess import CalledProcessError, PIPE, STDOUT, Popen
from pysam import Samfile
from pytadbit.utils.file_handling import mkdir, which, is_fastq
from pytadbit.utils.file_handling import magic_open, get_free_space_mb
from pytadbit.parsers.sam_parser import parse_gem_3c, merge_sort
from pytadbit.mapping.restriction_enzymes import religateds
from pytadbit.mapping.restriction_enzymes import RESTRICTION_ENZYMES
from pytadbit.mapping.restriction_enzymes import map_re_sites
from pytadbit.mapping.restriction_enzymes import iupac2regex
try:
basestring
except NameError:
basestring = str
def transform_fastq(fastq_path, out_fastq, trim=None, r_enz=None, add_site=True,
min_seq_len=15, fastq=True, verbose=True,
light_storage=False, **kwargs):
"""
Given a FASTQ file it can split it into chunks of a given number of reads,
trim each read according to a start/end positions or split them into
restriction enzyme fragments
:param True add_site: when splitting the sequence by ligated sites found,
removes the ligation site, and put back the original RE site.
"""
skip = kwargs.get('skip', False)
## define local functions to process reads and sequences
def _get_fastq_read_heavy(rlines):
"""
returns header and sequence of 1 FASTQ entry
Note: header also contains the sequence
"""
rlines = rlines.rstrip('\n').split()[0][1:]
seq = next(fhandler)
_ = next(fhandler) # lose qualities header but not needed
qal = next(fhandler) # lose qualities but not needed
# header now also contains original read
return (rlines.split('/',1)[0].split('~',1)[0] + ' ' + seq.strip() + ' ' + qal.strip(),
seq.strip(), qal.strip())
def _get_fastq_read_light(rlines):
"""
returns header and sequence of 1 FASTQ entry
Note: header also contains the sequence
"""
rlines = rlines.rstrip('\n').split()[0][1:]
seq = next(fhandler)
_ = next(fhandler) # lose qualities header but not needed
qal = next(fhandler) # lose qualities but not needed
return (rlines.split('/',1)[0].split('~',1)[0], seq.strip(), qal.strip())
def _get_map_read_heavy(line):
header = line.split('\t', 1)[0]
seq, qal = header.rsplit(' ', 2)[-2:]
return header, seq, qal
def _get_map_read_light(line):
header, seq, qal, _ = line.split('\t', 3)
return header, seq, qal
def _inverse_find(seq, pat):
try:
return pat.search(seq).start()
except AttributeError:
return float('inf')
def find_patterns(seq, patterns):
pos, pat = min((_inverse_find(seq, patterns[p]), p) for p in patterns)
return int(pos), pat
def _split_read_re(seq, qal, patterns, site, max_seq_len=None, cnt=0):
"""
Recursive generator that splits reads according to the
predefined restriction enzyme.
RE fragments yielded are followed and preceded by the RE site if a
ligation site was found after the fragment.
EXAMPLE:
seq = '-------oGATCo========oGATCGATCo_____________oGATCGATCo~~~~~~~~~~~~'
qal = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
should yield these fragments:
-------oGATCo========oGATC
xxxxxxxxxxxxxxxxxxxxxxHHHH
GATCo_____________oGATC
HHHHxxxxxxxxxxxxxxxHHHH
GATCo~~~~~~~~~~~~
HHHHxxxxxxxxxxxxx
:param seq: sequence of the read fragment
:param qal: quality of the sequence of the read fragment
:param patterns: list of patterns of the ligated cut sites
:param None max_seq_len: to control that all reads are bellow this
length
:param '' site: non-ligated cut site to replace ligation site
:param 0 cnt: to count number of fragments
:yields: seq fragments, their qualities and their count, or index
(higher than 0 if ligation sites are found)
"""
cnt += 1
try:
pos, (r_enz1, r_enz2) = find_patterns(seq, patterns)
except OverflowError:
if len(seq) == max_seq_len:
raise ValueError
if len(seq) > min_seq_len:
yield seq, qal, cnt
return
# add quality before corresponding to the space occupied by the cut-site
xqal1 = ('H' * len(site[r_enz1]))
xqal2 = ('H' * len(site[r_enz2]))
if pos < min_seq_len:
split_read(site[r_enz2] + seq[pos + len_relgs[(r_enz1, r_enz2)]:],
xqal2 + qal[pos + len_relgs[(r_enz1, r_enz2)]:],
patterns, no_site, max_seq_len, cnt=cnt)
else:
yield seq[:pos] + site[r_enz1], qal[:pos] + xqal1, cnt
new_pos = pos + len_relgs[(r_enz1, r_enz2)]
for sseq, sqal, cnt in split_read(site[r_enz2] + seq[new_pos:],
xqal2 + qal[new_pos:], patterns,
site, max_seq_len, cnt=cnt):
yield sseq, sqal, cnt
# Define function for stripping lines according to focus
if isinstance(trim, tuple):
beg, end = trim
beg -= 1
strip_line = lambda x: x[beg:end]
else:
strip_line = lambda x: x
# define function to split reads according to restriction enzyme sites
if isinstance(r_enz, basestring):
r_enzs = [r_enz]
elif isinstance(r_enz, list):
r_enzs = r_enz
else:
r_enzs = None
if r_enzs:
enzymes = {}
enz_patterns = {}
for r_enz in r_enzs:
enzymes[r_enz] = RESTRICTION_ENZYMES[r_enz].replace('|', '')
enz_patterns = religateds(r_enzs)
sub_enz_patterns = {}
len_relgs = {}
for r_enz1, r_enz2 in enz_patterns:
sub_enz_patterns[(r_enz1, r_enz2)] = (
enz_patterns[(r_enz1, r_enz2)][:len(enz_patterns[(r_enz1, r_enz2)])
// 2])
len_relgs[(r_enz1, r_enz2)] = len(enz_patterns[(r_enz1, r_enz2)])
print(' - splitting into restriction enzyme (RE) fragments using ligation sites')
print(' - ligation sites are replaced by RE sites to match the reference genome')
for r_enz1 in r_enzs:
for r_enz2 in r_enzs:
print(' * enzymes: %s & %s, ligation site: %s, RE site: %s & %s' % (
r_enz1, r_enz2, enz_patterns[(r_enz1, r_enz2)],
enzymes[r_enz1], enzymes[r_enz2]))
# replace pattern with regex to support IUPAC annotation
for ezp in enz_patterns:
enz_patterns[ezp] = re.compile(iupac2regex(enz_patterns[ezp]))
for ezp in sub_enz_patterns:
sub_enz_patterns[ezp] = iupac2regex(sub_enz_patterns[ezp])
split_read = _split_read_re
else:
enzymes = ''
enz_patterns = ''
sub_enz_patterns = ''
split_read = lambda x, y, z, after_z, after_after_z: (yield x, y , 1)
# function to yield reads from input file
if light_storage:
get_seq = _get_fastq_read_light if fastq else _get_map_read_light
insert_mark = insert_mark_light
else:
get_seq = _get_fastq_read_heavy if fastq else _get_map_read_heavy
insert_mark = insert_mark_heavy
## Start processing the input file
if verbose:
print('Preparing %s file' % ('FASTQ' if fastq else 'MAP'))
if fastq:
print(' - conversion to MAP format')
if trim:
print(' - trimming reads %d-%d' % tuple(trim))
counter = 0
if skip:
if fastq:
print(' ... skipping, only counting lines')
counter = sum(1 for _ in magic_open(fastq_path,
cpus=kwargs.get('nthreads')))
counter /= 4 if fastq else 1
print(' ' + fastq_path, counter, fastq)
return out_fastq, counter
# open input file
fhandler = magic_open(fastq_path, cpus=kwargs.get('nthreads'))
# create output file
out_name = out_fastq
out = open(out_fastq, 'w')
# iterate over reads and strip them
no_site = dict([(r_enz, '') for r_enz in enzymes])
site = enzymes if add_site else no_site
for header in fhandler:
header, seq, qal = get_seq(header)
counter += 1
# trim on wanted region of the read
seq = strip_line(seq)
qal = strip_line(qal)
# get the generator of restriction enzyme fragments
iter_frags = split_read(seq, qal, enz_patterns, site, len(seq))
# the first fragment should not be preceded by the RE site
try:
seq, qal, cnt = next(iter_frags)
except StopIteration:
# read full of ligation events, fragments not reaching minimum
continue
except ValueError:
# or not ligation site found, in which case we try with half
# ligation site in case there was a sequencing error (half ligation
# site is a RE site or nearly, and thus should not be found anyway)
iter_frags = split_read(seq, qal, sub_enz_patterns, no_site, len(seq))
try:
seq, qal, cnt = next(iter_frags)
except ValueError:
continue
except StopIteration:
continue
out.write(_map2fastq('\t'.join((insert_mark(header, cnt),
seq, qal, '0', '-\n'))))
# the next fragments should be preceded by the RE site
# continue
for seq, qal, cnt in iter_frags:
out.write(_map2fastq('\t'.join((insert_mark(header, cnt),
seq, qal, '0', '-\n'))))
out.close()
return out_name, counter
def insert_mark_heavy(header, num):
if num == 1 :
return header
h, s, q = header.rsplit(' ', 2)
return '%s~%d~ %s %s' % (h, num, s, q)
def insert_mark_light(header, num):
if num == 1 :
return header
return '%s~%d~' % (header, num)
def _map2fastq(read):
return '@{0}\n{1}\n+\n{2}\n'.format(*read.split('\t', 3)[:-1])
def _sam_filter(fnam, fastq_path, unmap_out, map_out):
"""
Divides reads in a map file in two categories: uniquely mapped, and not.
Writes them in two files
"""
try:
fhandler = Samfile(fnam)
except IOError:
raise Exception('ERROR: file "%s" not found' % fnam)
# getrname chromosome names
i = 0
crm_dict = {}
while True:
try:
crm_dict[i] = fhandler.getrname(i)
i += 1
except ValueError:
break
# iteration over reads
unmap_out = open(unmap_out, 'w')
map_out = open(map_out, 'w')
fastq_in = open(fastq_path , 'r')
for line in fhandler:
line_in = fastq_in.readline()
if line.is_unmapped or line.mapq < 4:
read = '%s\t%s\t%s\t%s\t%s\n' % (
line_in.split('\t', 1)[0].rstrip('\n')[1:],
line.seq, line.qual, '-', '-'
)
unmap_out.write(read)
else:
read = '%s\t%s\t%s\t%s\t%s:%s:%d:%d\n' % (
line.qname, line.seq, line.qual, '1',
crm_dict[line.tid],
'-' if line.is_reverse else '+', line.pos + 1, len(line.seq))
map_out.write(read)
for _ in range(3):
fastq_in.readline()
unmap_out.close()
map_out.close()
fastq_in.close()
def _gem_filter(fnam, unmap_out, map_out):
"""
Divides reads in a map file in two categories: uniquely mapped, and not.
Writes them in two files
Notes:
- GEM unique-maps can not be used as it gets rid of reads like 1:0:0:5
- not feasible with gt.filter
"""
fhandler = magic_open(fnam) if isinstance(fnam, basestring) else fnam
unmap_out = open(unmap_out, 'w')
map_out = open(map_out , 'w')
def _strip_read_name(line):
"""
remove original sequence from read name when read is mapped uniquely
"""
header, line = line.split('\t', 1)
return '\t'.join((header.rsplit(' ', 2)[0], line))
for line in fhandler:
matches = line.rsplit('\t', 2)[1]
bad = False
if matches != '1':
for m in matches.replace('+', ':').split(':'):
if m == '0':
continue
if m != '1':
bad = True
unmap_out.write(line)
break
break
else:
bad = True
unmap_out.write(line)
if not bad:
map_out.write(_strip_read_name(line))
unmap_out.close()
map_out.close()
def _bowtie2_mapping(bowtie2_index_path, fastq_path1, out_map_path, fastq_path2 = None,
bowtie2_binary='bowtie2', bowtie2_params=None, **kwargs):
"""
"""
bowtie2_index_path= os.path.abspath(os.path.expanduser(bowtie2_index_path))
fastq_path1 = os.path.abspath(os.path.expanduser(fastq_path1))
paired_map = False
if fastq_path2:
fastq_path2 = os.path.abspath(os.path.expanduser(fastq_path2))
paired_map = True
out_map_path = os.path.abspath(os.path.expanduser(out_map_path))
nthreads = kwargs.get('nthreads' , 8)
# check that we have the GEM binary:
bowtie2_binary = which(bowtie2_binary)
if not bowtie2_binary:
raise Exception('\n\nERROR: %s binary not found'%bowtie2_binary)
# mapping
print('TO %s'%bowtie2_binary, fastq_path1, fastq_path2)
bowtie2_cmd = [
bowtie2_binary, '-x', bowtie2_index_path,
'-p', str(nthreads), '--reorder','-k','1','-S',
out_map_path]
if paired_map:
bowtie2_cmd += ['-1',fastq_path1,'-2',fastq_path2]
else:
bowtie2_cmd += ['-U', fastq_path1]
if bowtie2_params:
if isinstance(bowtie2_params, dict):
for bow_param in bowtie2_params:
bowtie2_cmd.append('-'+bow_param)
if bowtie2_params[bow_param]:
bowtie2_cmd.append(bowtie2_params[bow_param])
elif isinstance(bowtie2_params, list):
bowtie2_cmd += bowtie2_params
elif bowtie2_binary == 'bowtie2':
bowtie2_cmd.append('--very-sensitive')
print(' '.join(bowtie2_cmd))
try:
# check_call(gem_cmd, stdout=PIPE, stderr=PIPE)
out, err = Popen(bowtie2_cmd, stdout=PIPE, stderr=PIPE,
universal_newlines=True).communicate()
except CalledProcessError as e:
print(out)
print(err)
raise Exception(e.output)
def _gem_mapping(gem_index_path, fastq_path, out_map_path, fastq_path2 = None,
r_enz=None, gem_binary='gem-mapper', gem_version=2, compress=False,
gem_params=None, **kwargs):
"""
:param None focus: trims the sequence in the input FASTQ file according to a
(start, end) position, or the name of a restriction enzyme. By default it
uses the full sequence.
:param 33 quality: set it to 'ignore' in order to speed-up the mapping
"""
gem_index_path = os.path.abspath(os.path.expanduser(gem_index_path))
fastq_path = os.path.abspath(os.path.expanduser(fastq_path))
out_map_path = os.path.abspath(os.path.expanduser(out_map_path))
nthreads = kwargs.get('nthreads' , 8)
max_edit_distance = kwargs.get('max_edit_distance' , 0.04)
mismatches = kwargs.get('mismatches' , 0.04)
# mapping
print('TO GEM', gem_version, fastq_path)
kgt = kwargs.get
if gem_version == 2:
gem_cmd = [
gem_binary, '-I', gem_index_path,
'-q' , kgt('q', 'offset-33' ),
'-m' , kgt('m', str(max_edit_distance ) ),
'-s' , kgt('s', kgt('strata-after-best', '0') ),
'--allow-incomplete-strata' , kgt('allow-incomplete-strata', '0.00' ),
'--granularity' , kgt('granularity', '10000' ),
'--max-decoded-matches' , kgt('max-decoded-matches', kgt('d', '1')),
'--min-decoded-strata' , kgt('min-decoded-strata', kgt('D', '0') ),
'--min-insert-size' , kgt('min-insert-size', '0' ),
'--max-insert-size' , kgt('max-insert-size', '0' ),
'--min-matched-bases' , kgt('min-matched-bases', '0.8' ),
'--gem-quality-threshold' , kgt('gem-quality-threshold', '26' ),
'--max-big-indel-length' , kgt('max-big-indel-length', '15' ),
'--mismatch-alphabet' , kgt('mismatch-alphabet', 'ACGT' ),
'-E' , kgt('E', '0.30' ),
'--max-extendable-matches' , kgt('max-extendable-matches', '20' ),
'--max-extensions-per-match', kgt('max-extensions-per-match', '1' ),
'-e' , kgt('e', str(mismatches) ),
'-T' , str(nthreads),
'-i' , fastq_path,
'-o', out_map_path.replace('.map', '')]
if 'paired-end-alignment' in kwargs or 'p' in kwargs:
gem_cmd.append('--paired-end-alignment')
if 'map-both-ends' in kwargs or 'b' in kwargs:
gem_cmd.append('--map-both-ends')
if 'fast-mapping' in kwargs:
gem_cmd.append('--fast-mapping')
if 'unique-mapping' in kwargs:
gem_cmd.append('--unique-mapping')
if 'unique-pairing' in kwargs:
gem_cmd.append('--unique-pairing')
# check kwargs
for kw in kwargs:
if not kw in ['nthreads', 'max_edit_distance',
'mismatches', 'max_reads_per_chunk',
'out_files', 'temp_dir', 'skip', 'q', 'm', 's',
'strata-after-best', 'allow-incomplete-strata',
'granularity', 'max-decoded-matches',
'min-decoded-strata', 'min-insert-size',
'max-insert-size', 'min-matched-bases',
'gem-quality-threshold', 'max-big-indel-length',
'mismatch-alphabet', 'E', 'max-extendable-matches',
'max-extensions-per-match', 'e', 'paired-end-alignment',
'p', 'map-both-ends', 'fast-mapping', 'unique-mapping',
'unique-pairing', 'suffix']:
warn('WARNING: %s not in usual keywords, misspelled?' % kw)
else:
gem_cmd = [
gem_binary, '-I', gem_index_path,
'-t' , str(nthreads),
'-F' , 'SAM',
'-o', out_map_path]
if compress:
gem_cmd += ['--gzip-output']
if fastq_path2:
if not r_enz:
raise Exception('ERROR: need enzyme name to fragment.')
print('Using GEM', gem_version, 'with 3c mapping')
gem_cmd += ['--i1', fastq_path, '--i2', fastq_path2, '--3c']
if isinstance(r_enz, basestring):
gem_cmd += ['--restriction-enzyme', r_enz]
elif isinstance(r_enz, list):
for r_z in r_enz:
gem_cmd += ['--restriction-enzyme', r_z]
else:
gem_cmd += ['-i', fastq_path]
if gem_params:
if isinstance(gem_params, dict):
for gem_param in gem_params:
gem_cmd.append('-'+gem_param)
if gem_params[gem_param]:
gem_cmd.append(gem_params[gem_param])
elif isinstance(gem_params, list):
gem_cmd += gem_params
print(' '.join(gem_cmd))
try:
# check_call(gem_cmd, stdout=PIPE, stderr=PIPE)
out, err = Popen(gem_cmd, stdout=PIPE, stderr=PIPE,
universal_newlines=True).communicate()
except CalledProcessError as e:
print(out)
print(err)
raise Exception(e.output)
[docs]def full_mapping(mapper_index_path, fastq_path, out_map_dir, mapper='gem',
r_enz=None, frag_map=True, min_seq_len=15, windows=None,
add_site=True, clean=False, get_nread=False,
mapper_binary=None, mapper_params=None, **kwargs):
"""
Maps FASTQ reads to an indexed reference genome. Mapping can be done either
without knowledge of the restriction enzyme used, or for experiments
performed without one, like Micro-C (iterative mapping), or using the
ligation sites created from the digested ends (fragment-based mapping).
:param mapper_index_path: path to index file created from a reference genome
using gem-index tool or bowtie2-build
:param fastq_path: PATH to FASTQ file, either compressed or not.
:param out_map_dir: path to a directory where to store mapped reads in MAP
format .
:param None r_enz: name of the restriction enzyme used in the experiment e.g.
HindIII. This is optional if frag_map option is False
:param True frag_map: two step mapper, first full length is mapped, then
remaining, unmapped reads, are divided into restriction-enzyme fragments
andeach is mapped.
:param True add_site: when splitting the sequence by ligated sites found,
removes the ligation site, and put back the original RE site.
:param 15 min_seq_len: minimum size of a fragment to map
:param None windows: tuple of ranges for beginning and end of the
mapping. This parameter allows to do classical iterative mapping, e.g.
windows=((1,25),(1,30),(1,35),(1,40),(1,45),(1,50))
A unique window can also be passed, for trimming, like this:
windows=((1,101),)
:param False clean: remove intermediate files created in temp_dir
:param 4 nthreads: number of threads to use for mapping (number of CPUs)
:param 0.04 max_edit_distance: The maximum number of edit operations allowed
while verifying candidate matches by dynamic programming.
:param 0.04 mismatches: The maximum number of nucleotide substitutions
allowed while mapping each k-mer. It is always guaranteed that, however
other options are chosen, all the matches up to the specified number of
substitutions will be found by the program.
:param /tmp temp_dir: important to change. Intermediate FASTQ files will be
written there.
:param False get_nreads: returns a list of lists where each element contains
a path and the number of reads processed
:param gem-mapper mapper_binary: path to the binary mapper
:param None mapper_params: extra parameters for the mapper
:returns: a list of paths to generated outfiles. To be passed to
:func:`pytadbit.parsers.map_parser.parse_map`
"""
skip = kwargs.get('skip', False)
suffix = kwargs.get('suffix', '')
suffix = ('_' * (suffix != '')) + suffix
nthreads = kwargs.get('nthreads', 8)
outfiles = []
temp_dir = os.path.abspath(os.path.expanduser(
kwargs.get('temp_dir', gettempdir())))
if mapper == 'gem':
gem_version = None
# check that we have the GEM binary:
gem_binary = mapper_binary or 'gem-mapper'
gem_binary = which(gem_binary)
if not gem_binary:
raise Exception('\n\nERROR: GEM binary not found, install it from:'
'\nhttps://sourceforge.net/projects/gemlibrary/files/gem-library/Binary%20pre-release%202/'
'\n - Download the GEM-binaries-Linux-x86_64-core_i3 if'
'have a recent computer, the '
'GEM-binaries-Linux-x86_64-core_2 otherwise\n - '
'Uncompress with "tar xjvf GEM-binaries-xxx.tbz2"\n - '
'Copy the binary gem-mapper to /usr/local/bin/ for '
'example (somewhere in your PATH).\n\nNOTE: GEM does '
'not provide any binary for MAC-OS.')
try:
out, err = Popen([gem_binary,'--version'], stdout=PIPE, stderr=STDOUT,
universal_newlines=True).communicate()
gem_version = int(out[1])
except ValueError as e:
gem_version = 2
print('Falling to gem v2')
if mapper_params and isinstance(mapper_params, dict):
kwargs.update(mapper_params)
# create directories
for rep in [temp_dir, out_map_dir]:
mkdir(rep)
# check space
fspace = int(get_free_space_mb(temp_dir, div=3))
if fspace < 200:
warn('WARNING: only %d Gb left on tmp_dir: %s\n' % (fspace, temp_dir))
# iterative mapping
base_name = os.path.split(fastq_path)[-1].replace('.gz', '')
base_name = '.'.join(base_name.split('.')[:-1])
input_reads = fastq_path
if windows is None:
light_storage = True
windows = (None, )
elif isinstance(windows[0], int):
# if windows starts at zero we do not need to store all the sequence
# otherwise we need it because sequence can be trimmed two times
# in fragment based mapping
light_storage = True if not windows[0] else False
windows = [tuple(windows)]
else:
# ensure that each element is a tuple, not a list
windows = [tuple(win) for win in windows]
# in this case we will need to keep the information about original
# sequence at any point, light storage is thus not possible.
light_storage = False
for win in windows:
# Prepare the FASTQ file and iterate over them
curr_map, counter = transform_fastq(
input_reads, mkstemp(prefix=base_name + '_', dir=temp_dir)[1],
fastq=is_fastq(input_reads),
min_seq_len=min_seq_len, trim=win, skip=skip, nthreads=nthreads,
light_storage=light_storage)
# clean
if input_reads != fastq_path and clean:
print(' x removing original input %s' % input_reads)
os.system('rm -f %s' % (input_reads))
# First mapping, full length
if not win:
beg, end = 1, 'end'
else:
beg, end = win
out_map_path = curr_map + '_full_%s-%s%s.map' % (beg, end, suffix)
if end:
print('Mapping reads in window %s-%s%s...' % (beg, end, suffix))
else:
print('Mapping full reads...', curr_map)
if not skip:
if mapper == 'gem':
_gem_mapping(mapper_index_path, curr_map, out_map_path,
gem_binary=gem_binary, gem_version=gem_version,
gem_params=mapper_params, **kwargs)
# parse map file to extract not uniquely mapped reads
print('Parsing result...')
if gem_version >= 3:
_sam_filter(out_map_path, curr_map,
curr_map + '_filt_%s-%s%s.map' % (beg, end, suffix),
os.path.join(out_map_dir,
base_name + '_full_%s-%s%s.map' % (beg, end, suffix)))
else:
_gem_filter(out_map_path,
curr_map + '_filt_%s-%s%s.map' % (beg, end, suffix),
os.path.join(out_map_dir,
base_name + '_full_%s-%s%s.map' % (
beg, end, suffix)))
elif mapper == 'bowtie2' or mapper == 'hisat2':
_bowtie2_mapping(mapper_index_path, curr_map, out_map_path,
bowtie2_binary=(mapper_binary if mapper_binary else mapper),
bowtie2_params=mapper_params, **kwargs)
# parse map file to extract not uniquely mapped reads
print('Parsing result...')
_sam_filter(out_map_path, curr_map,
curr_map + '_filt_%s-%s%s.map' % (beg, end, suffix),
os.path.join(out_map_dir,
base_name + '_full_%s-%s%s.map' % (beg, end, suffix)))
else:
raise Exception('ERROR: unknown mapper.')
# clean
if clean:
print(' x removing %s input %s' % (mapper.upper(),curr_map))
os.system('rm -f %s' % (curr_map))
print(' x removing map %s' % out_map_path)
os.system('rm -f %s' % (out_map_path))
# for next round, we will use remaining unmapped reads
input_reads = curr_map + '_filt_%s-%s%s.map' % (beg, end, suffix)
outfiles.append(
(os.path.join(out_map_dir,
base_name + '_full_%s-%s%s.map' % (beg, end, suffix)),
counter))
# map again splitting unmapped reads into RE fragments
# (no need to trim this time)
if frag_map:
if not r_enz:
raise Exception('ERROR: need enzyme name to fragment.')
frag_map, counter = transform_fastq(
input_reads, mkstemp(prefix=base_name + '_', dir=temp_dir)[1],
min_seq_len=min_seq_len, trim=win, fastq=False, r_enz=r_enz,
add_site=add_site, skip=skip, nthreads=nthreads,
light_storage=light_storage)
# clean
if clean:
print(' x removing pre-%s input %s' % (mapper.upper(),input_reads))
os.system('rm -f %s' % (input_reads))
if not win:
beg, end = 1, 'end'
else:
beg, end = win
out_map_path = frag_map + '_frag_%s-%s%s.map' % (beg, end, suffix)
if not skip:
if mapper == 'gem':
print('Mapping fragments of remaining reads...')
_gem_mapping(mapper_index_path, frag_map, out_map_path,
gem_binary=gem_binary, gem_version=gem_version,
**kwargs)
print('Parsing result...')
# check if output is sam format for gem3
if gem_version >= 3:
_sam_filter(out_map_path, frag_map,
curr_map + '_fail%s.map' % (suffix),
os.path.join(out_map_dir,
base_name + '_frag_%s-%s%s.map' % (beg, end, suffix)))
else:
_gem_filter(out_map_path, curr_map + '_fail%s.map' % (suffix),
os.path.join(out_map_dir,
base_name + '_frag_%s-%s%s.map' % (beg, end, suffix)))
elif mapper == 'bowtie2' or mapper == 'hisat2':
print('Mapping fragments of remaining reads...')
_bowtie2_mapping(mapper_index_path, frag_map, out_map_path,
bowtie2_binary=(mapper_binary if mapper_binary else mapper),
bowtie2_params=mapper_params, **kwargs)
print('Parsing result...')
_sam_filter(out_map_path, frag_map,
curr_map + '_fail%s.map' % (suffix),
os.path.join(out_map_dir,
base_name + '_frag_%s-%s%s.map' % (beg, end, suffix)))
else:
raise Exception('ERROR: unknown mapper.')
# clean
if clean:
print(' x removing %s input %s' % (mapper.upper(),frag_map))
os.system('rm -f %s' % (frag_map))
print(' x removing failed to map ' + curr_map + '_fail%s.map' % (suffix))
os.system('rm -f %s' % (curr_map + '_fail%s.map' % (suffix)))
print(' x removing tmp mapped %s' % out_map_path)
os.system('rm -f %s' % (out_map_path))
outfiles.append((os.path.join(out_map_dir,
base_name + '_frag_%s-%s%s.map' % (beg, end, suffix)),
counter))
if clean:
os.system('rm -rf %s' % (temp_dir))
if get_nread:
return outfiles
return [out for out, _ in outfiles]
def fast_fragment_mapping(mapper_index_path, fastq_path1, fastq_path2, r_enz,
genome_seq, out_map, clean=True, get_nread=False,
mapper_binary=None, mapper_params=None,
samtools = 'samtools', **kwargs):
"""
Maps FASTQ reads to an indexed reference genome with the knowledge of
the restriction enzyme used (fragment-based mapping).
:param mapper_index_path: path to index file created from a reference genome
using gem-index tool, bowtie2-build or hisat2-build
:param fastq_path1: PATH to FASTQ file of read 1, either compressed or not.
:param fastq_path2: PATH to FASTQ file of read 2, either compressed or not.
:param out_map_dir: path to outfile tab separated format containing mapped
read information.
:param r_enz: name of the restriction enzyme used in the experiment e.g.
HindIII.
:param genome_seq: a dictionary generated by :func:`pyatdbit.parser.genome_parser.parse_fasta`.
containing the genomic sequence
:param False clean: remove intermediate files created in temp_dir
:param False get_nread: returns a list of lists where each element contains
a path and the number of reads processed
:param 4 nthreads: number of threads to use for mapping (number of CPUs)
:param /tmp temp_dir: important to change. Intermediate FASTQ files will be
written there.
:param gem-mapper mapper_binary: path to the binary mapper
:param None mapper_params: extra parameters for the mapper
:param samtools samtools: path to samtools binary.
:returns: outfile with the intersected read pairs
"""
suffix = kwargs.get('suffix', '')
suffix = ('_' * (suffix != '')) + suffix
nthreads = kwargs.get('nthreads', 8)
samtools = which(samtools)
# check out folder
if not os.path.isdir(os.path.dirname(os.path.abspath(out_map))):
raise Exception('\n\nERROR: Path to store the output does not exist.\n')
temp_dir = os.path.abspath(os.path.expanduser(
kwargs.get('temp_dir', gettempdir())))
gem_version = None
# check that we have the GEM binary:
gem_binary = mapper_binary or 'gem-mapper'
gem_binary = which(gem_binary)
if not gem_binary:
raise Exception('\n\nERROR: GEM v3 binary not found, install it from:'
'\nhttps://github.com/smarco/gem3-mapper'
'Copy the binary gem-mapper to /usr/local/bin/ for '
'example (somewhere in your PATH).\n')
try:
out, err = Popen([gem_binary,'--version'], stdout=PIPE, stderr=STDOUT,
universal_newlines=True).communicate()
gem_version = int(out[1])
except ValueError as e:
gem_version = 2
print('Falling to gem v2')
if gem_version < 3:
raise Exception('\n\nERROR: GEM v3 binary not found, install it from:'
'\nhttps://github.com/smarco/gem3-mapper'
'Copy the binary gem-mapper to /usr/local/bin/ for '
'example (somewhere in your PATH).\n')
if mapper_params:
kwargs.update(mapper_params)
# create directories
for rep in [temp_dir]:
mkdir(rep)
# check space
fspace = int(get_free_space_mb(temp_dir, div=3))
if fspace < 200:
warn('WARNING: only %d Gb left on tmp_dir: %s\n' % (fspace, temp_dir))
# iterative mapping
base_name1 = os.path.split(fastq_path1)[-1].replace('.gz', '')
base_name1 = '.'.join(base_name1.split('.')[:-1])
curr_map1, _ = transform_fastq(
fastq_path1, mkstemp(prefix=base_name1 + '_', dir=temp_dir)[1],
fastq=is_fastq(fastq_path1), nthreads=nthreads, light_storage=True)
base_name2 = os.path.split(fastq_path2)[-1].replace('.gz', '')
base_name2 = '.'.join(base_name2.split('.')[:-1])
curr_map2, count_fastq = transform_fastq(
fastq_path2, mkstemp(prefix=base_name2 + '_', dir=temp_dir)[1],
fastq=is_fastq(fastq_path1), nthreads=nthreads, light_storage=True)
out_map_path = curr_map1 + '_frag%s.map' % (suffix)
print('Mapping fragments of remaining reads...')
_gem_mapping(mapper_index_path, curr_map1, out_map_path,fastq_path2=curr_map2,
r_enz=r_enz, gem_binary=gem_binary, gem_version=gem_version, **kwargs)
# clean
if clean:
print(' x removing GEM 3 input %s' % (curr_map1))
os.system('rm -f %s' % (curr_map1))
print(' x removing GEM 3 input %s' % (curr_map2))
os.system('rm -f %s' % (curr_map2))
#sort sam file
os.system(samtools + ' sort -n -O SAM -@ %d -T %s -o %s %s'
% (nthreads, out_map_path, out_map_path, out_map_path))
genome_lengths = dict((crm, len(genome_seq[crm])) for crm in genome_seq)
frag_chunk = kwargs.get('frag_chunk', 100000)
frags = map_re_sites(r_enz, genome_seq, frag_chunk=frag_chunk)
if samtools and nthreads > 1:
print('Splitting sam file')
# headers
for i in range(nthreads):
os.system(samtools + ' view -H -O SAM %s > "%s_%d"'
% (out_map_path, out_map_path, (i+1)))
chunk_lines = int((count_fastq*2.3)/nthreads) # estimate lines in sam with reads and frags
os.system(samtools + ''' view -O SAM %s | awk -v n=%d -v FS="\\t" '
BEGIN { part=0; line=n }
{ if( line>=n && $1!=last_read ) {part++; line=1; print $0 >> "%s_"part }
else { print $0 >> "%s_"part; line++; }
last_read = $1;
}'
''' % (out_map_path, chunk_lines, out_map_path, out_map_path))
if clean:
print(' x removing tmp mapped %s' % out_map_path)
os.system('rm -f %s' % (out_map_path))
print('Parsing results...')
kwargs['nthreads'] = 1
procs = []
pool = mu.Pool(nthreads)
for i in range(nthreads):
frags_shared = copy.deepcopy(frags)
procs.append(pool.apply_async(
parse_gem_3c, args=('%s_%d' % (out_map_path,(i+1)),
'%s_parsed_%d' % (out_map_path,(i+1)),
copy.deepcopy(genome_lengths), frags_shared,
False, True), kwds=kwargs))
#results.append('%s_parsed_%d' % (out_map_path,(i+1)))
pool.close()
pool.join()
results = [proc.get() for proc in procs if proc.get()]
if clean:
for i in range(nthreads):
print(' x removing tmp mapped %s_%d' % (out_map_path,(i+1)))
os.system('rm -f %s_%d' % (out_map_path,(i+1)))
#Final sort and merge
nround = 0
while len(results) > 1:
nround += 1
num_procs = min(nthreads,int(len(results)/2))
pool = mu.Pool(num_procs)
procs = [pool.apply_async(
merge_sort,
(results.pop(0), results.pop(0), out_map_path+'_%d' % nround, i, True)
) for i in range(num_procs)]
pool.close()
pool.join()
results = [proc.get() for proc in procs if proc.get()]
map_out = open(out_map, 'w')
tmp_reads_fh = open(results[0],'r')
for crm in genome_seq:
map_out.write('# CRM %s\t%d\n' % (crm, len(genome_seq[crm])))
for read_line in tmp_reads_fh:
read = read_line.split('\t')
map_out.write('\t'.join([read[0]]+read[2:8]+read[9:]))
map_out.close()
if clean:
print(' x removing tmp mapped %s' % results[0])
os.system('rm -f %s' % (results[0]))
else:
print('Parsing result...')
parse_gem_3c(out_map_path, out_map, genome_lengths, frags, verbose=False,
tmp_format=False, **kwargs)
# clean
if clean:
print(' x removing tmp mapped %s' % out_map_path)
os.system('rm -f %s' % (out_map_path))
if get_nread:
return [(out_map,count_fastq)]
return out_map
