RIC-seq-based Prediction of Complementary RNA-RNA Contacts

public 1yr ago 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output, operation

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

Prediction of complementary regions between RNA-RNA contacts derived from RIC-seq data. Developed by Sergei Margasyuk (smargasyuk@gmail.com) and Dmitri Pervouchine (pervouchine@gmail.com).

Description

This package contains a pipeline for

Code Snippets

    shell: """
cat {input} |\
awk -v 'OFS=\t' '($10>{params.min_support})&&($22>{params.min_support})' |\
awk -v 'OFS=\t' '{{print $1,$15,$2,$13,$6,$17,$7,$8,$9,$10,$19,$20,$21,$22,$11,$12,$25}}' > {output}
"""

SnakeMake From line 7 of rules/fold.smk

    shell: """
paste <(cut -f1-6 {input.pairs}) {input.pairs} |\
awk -v 'OFS=\t' '($3-$2<={params.max_radius}) && ($3-$2>={params.min_radius})' |\
sort-bed - |\
bedops -n 1 - <(sort-bed {input.rmsk}) |\
cut -f4- |\
awk -v 'OFS=\t' '($3-$2<={params.max_radius}) && ($3-$2>={params.min_radius})' |\
sort-bed - |\
bedops -n 1 - <(sort-bed {input.rmsk}) |\
cut -f 4- > {output}
"""

SnakeMake Bedops From line 23 of rules/fold.smk

    shell: """
cat {input.pairs} |\
awk -v 'OFS=\t' '{{print $1,$2,$3,$17,".",$15}}' |\
bedtools getfasta -s -bed stdin -fi {input.genome_fasta} > {output.seq_L}

cat {input.pairs} |\
awk -v 'OFS=\t' '{{print $4,$5,$6,$17,".",$15}}' |\
bedtools getfasta -s -bed stdin -fi {input.genome_fasta} > {output.seq_R}
"""

SnakeMake BEDTools From line 43 of rules/fold.smk

    shell: """
paste  <(grep -v '>' {input.seq_L}) <(grep -v '>' {input.seq_R}) |\
grep -v "^$" |\
python workflow/scripts/PrePH/src/fold2.py -e {params.max_energy} -a 8 -u False -j{threads} \
> {output}
"""

SnakeMake From line 62 of rules/fold.smk

    shell: """
paste {input.pairs} {input.preph} <(grep -v '>' {input.seq_L}) <(grep -v '>' {input.seq_R}) |\
awk -v 'OFS=\t' '{{gsub(/[\(\)]/, "", $22)}}1' |\
awk -v 'OFS=\t' 'NF==24' |\
awk -v 'OFS=\t' '{{print $1,$2,$6,$17,$10+$14,$15,$2,$6,"0,0,0",2,$3-$2","$6-$5,0","$5-$2,$0}}'|\
sort-bed - > {output}
"""

SnakeMake From line 77 of rules/fold.smk

run:
    t0 = pd.read_table(str(input), header=None)
    t3 = postprocess_preph_table(t0)
    t3.to_csv(str(output.tsv), index=None, sep='\t')

SnakeMake From line 90 of rules/fold.smk

    shell: """
mkdir -p $(dirname {output.starch})
gtfToGenePred {input} {output.gp1} -geneNameAsName2 -ignoreGroupsWithoutExons -genePredExt
awk '{{print $12"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12"\t"$13"\t"$14"\t"$15"\t"}}' {output.gp1} > {output.gp2}
genePredToBed {output.gp2} {output.bed1}
sort-bed {output.bed1} | starch -  > {output.starch}
"""

SnakeMake gtftogenepred From line 9 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
unpigz -c {input} |\
awk -v 'OFS=\t' '{{print $0,"{wildcards.project}","{wildcards.sample}","{wildcards.t}"}}' |\
pigz - > {output}    
"""

SnakeMake From line 23 of rules/nested.smk

    shell: """
unpigz -c {input} | pigz - > {output}    
"""

SnakeMake From line 34 of rules/nested.smk

    shell: """
unpigz -c {input} |\
awk -v 'OFS=\t' '$1==$5' |\
awk -v 'OFS=\t' '{{if ($2>$6){{s1=$2;s2=$3;s3=$4;$2=$6;$3=$7;$4=$8;$6=s1;$7=s2;$8=s3}};print}}' |\
awk -v 'OFS=\t' '$6-$2>{params.min_spread}' |\
pigz > {output}
"""

SnakeMake From line 44 of rules/nested.smk

    shell: """
unpigz -c {input} |\
cut -f1,2,3,5,6,7,11 |\
sort --buffer-size={resources.mem_mb}M --parallel={threads} |\
uniq -c |\
sed -E 's/^ *//; s/ /\t/' |\
awk -v 'OFS=\t' '{{print $2,$3,$4,$5,$6,$7,$8,$1}}' |\
pigz > {output}
"""

SnakeMake From line 59 of rules/nested.smk

    shell: """
unpigz -c {input} |\
awk -v 'OFS=\t' '{{print $0,$0}}' |\
sort-bed --tmpdir {resources.tmpdir} --max-mem {resources.mem_mb}M - |\
bedtools cluster -d {params.clustering_distance} -i stdin |\
cut -f 4- |\
sort-bed --tmpdir {resources.tmpdir} --max-mem {resources.mem_mb}M - |\
bedtools cluster -d {params.clustering_distance} -i stdin |\
cut -f 6- |\
awk -v 'OFS=\t' '{{print $0,"id_"NR}}' |\
sort --parallel={threads}  --buffer-size={resources.mem_mb}M -k9,9 -k10,10 |\
datamash groupby 9,10 first 1 min 2 max 3 first 4 min 5 max 6 unique 7 countunique 7 countunique 11 sum 8 |\
cut -f 3- |\
pigz > {output}
"""

SnakeMake BEDTools From line 78 of rules/nested.smk

    shell: """
unpigz -c {input} |\
awk -v 'OFS=\t' '{{print $1,$2,$6,".", ".", "+",$0}}' |\
sort-bed --tmpdir {resources.tmpdir} --max-mem {resources.mem_mb}M - |\
starch - > {output}
"""

SnakeMake From line 100 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
bedmap --chrom {wildcards.chrom} --fraction-ref 1 --echo --echo-map --delim "\t" --multidelim "\t" --skip-unmapped {input.contacts} {input.annotation} |\
awk -v "OFS=\t" '{{$17=$22;$18=$20; print $0}}' | cut -f7-18 |\
awk -v "OFS=\t" '{{print $0,"id_"NR}}' |\
pigz > {output}
"""

SnakeMake From line 114 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
unpigz -c {input} |\
awk -v 'OFS=\t' '{{print $1,$2,$3,$13}}' |\
sort-bed - > {output}  
""" 

SnakeMake From line 126 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
unpigz -c {input} |\
awk -v 'OFS=\t' '{{print $4,$5,$6,$13}}' |\
sort-bed - > {output}  
""" 

SnakeMake From line 137 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
cat {input} |\
awk -v 'OFS=\t' '{{$2=$3+{params.min_radius};$3=$3+{params.max_radius};print $0}}' | sort-bed -  > {output}  
""" 

SnakeMake From line 151 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
cat {input} |\
awk -v 'OFS=\t' '{{$3=$2-{params.min_radius};$2=$2-{params.max_radius};print $0}}' | sort-bed -  > {output}  
""" 

SnakeMake From line 164 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
intersectBed -wa -wb -sorted \
-a {input.handles} \
-b {input.windows} |\
awk -v 'OFS=\t' '{{print "{wildcards.chrom}_"$4"_"$8,$4,$8}}' |\
sort --parallel={threads}  --buffer-size={resources.mem_mb}M  -T $TEMPDIR |\
pigz - > {output}
"""

SnakeMake BEDTools From line 179 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
comm -12 <(unpigz -c {input.left}) <(unpigz -c {input.right}) | awk -v 'OFS=\t' '$2!=$3' > {output}
"""

SnakeMake From line 195 of rules/nested.smk

    shell: """
mkdir -p $(dirname {output})
join -1 2 -2 13 -t$'\t' <(sort -k2,2 {input.pairs}) <(unpigz -c {input.contacts} | sort -k13,13) > {output.l}
join -1 3 -2 13 -t$'\t' <(sort -k3,3 {output.l})    <(unpigz -c {input.contacts} | sort -k13,13) | cut -f 4- > {output.r}
"""

SnakeMake From line 208 of rules/nested.smk

    shell: """
awk -v 'OFS=\t' '{{print $0,"id_"NR}}' {input} > {output}    
"""

SnakeMake From line 218 of rules/nested.smk

run:
    for f1, f2 in zip(input["pipeline_out"], input["ref_out"]):
        if not filecmp.cmp(f1, f2, shallow = False):
            print(f"Files {f1} and {f2} are different")
            raise ValueError

SnakeMake From line 8 of rules/test.smk

from numpy import (
    argmin,
    unravel_index,
    full,
    empty,
    load,
    set_printoptions,
    argwhere,
    argsort,
)
from math import ceil
import binascii, itertools, sys, getopt, os
from functools import partial
from sys import getsizeof
import sys
from multiprocessing import Pool
import argparse

inf = float("inf")
os.chdir(os.path.dirname(os.path.abspath(__file__)))

# Dictionary for nts (used in 1x1, 2x1, 2x2 loops in last 2 dims)
Dic_nt = {"@": 0, "A": 1, "C": 2, "G": 3, "T": 4}
# Dictionary of basepairs (used in 1x1, 2x1, 2x2 loops in first 2 dims)
Dic_bp = {"CG": 0, "GC": 1, "GT": 2, "TG": 3, "AT": 4, "TA": 5}
stacking_matrix = load("../lib//stacking_matrix.npy")
bulge_list = load("../lib/bulge_list.npy")
intl11_matrix = load("../lib/intl11_matrix.npy")
intl12_matrix = load("../lib/intl12_matrix.npy")
intl22_matrix = load("../lib/intl22_matrix.npy")

# Adding for long bulges
TerminalAU = 50


def Seq_to_bin(seq):  # Uses bitwise shift to make bin from sequence
    Dict = {"A": 0b0, "T": 0b10, "G": 0b11, "C": 0b1}
    s = 0
    for char in seq:
        try:
            s = s << 2 | Dict[char]
        except KeyError:
            s = False
    return s


def Index_seq(seq, k):  # Uses bitwise shift to divide bin into kmers
    seq_bin = Seq_to_bin(seq)
    if seq_bin is False:
        return False
    else:
        seq_indxd_tmp = []
        mask = 2 ** (k * 2) - 1
        for i in range(len(seq) - (k - 1)):
            seq_indxd_tmp.append(mask & seq_bin)
            seq_bin >>= 2
        return seq_indxd_tmp


def Initiate_with_kmers(
    seq, seq_compl, seq_indxd_tmp, seq_compl_indxd_tmp, kmers_stacking_matrix, k
):
    min_energy = end_pos_i = end_pos_j = start_pos_i = start_pos_j = 0
    # seq_indxd_tmp = Index_seq(seq, k=k)
    # seq_compl_indxd_tmp = Index_seq(seq_compl[::-1], k=k)
    if (seq_indxd_tmp == False) | (seq_compl_indxd_tmp == False):
        return (
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
            False,
        )
    else:
        seq_length = len(seq)
        seq = seq + "$" * (k + 2)  # (horizontally) _
        seq_compl = "$" * (k + 2) + seq_compl  # (vertically) |
        seq_compl_length = len(seq_compl)
        D = full((len(seq_compl), len(seq)), inf)  # distance matrix
        zero_coords = empty((), dtype=object)
        zero_coords[()] = (0, 0)
        B = full(
            (len(seq_compl), len(seq)), zero_coords, dtype=object
        )  # backtracker matrix
        S = empty(
            [len(seq_compl), len(seq)], dtype="S" + str(len(seq) + len(seq_compl))
        )  # dot bracket structure matrix
        for I, kmer_i in enumerate(seq_compl_indxd_tmp):
            i = seq_compl_length - I - k
            for J, kmer_j in enumerate(seq_indxd_tmp):
                j = seq_length - J - 1
                D[i][j] = kmers_stacking_matrix[kmer_j, kmer_i]
                if D[i][j] != inf:
                    B[i][j] = (i + k - 1, j - k + 1)
                    S[i][j] = "(" * k + ")" * k
                    if D[i][j] < min_energy:
                        min_energy = D[i][j]
                        end_pos_i = i
                        end_pos_j = j
                        start_pos_i = i + k - 1
                        start_pos_j = j - k + 1
        return (
            D,
            B,
            S,
            min_energy,
            end_pos_i,
            end_pos_j,
            start_pos_i,
            start_pos_j,
            seq_indxd_tmp,
            seq_compl_indxd_tmp,
            seq_length,
            seq_compl_length,
            seq,
            seq_compl,
        )


def End_coords(argmin_, i, j, old_end_pos_i, old_end_pos_j, k=3):
    Dict_end_coords = {
        1: (-1, +1),
        2: (-k, +1 + k),
        3: (-1 - k, k),
        4: (-k, +2 + k),
        5: (-2 - k, k),
        6: (-1 - k, +1 + k),
        7: (-1 - k, +2 + k),
        8: (-2 - k, +1 + k),
        9: (-2 - k, +2 + k),
    }
    if argmin_ == 0:
        return (old_end_pos_i, old_end_pos_j)
    else:
        add_i, add_j = Dict_end_coords[argmin_]
        return (i + add_i, j + add_j)


def Start_coords(argmin_, backtrack, old_start_pos_i, old_start_pos_j):
    if argmin_ == 0:
        return (old_start_pos_i, old_start_pos_j)
    else:
        return (backtrack[0], backtrack[1])


def Backtrack(argmin_, old_coords, new_coords):
    if argmin_ == 0:
        return (0, 0)
    elif argmin_ == 1:
        return old_coords
    else:
        return new_coords


def Check_ranges_overlap(x1, x2, y1, y2):
    return (x1 < y2) & (y1 < x2)


# segment1[start[i,j],end[i,j]], segment2 - part of the square [start[i,j],end[i,j]]
def Check_segments_intersection(segment1, segment2, slope):
    a1 = float(-segment1[0][0] + segment1[1][0]) / (segment1[0][1] - segment1[1][1])
    b1 = -segment1[0][0] - a1 * segment1[0][1]
    if slope == "v":
        # check vertical
        if segment1[0][1] <= segment2[0][1] & segment2[1][1] <= segment1[1][1]:
            intersection_i = -(a1 * segment2[0][1] + b1)
            if (intersection_i > min(segment1[0][0], segment2[0][0])) | (
                intersection_i < max(segment1[1][0], segment2[1][0])
            ):
                return False  # intersection is out of bound
            return (int(intersection_i), segment2[0][1])
        else:
            return False
    # check horizontal
    elif segment1[0][0] >= segment2[0][0] & segment2[1][0] >= segment1[1][0]:
        intersection_j = (-segment2[0][0] - b1) / a1
        if (intersection_j < max(segment1[0][1], segment2[0][1])) | (
            intersection_j > min(segment1[1][1], segment2[1][1])
        ):
            return False  # intersection is out of bound
        return (segment2[0][0], int(intersection_j))
    else:
        return False


def FindMinEnLocAlkmer(
    seq,
    seq_compl,
    seq_indxd,
    seq_compl_indxd,
    k,
    energy_threshold,
    handle_length_threshold,
    need_suboptimal,
    kmers_stacking_matrix,
):
    if (seq_indxd == False) | (seq_compl_indxd == False):
        return 0
    else:
        seq_indxd_tmp = seq_indxd[:]
        seq_compl_indxd_tmp = seq_compl_indxd[:]
        (
            D,
            B,
            S,
            min_energy,
            end_pos_i,
            end_pos_j,
            start_pos_i,
            start_pos_j,
            seq_indxd_tmp,
            seq_compl_indxd_tmp,
            seq_length,
            seq_compl_length,
            seq,
            seq_compl,
        ) = Initiate_with_kmers(
            seq, seq_compl, seq_indxd_tmp, seq_compl_indxd_tmp, kmers_stacking_matrix, k
        )
        only_optimal = True
        if (min_energy != 0) & (seq != False):
            only_optimal = False
            alignments = []
            seq_indxd_tmp.extend([kmers_stacking_matrix.shape[0] - 1] * (k + 2))
            seq_compl_indxd_tmp.extend([kmers_stacking_matrix.shape[0] - 1] * (k + 2))
            # go through matrices and fill them in
            for i in range(len(seq_compl) - k, k + 2, -1):
                I = seq_compl_length - k - i
                for j in range(k - 1, len(seq) - k - 3):
                    J = seq_length - j - 1
                    if D[i][j] != inf and D[i][j] != 0:  # found kmer stacking
                        S_head = S[i][j][: S[i][j].find(")")]
                        S_tail = S[i][j][S[i][j].find(")") :]
                        # stem
                        new_en = (
                            D[i][j]
                            + stacking_matrix[
                                Dic_bp.get(seq_compl[i - 1] + seq[j + 1], 6)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 6)]
                        )
                        argmin_ = argmin([0, D[i - 1][j + 1], new_en])
                        B[i - 1][j + 1] = Backtrack(argmin_, B[i - 1][j + 1], B[i][j])
                        D[i - 1][j + 1] = [0, D[i - 1][j + 1], new_en][argmin_]
                        S[i - 1][j + 1] = [
                            "*",
                            S[i - 1][j + 1],
                            S_head + "()" + S_tail,
                        ][argmin_]
                        # bulge01 (seq has 1 more nt)
                        new_en = (
                            D[i][j]
                            + bulge_list[1]
                            + stacking_matrix[
                                Dic_bp.get(seq_compl[i - 1] + seq[j + 2], 6)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 6)]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 1 - k], seq_compl_indxd_tmp[I + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - k][j + 1 + k], new_en])
                        B[i - k][j + 1 + k] = Backtrack(
                            argmin_, B[i - k][j + 1 + k], B[i][j]
                        )
                        D[i - k][j + 1 + k] = [0, D[i - k][j + 1 + k], new_en][argmin_]
                        S[i - k][j + 1 + k] = [
                            "*",
                            S[i - k][j + 1 + k],
                            S_head + "." + "(" * k + ")" * k + S_tail,
                        ][argmin_]
                        # bulge10 (seq_compl has 1 more nt)
                        new_en = (
                            D[i][j]
                            + bulge_list[1]
                            + stacking_matrix[
                                Dic_bp.get(seq_compl[i - 2] + seq[j + 1], 6)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 6)]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - k], seq_compl_indxd_tmp[I + 1 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 1 - k][j + k], new_en])
                        B[i - 1 - k][j + k] = Backtrack(
                            argmin_, B[i - 1 - k][j + k], B[i][j]
                        )
                        D[i - 1 - k][j + k] = [0, D[i - 1 - k][j + k], new_en][argmin_]
                        S[i - 1 - k][j + k] = [
                            "*",
                            S[i - 1 - k][j + k],
                            S_head + "(" * k + ")" * k + "." + S_tail,
                        ][argmin_]
                        # bulge 02
                        new_en = (
                            D[i][j]
                            + bulge_list[2]
                            + (
                                TerminalAU
                                if Dic_bp.get(seq_compl[i - 1] + seq[j + 3], 6) > 1
                                or Dic_bp.get(seq[j] + seq_compl[i], 6) > 1
                                else 0
                            )
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 2 - k], seq_compl_indxd_tmp[I + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - k][j + 2 + k], new_en])
                        B[i - k][j + 2 + k] = Backtrack(
                            argmin_, B[i - k][j + 2 + k], B[i][j]
                        )
                        D[i - k][j + 2 + k] = [0, D[i - k][j + 2 + k], new_en][argmin_]
                        S[i - k][j + 2 + k] = [
                            "*",
                            S[i - k][j + 2 + k],
                            S_head + ".." + "(" * k + ")" * k + S_tail,
                        ][argmin_]
                        # bulge 20
                        new_en = (
                            D[i][j]
                            + bulge_list[2]
                            + (
                                TerminalAU
                                if Dic_bp.get(seq_compl[i - 3] + seq[j + 1], 6) > 1
                                or Dic_bp.get(seq[j] + seq_compl[i], 6) > 1
                                else 0
                            )
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - k], seq_compl_indxd_tmp[I + 2 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 2 - k][j + k], new_en])
                        B[i - 2 - k][j + k] = Backtrack(
                            argmin_, B[i - 2 - k][j + k], B[i][j]
                        )
                        D[i - 2 - k][j + k] = [0, D[i - 2 - k][j + k], new_en][argmin_]
                        S[i - 2 - k][j + k] = [
                            "*",
                            S[i - 2 - k][j + k],
                            S_head + "(" * k + ")" * k + ".." + S_tail,
                        ][argmin_]
                        # loop11
                        new_en = (
                            D[i][j]
                            + intl11_matrix[
                                Dic_bp.get(seq_compl[i - 2] + seq[j + 2], 7)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 7)][
                                Dic_nt.get(seq_compl[i - 1], 5)
                            ][
                                Dic_nt.get(seq[j + 1], 5)
                            ]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 1 - k], seq_compl_indxd_tmp[I + 1 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 1 - k][j + 1 + k], new_en])
                        B[i - 1 - k][j + 1 + k] = Backtrack(
                            argmin_, B[i - 1 - k][j + 1 + k], B[i][j]
                        )
                        D[i - 1 - k][j + 1 + k] = [0, D[i - 1 - k][j + 1 + k], new_en][
                            argmin_
                        ]
                        S[i - 1 - k][j + 1 + k] = [
                            "*",
                            S[i - 1 - k][j + 1 + k],
                            S_head + "." + "(" * k + ")" * k + "." + S_tail,
                        ][argmin_]
                        # loop12
                        new_en = (
                            D[i][j]
                            + intl12_matrix[
                                Dic_bp.get(seq_compl[i - 2] + seq[j + 3], 7)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 7)][
                                Dic_nt.get(seq_compl[i - 1], 5)
                            ][
                                Dic_nt.get(seq[j + 1], 5)
                            ][
                                Dic_nt.get(seq[j + 2], 5)
                            ]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 2 - k], seq_compl_indxd_tmp[I + 1 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 1 - k][j + 2 + k], new_en])
                        B[i - 1 - k][j + 2 + k] = Backtrack(
                            argmin_, B[i - 1 - k][j + 2 + k], B[i][j]
                        )
                        D[i - 1 - k][j + 2 + k] = [0, D[i - 1 - k][j + 2 + k], new_en][
                            argmin_
                        ]
                        S[i - 1 - k][j + 2 + k] = [
                            "*",
                            S[i - 1 - k][j + 2 + k],
                            S_head + ".." + "(" * k + ")" * k + "." + S_tail,
                        ][argmin_]
                        # loop21
                        new_en = (
                            D[i][j]
                            + intl12_matrix[Dic_bp.get(seq[j] + seq_compl[i], 7)][
                                Dic_bp.get(seq_compl[i - 3] + seq[j + 2], 7)
                            ][Dic_nt.get(seq[j + 1], 5)][
                                Dic_nt.get(seq_compl[i - 2], 5)
                            ][
                                Dic_nt.get(seq_compl[i - 1], 5)
                            ]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 1 - k], seq_compl_indxd_tmp[I + 2 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 2 - k][j + 1 + k], new_en])
                        B[i - 2 - k][j + 1 + k] = Backtrack(
                            argmin_, B[i - 2 - k][j + 1 + k], B[i][j]
                        )
                        D[i - 2 - k][j + 1 + k] = [0, D[i - 2 - k][j + 1 + k], new_en][
                            argmin_
                        ]
                        S[i - 2 - k][j + 1 + k] = [
                            "*",
                            S[i - 2 - k][j + 1 + k],
                            S_head + "." + "(" * k + ")" * k + ".." + S_tail,
                        ][argmin_]
                        # loop22
                        new_en = (
                            D[i][j]
                            + intl22_matrix[
                                Dic_bp.get(seq_compl[i - 3] + seq[j + 3], 7)
                            ][Dic_bp.get(seq[j] + seq_compl[i], 7)][
                                Dic_nt.get(seq_compl[i - 2], 5)
                            ][
                                Dic_nt.get(seq_compl[i - 1], 5)
                            ][
                                Dic_nt.get(seq[j + 1], 5)
                            ][
                                Dic_nt.get(seq[j + 2], 5)
                            ]
                            + kmers_stacking_matrix[
                                seq_indxd_tmp[J - 2 - k], seq_compl_indxd_tmp[I + 2 + k]
                            ]
                        )
                        argmin_ = argmin([0, D[i - 2 - k][j + 2 + k], new_en])
                        B[i - 2 - k][j + 2 + k] = Backtrack(
                            argmin_, B[i - 2 - k][j + 2 + k], B[i][j]
                        )
                        D[i - 2 - k][j + 2 + k] = [0, D[i - 2 - k][j + 2 + k], new_en][
                            argmin_
                        ]
                        S[i - 2 - k][j + 2 + k] = [
                            "*",
                            S[i - 2 - k][j + 2 + k],
                            S_head + ".." + "(" * k + ")" * k + ".." + S_tail,
                        ][argmin_]
                        # check if found min Energy
                        argmin_ = argmin(
                            [
                                min_energy,
                                D[i - 1][j + 1],
                                D[i - k][j + 1 + k],
                                D[i - 1 - k][j + k],
                                D[i - k][j + 2 + k],
                                D[i - 2 - k][j + k],
                                D[i - 1 - k][j + 1 + k],
                                D[i - 1 - k][j + 2 + k],
                                D[i - 2 - k][j + 1 + k],
                                D[i - 2 - k][j + 2 + k],
                            ]
                        )
                        end_pos_i, end_pos_j = End_coords(
                            argmin_, i, j, end_pos_i, end_pos_j, k=k
                        )
                        start_pos_i, start_pos_j = Start_coords(
                            argmin_, B[i][j], start_pos_i, start_pos_j
                        )
                        min_energy = [
                            min_energy,
                            D[i - 1][j + 1],
                            D[i - k][j + 1 + k],
                            D[i - 1 - k][j + k],
                            D[i - k][j + 2 + k],
                            D[i - 2 - k][j + k],
                            D[i - 1 - k][j + 1 + k],
                            D[i - 1 - k][j + 2 + k],
                            D[i - 2 - k][j + 1 + k],
                            D[i - 2 - k][j + 2 + k],
                        ][argmin_]
            # save optimal structure
            if (
                (min_energy / 100.0 <= energy_threshold)
                & (end_pos_j - start_pos_j + 1 >= handle_length_threshold)
                & (start_pos_i - end_pos_i + 1 >= handle_length_threshold)
            ):
                alignments.append(
                    (
                        min_energy / 100,
                        start_pos_j,
                        end_pos_j,
                        end_pos_i - k - 2,
                        start_pos_i - k - 2,
                        seq[start_pos_j : end_pos_j + 1],
                        seq_compl[end_pos_i : start_pos_i + 1],
                        S[end_pos_i, end_pos_j],
                    )
                )
                if (only_optimal == False) & (need_suboptimal == True):
                    min_energy_potential = 0
                    while True:
                        # LOOK FOR SUBOPTIMAL STRUCTURES
                        ## put zeros in square from start to end of the best alignment
                        len_zero_i = start_pos_i - end_pos_i + 1
                        len_zero_j = end_pos_j - start_pos_j + 1
                        zero_matrix_D = full((len_zero_i, len_zero_j), 0)
                        D[
                            end_pos_i : start_pos_i + 1, start_pos_j : end_pos_j + 1
                        ] = zero_matrix_D
                        ## find next suboptimal structure
                        argmin_ = D.argmin()
                        end_pos_i_potential, end_pos_j_potential = unravel_index(
                            argmin_, D.shape
                        )
                        start_pos_i_potential, start_pos_j_potential = Start_coords(
                            1, B[end_pos_i_potential][end_pos_j_potential], 0, 0
                        )
                        min_energy_potential = D[end_pos_i_potential][
                            end_pos_j_potential
                        ]
                        ## check it's energy is higher than threshold
                        if min_energy_potential / 100.0 > energy_threshold:
                            break
                            ## check it doesn't overlap any previous alignments
                        overlap = False
                        for alignment in alignments:
                            start_pos_j = alignment[1]
                            end_pos_j = alignment[2] + k + 2
                            end_pos_i = alignment[3]
                            start_pos_i = alignment[4] + k + 2
                            if Check_ranges_overlap(
                                start_pos_j,
                                end_pos_j,
                                start_pos_j_potential,
                                end_pos_j_potential,
                            ) & Check_ranges_overlap(
                                end_pos_i_potential,
                                start_pos_i_potential,
                                end_pos_i,
                                start_pos_i,
                            ):
                                overlap = True
                                break
                                ## if it overlaps, delete it from D matrix
                        if overlap:
                            D[end_pos_i_potential][end_pos_j_potential] = 0
                        ## if it doesn't overlap, save the suboptimal structure
                        else:
                            min_energy = min_energy_potential
                            start_pos_i = start_pos_i_potential
                            end_pos_i = end_pos_i_potential
                            start_pos_j = start_pos_j_potential
                            end_pos_j = end_pos_j_potential
                            if (
                                end_pos_j - start_pos_j + 1 >= handle_length_threshold
                            ) & (
                                start_pos_i - end_pos_i + 1 >= handle_length_threshold
                            ):  # check the al is long enough
                                alignments.append(
                                    (
                                        min_energy / 100,
                                        start_pos_j,
                                        end_pos_j,
                                        end_pos_i - k - 2,
                                        start_pos_i - k - 2,
                                        seq[start_pos_j : end_pos_j + 1],
                                        seq_compl[end_pos_i : start_pos_i + 1],
                                        S[end_pos_i, end_pos_j],
                                    )
                                )
                return alignments
            else:
                return 0
        else:
            return 0


def process_pair(
    seq,
    seq_compl,
    k,
    energy_threshold,
    handle_length_threshold,
    need_suboptimal,
    kmers_stacking_matrix,
):
    seq_indxd = Index_seq(seq, k)
    seq_compl_indxd = Index_seq(seq_compl, k)
    res = FindMinEnLocAlkmer(
        seq,
        seq_compl,
        seq_indxd,
        seq_compl_indxd,
        k,
        energy_threshold,
        handle_length_threshold,
        need_suboptimal,
        kmers_stacking_matrix,
    )
    return res


def process_pair_full(args):
    res = process_pair(*args)
    if res != 0:
        res = [str(el) for el in res[0]]
        res_fmt = "\t".join(
            [
                res[7],
                res[1] + "," + res[2],
                ":",
                res[3] + "," + res[4],
                "(" + res[0] + ")",
            ]
        )
        return res_fmt
    else:
        return ""


def main(argv):

    parser = argparse.ArgumentParser(description="Process some integers.")
    parser.add_argument("-k", "--k", dest="k", type=int, default=5, help="kmer_length")
    parser.add_argument(
        "-a",
        "--handle_len_min",
        dest="handle_length_threshold",
        type=int,
        default=10,
        help="handle_len_min",
    )
    parser.add_argument(
        "-e",
        "--energy_max",
        dest="energy_threshold",
        type=float,
        default=-15,
        help="energy_max, kcal/mol",
    )
    parser.add_argument(
        "-u",
        "--need_subopt",
        dest="need_suboptimal",
        type=bool,
        default=True,
        help="need_subopt",
    )
    parser.add_argument(
        "-d",
        "--gt_threshold",
        dest="GT_threshold",
        type=int,
        default=2,
        help="gt_threshold",
    )
    parser.add_argument(
        "-j",
        dest="nproc",
        type=int,
        default=4,
        help="number of processes",
    )

    args = parser.parse_args()

    kmers_stacking_matrix = load(
        "../data/"
        + str(args.k)
        + str(args.GT_threshold)
        + "mers_stacking_energy_binary.npy"
    )

    pool = Pool(args.nproc)
    pairs = [line.strip().split("\t") for line in sys.stdin if line.strip()]
    content = [
        (
            pair[0],
            pair[1],
            args.k,
            args.energy_threshold,
            args.handle_length_threshold,
            args.need_suboptimal,
            kmers_stacking_matrix,
        )
        for pair in pairs
    ]

    for r in pool.map(process_pair_full, content):
        print(r)


if __name__ == "__main__":
    main(sys.argv[1:])