def alignment_to_contacts(<EOL>sam_merged,<EOL>assembly,<EOL>output_dir,<EOL>output_file_network=DEFAULT_NETWORK_FILE_NAME,<EOL>output_file_chunk_data=DEFAULT_CHUNK_DATA_FILE_NAME,<EOL>parameters=DEFAULT_PARAMETERS,<EOL>):
all_contacts = collections.Counter()<EOL>all_chunks = collections.Counter()<EOL>chunk_size = int(parameters["<STR_LIT>"])<EOL>mapq_threshold = int(parameters["<STR_LIT>"])<EOL>size_chunk_threshold = int(parameters["<STR_LIT>"])<EOL>read_size = int(parameters["<STR_LIT>"])<EOL>self_contacts = parameters["<STR_LIT>"]<EOL>normalized = parameters["<STR_LIT>"]<EOL>logger.info("<STR_LIT>")<EOL>chunk_complete_data = dict()<EOL>global_id = <NUM_LIT:1><EOL>for record in SeqIO.parse(assembly, "<STR_LIT>"):<EOL><INDENT>length = len(record.seq)<EOL>n_chunks = length // chunk_size<EOL>n_chunks += (length % chunk_size) >= size_chunk_threshold<EOL>for i in range(n_chunks):<EOL><INDENT>if (i + <NUM_LIT:1>) * chunk_size <= length:<EOL><INDENT>size = chunk_size<EOL><DEDENT>else:<EOL><INDENT>size = length % chunk_size<EOL><DEDENT>chunk_name = "<STR_LIT>".format(record.id, i)<EOL>chunk_complete_data[chunk_name] = {<EOL>"<STR_LIT:id>": global_id,<EOL>"<STR_LIT>": <NUM_LIT:0>,<EOL>"<STR_LIT:size>": size,<EOL>"<STR_LIT>": <NUM_LIT:0>,<EOL>}<EOL>global_id += <NUM_LIT:1><EOL><DEDENT><DEDENT>logger.info("<STR_LIT>")<EOL>current_read = None<EOL>with pysam.AlignmentFile(sam_merged, "<STR_LIT:rb>") as alignment_merged_handle:<EOL><INDENT>names = alignment_merged_handle.references<EOL>lengths = alignment_merged_handle.lengths<EOL>names_and_lengths = {<EOL>name: length for name, length in zip(names, lengths)<EOL>}<EOL>logger.info("<STR_LIT>")<EOL>while "<STR_LIT>":<EOL><INDENT>try:<EOL><INDENT>my_read = next(alignment_merged_handle)<EOL>if current_read is None:<EOL><INDENT>current_read = my_read<EOL>continue<EOL><DEDENT>elif current_read.query_name != my_read.query_name:<EOL><INDENT>current_read = my_read<EOL>continue<EOL><DEDENT>read_forward, read_reverse = current_read, my_read<EOL><DEDENT>except StopIteration:<EOL><INDENT>break<EOL><DEDENT>read_name_forward = read_forward.query_name<EOL>read_name_reverse = read_reverse.query_name<EOL>flag_forward, flag_reverse = read_forward.flag, read_reverse.flag<EOL>try:<EOL><INDENT>assert read_name_forward == read_name_reverse<EOL><DEDENT>except AssertionError:<EOL><INDENT>logger.error(<EOL>"<STR_LIT>" "<STR_LIT>",<EOL>read_name_forward,<EOL>read_name_reverse,<EOL>)<EOL>raise<EOL><DEDENT>def is_unmapped(flag):<EOL><INDENT>return np.base_repr(flag, padding=<NUM_LIT:3>)[-<NUM_LIT:3>] == "<STR_LIT:1>"<EOL><DEDENT>if is_unmapped(flag_forward) or is_unmapped(flag_reverse):<EOL><INDENT>continue<EOL><DEDENT>contig_name_forward = read_forward.reference_name<EOL>contig_name_reverse = read_reverse.reference_name<EOL>len_contig_for = names_and_lengths[contig_name_forward]<EOL>len_contig_rev = names_and_lengths[contig_name_reverse]<EOL>position_forward = read_forward.reference_start<EOL>position_reverse = read_reverse.reference_start<EOL>mapq_forward = read_forward.mapping_quality<EOL>mapq_reverse = read_reverse.mapping_quality<EOL>mapq_test = min(mapq_forward, mapq_reverse) > mapq_threshold<EOL>min_length = min(len_contig_for, len_contig_rev)<EOL>length_test = min_length > size_chunk_threshold<EOL>if position_forward < chunk_size * (len_contig_for // chunk_size):<EOL><INDENT>current_chunk_forward_size = chunk_size<EOL><DEDENT>else:<EOL><INDENT>current_chunk_forward_size = len_contig_for % chunk_size<EOL><DEDENT>if position_reverse < chunk_size * (len_contig_rev // chunk_size):<EOL><INDENT>current_chunk_reverse_size = chunk_size<EOL><DEDENT>else:<EOL><INDENT>current_chunk_reverse_size = len_contig_rev % chunk_size<EOL><DEDENT>min_chunk_size = min(<EOL>current_chunk_forward_size, current_chunk_reverse_size<EOL>)<EOL>chunk_test = min_chunk_size >= size_chunk_threshold<EOL>if mapq_test and length_test and chunk_test:<EOL><INDENT>chunk_forward = position_forward // chunk_size<EOL>chunk_reverse = position_reverse // chunk_size<EOL>chunk_name_forward = "<STR_LIT>".format(<EOL>contig_name_forward, chunk_forward<EOL>)<EOL>chunk_name_reverse = "<STR_LIT>".format(<EOL>contig_name_reverse, chunk_reverse<EOL>)<EOL>if self_contacts or chunk_name_forward != chunk_name_reverse:<EOL><INDENT>contact = tuple(<EOL>sorted((chunk_name_forward, chunk_name_reverse))<EOL>)<EOL>all_contacts[contact] += <NUM_LIT:1><EOL>chunk_key_forward = (<EOL>chunk_name_forward,<EOL>current_chunk_forward_size,<EOL>)<EOL>all_chunks[chunk_key_forward] += <NUM_LIT:1><EOL>chunk_key_reverse = (<EOL>chunk_name_reverse,<EOL>current_chunk_reverse_size,<EOL>)<EOL>all_chunks[chunk_key_reverse] += <NUM_LIT:1><EOL><DEDENT><DEDENT><DEDENT><DEDENT>logger.info("<STR_LIT>")<EOL>output_chunk_data_path = os.path.join(output_dir, output_file_chunk_data)<EOL>with open(output_chunk_data_path, "<STR_LIT:w>") as chunk_data_file_handle:<EOL><INDENT>for name in sorted(chunk_complete_data.keys()):<EOL><INDENT>chunk_data = chunk_complete_data[name]<EOL>size = chunk_data["<STR_LIT:size>"]<EOL>chunk = (name, chunk_data["<STR_LIT:size>"])<EOL>hit = all_chunks[chunk]<EOL>coverage = hit * read_size * <NUM_LIT:1.0> / size<EOL>try:<EOL><INDENT>chunk_complete_data[name]["<STR_LIT>"] = hit<EOL>chunk_complete_data[name]["<STR_LIT>"] = coverage<EOL><DEDENT>except KeyError:<EOL><INDENT>logger.error(<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>)<EOL>raise<EOL><DEDENT>idx = chunk_complete_data[name]["<STR_LIT:id>"]<EOL>line = "<STR_LIT>".format(<EOL>idx, name, hit, size, coverage<EOL>)<EOL>chunk_data_file_handle.write(line)<EOL><DEDENT><DEDENT>logger.info("<STR_LIT>")<EOL>output_network_path = os.path.join(output_dir, output_file_network)<EOL>with open(output_network_path, "<STR_LIT:w>") as network_file_handle:<EOL><INDENT>for chunks in sorted(all_contacts.keys()):<EOL><INDENT>chunk_name1, chunk_name2 = chunks<EOL>contact_count = all_contacts[chunks]<EOL>if normalized:<EOL><INDENT>coverage1 = chunk_complete_data[chunk_name1]["<STR_LIT>"]<EOL>coverage2 = chunk_complete_data[chunk_name2]["<STR_LIT>"]<EOL>mean_coverage = np.sqrt(coverage1 * coverage2)<EOL>effective_count = contact_count * <NUM_LIT:1.0> / mean_coverage<EOL><DEDENT>else:<EOL><INDENT>effective_count = contact_count<EOL><DEDENT>try:<EOL><INDENT>idx1 = chunk_complete_data[chunk_name1]["<STR_LIT:id>"]<EOL>idx2 = chunk_complete_data[chunk_name2]["<STR_LIT:id>"]<EOL>line = "<STR_LIT>".format(idx1, idx2, effective_count)<EOL>network_file_handle.write(line)<EOL><DEDENT>except KeyError as e:<EOL><INDENT>logger.warning("<STR_LIT>", e)<EOL><DEDENT><DEDENT><DEDENT>return chunk_complete_data, all_contacts<EOL>
Generates a network file (in edgelist form) from an
    alignment in sam or bam format. Contigs are virtually split into
    'chunks' of nearly fixed size (by default between 500 and 1000 bp)
    to reduce size bias. The chunks are the network nodes and the edges
    are the contact counts.

    The network is in a strict barebone form so that it can be reused and
    imported quickly into other applications etc. Verbose information about
    every single node in the network is written on a 'chunk data' file,
    by default called 'idx_contig_hit_size_cov.txt'

    Parameters
    ----------
    sam_merged : file, str or pathlib.Path
        The alignment file in SAM/BAM format to be processed.
    assembly : file, str or pathlib.Path
        The initial assembly acting as the alignment file's reference genome.
    output_dir : str or pathlib.Path
        The output directory to write the network and chunk data into.
    output_dir_file_network : str or pathlib.Path, optional
        The specific file name for the output network file. Default is
        network.txt
    output_file_chunk_data : str or pathlib.Path, optional
        The specific file name for the output chunk data file. Default is
        idx_contig_hit_size_cov.txt
    parameters : dict, optional
        A dictionary of parameters for converting the alignment file into a
        network. These are:
        -size_chunk_threshold: the size (in bp) under which chunks are
        discarded. Default is 500.
        -mapq_threshold: the mapping quality under which alignments are
        discarded. Default is 10.
        -chunk_size: the default chunk size (in bp) when applicable, save
        smaller contigs or tail-ends. Default is 1000.
        -read_size: the size of reads used for mapping. Default is 65.
        -self_contacts: whether to count alignments between a chunk and
        itself. Default is False.
        -normalized: whether to normalize contacts by their coverage.
        Default is False.

    Returns
    -------
    chunk_complete_data : dict
        A dictionary where the keys are chunks in (contig, position) form and
        the values are their id, name, total contact count, size and coverage.
    all_contacts : dict
        A counter dictionary where the keys are chunk pairs and the values are
        their contact count.