Experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "959bbb85",
   "metadata": {},
   "source": [
    "# Contents\n",
    "\n",
    "* [Read Filtering](#Read-Filtering)\n",
    "* [Result Evaluation](#Result-Evaluation)\n",
    "\n",
    "# Read Filtering\n",
    "\n",
    "The code in this section describes how the simulated set of reads is filtered to exclude all reads for which edlib finds more than 20 mapping positions.\n",
    "\n",
    "###### Requirements:\n",
    "* Set of simulated reads (`simulations/reads/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10.fasta`)\n",
    "* edlib mapping results for whole read set (`simulations/edlibMappings/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_ri0-69400.er`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e9ea26e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "\n",
    "edlibResultPath = \"simulations/edlibMappings/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0\" + \\\n",
    "\"009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_ri0-69400.er\"\n",
    "rarelyMappedReadIds = {}\n",
    "readId = -1\n",
    "resultCounter = 0\n",
    "\n",
    "for l in open(edlibResultPath, 'r'):\n",
    "    if l.startswith(\"simulations/edlibMappings/\"):\n",
    "        if readId > -1 and resultCounter <= 20:\n",
    "            rarelyMappedReadIds[readId] = None\n",
    "            \n",
    "        readId = int(l.split(\"_ri\")[1].split(\".er\")[0])\n",
    "        resultCounter = 0\n",
    "    else:\n",
    "        resultCounter += 1\n",
    "        \n",
    "if resultCounter <= 20:\n",
    "    rarelyMappedReadIds[readId] = None\n",
    "    \n",
    "allReadsPath = \"simulations/reads/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.000909090909\" + \\\n",
    "\"0909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10.fasta\"\n",
    "    \n",
    "rarelyMappedReadRecords = [r for r in SeqIO.parse(allReadsPath, \"fasta\") if int(r.id) in rarelyMappedReadIds]\n",
    "rarelyMappedReadsFilePath = \"simulations/reads/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0\" + \\\n",
    "\".0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20.fasta\"\n",
    "SeqIO.write(rarelyMappedReadRecords, open(rarelyMappedReadsFilePath, 'w'), \"fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad378f57",
   "metadata": {},
   "source": [
    "# Result Evaluation\n",
    "\n",
    "The code in this section allows to perform a mapping accuracy evaluation using BLAST confirmed mapping results as the ground truth.\n",
    "\n",
    "###### Requirements:\n",
    "* T2T assembly of human chromosome Y (`simulations/genomes/t2thumanChrY.fasta`)\n",
    "* Filtered set of simulated reads (`simulations/reads/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20.fasta`)\n",
    "* List of k-mers occurring more than 100 times (`highAbundKmersMiniK15w10Lrgr100BtStrnds.txt`)\n",
    "* edlib mapping results for whole read set (`simulations/edlibMappings/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_ri0-69400.er`)\n",
    "* eskemap mapping results for filtered read set (`simulations/homologies/homologies_t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20_k15_w10_c1_u1_de0.08758516_in-231.1585158515809.txt`)\n",
    "* minimap2 mapping results (`simulations/minimap2Res/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20_k15.paf.gz`)\n",
    "* Winnowmap2 mapping results (`simulations/Winnowmap2Res/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20_k15.paf.gz`)\n",
    "* BLAST results for mapping results (`simulations/blastRes/*_e0.01.tsv`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2263f490",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "from collections import deque\n",
    "import gzip\n",
    "\n",
    "NT_IN_BITS = {'A': 0, 'C': 1, 'G': 2, 'T': 3}\n",
    "\n",
    "#Read FASTA file\n",
    "def readFasta(path):\n",
    "    return [r for r in SeqIO.parse(open(path, 'r'), \"fasta\")]\n",
    "\n",
    "#Get a subset of all edlib results that is based on a maximum number of results per read\n",
    "def createEdlibResSubsetBasedOnMaxResNbPerRead(allEdlibRes, maxResNb):\n",
    "    subSet = {}\n",
    "    \n",
    "    for r in allEdlibRes:\n",
    "        if len(allEdlibRes[r]) <= maxResNb:\n",
    "            subSet[r] = list(allEdlibRes[r])\n",
    "            \n",
    "    return subSet\n",
    "\n",
    "#This function filters edlib's results based on a distance threshold\n",
    "#Parameters:\n",
    "#    res: A dictionary of edlib results with integer-based read ids as keys and a list representing ...\n",
    "#    thres: Relative distance threshold\n",
    "def filterByDistanceThreshold(res, thres):\n",
    "    filteredRes = {}\n",
    "    \n",
    "    for r in res:\n",
    "        filteredRes[r] = []\n",
    "        absThres = thres * len(rdSeqs[f\"s_{r}\"].seq)\n",
    "        \n",
    "        for rs in res[r]:\n",
    "            if int(rs[-1]) <= absThres:\n",
    "                filteredRes[r].append(rs)\n",
    "        \n",
    "    return filteredRes\n",
    "\n",
    "#This function calculates the hash of a bitwise k-mer representation. The function is influenced by the code of \"The\n",
    "#minimizer Jaccard estimator is biased and inconsistent.\" from Belbasi et al. (function \"minimap2_hash(seed,v,mask)\"\n",
    "#in file \"minimap2_hash_uncompiled.py\").\n",
    "def getHash(kmer, mask):\n",
    "    u = kmer & mask\n",
    "    u = ((~u) + (u << 21)) & mask # u = (u<<21)-(u+1) = 77594587*u-1\n",
    "    u = u ^ (u >> 24)\n",
    "    u = ((u + (u << 3)) + (u << 8)) & mask # u *= 265\n",
    "    u = u ^ (u >> 14)\n",
    "    u = ((u + (u << 2)) + (u << 4)) & mask # u *= 21\n",
    "    u = u ^ (u >> 28)\n",
    "    u = (u + (u << 31)) & mask # u *= 2147483649\n",
    "\n",
    "    return u\n",
    "\n",
    "#This function calculates the minimizer sketch of a sequence. It is influenced by the code of \"The minimizer Jaccard\n",
    "#estimator is biased and inconsistent.\" from Belbasi et al. (function \"winnowed_minimizers_linear(perm,windowSize)\" \n",
    "#in file \"winnowed_minimizers.py\").\n",
    "def calcMiniSketch(seq, k, w):\n",
    "    sketch = []\n",
    "    #A deque to store k-mer hashes inside the current window\n",
    "    windowKmers = deque()\n",
    "    mask = (4 ** k) - 1\n",
    "    lastIdx = -1\n",
    "\n",
    "    for i in range(len(seq) - k + 1):\n",
    "        kmerBits = 0\n",
    "        kmerBitsRevComp = 0\n",
    "        windowBorder = i - (w - 1)\n",
    "        \n",
    "        #Get bit representation of k-mer\n",
    "        for c in seq[i:i+k]:\n",
    "            kmerBits = (kmerBits << 2) + NT_IN_BITS[c]\n",
    "\n",
    "        #Get bit representation of k-mer's reverse complement\n",
    "        for c in str(Seq(seq[i:i+k]).reverse_complement()):\n",
    "            kmerBitsRevComp = (kmerBitsRevComp << 2) + NT_IN_BITS[c]\n",
    "\n",
    "        #If a k-mer is its own reverse complement we skip it\n",
    "        if kmerBits == kmerBitsRevComp:\n",
    "            continue\n",
    "\n",
    "        #Depending on which hash is smaller we consider either a k-mer or its reverse complement per position\n",
    "        if kmerBits < kmerBitsRevComp:\n",
    "            kmer = (i, kmerBits, getHash(kmerBits, mask))\n",
    "        else:\n",
    "            #A k-mer is a pair of k-mer's start position and its hash\n",
    "            kmer = (i, kmerBitsRevComp, getHash(kmerBitsRevComp, mask))\n",
    "            \n",
    "        #Remove all k-mers with a hash value larger than the newly calculated one\n",
    "        while (len(windowKmers) > 0) and (windowKmers[-1][2] > kmer[2]):\n",
    "            windowKmers.pop()\n",
    "\n",
    "        #Save new k-mer as window k-mer\n",
    "        windowKmers.append(kmer)\n",
    "\n",
    "        #Remove k-mer if it is not any longer inside the window\n",
    "        while (len(windowKmers) > 0) and (windowKmers[0][0] < windowBorder):\n",
    "            windowKmers.popleft()\n",
    "\n",
    "        #As soon as we have seen a first full window of k-mers choose a minimizer\n",
    "        if (windowBorder >= 0) and (len(windowKmers) > 0):      \n",
    "            #We do not choose the same minimizer for a second time\n",
    "            if lastIdx != windowKmers[0][0]:\n",
    "                sketch.append((windowKmers[0][0]+k-1, windowKmers[0][1], windowKmers[0][2]))\n",
    "                lastIdx = windowKmers[0][0]\n",
    "                \n",
    "            while len(windowKmers) > 1 and windowKmers[0][1] == windowKmers[1][1]:\n",
    "                windowKmers.popleft()\n",
    "                sketch.append((windowKmers[0][0]+k-1, windowKmers[0][1], windowKmers[0][2]))    \n",
    "                lastIdx = windowKmers[0][0] \n",
    "\n",
    "    #If our sequence was too small to get a full window of k-mers to consider take the smallest one found so far\n",
    "    if windowBorder < 0 and len(windowKmers) > 0:\n",
    "        sketch.append((windowKmers[0][0]+k-1, windowKmers[0][1], windowKmers[0][2]))\n",
    "        \n",
    "        while len(windowKmers) > 1 and windowKmers[0][1] == windowKmers[1][1]:\n",
    "            windowKmers.popleft()\n",
    "            sketch.append((windowKmers[0][0]+k-1, windowKmers[0][1], windowKmers[0][2]))\n",
    "\n",
    "    return sketch\n",
    "\n",
    "def loadFindThomsRes(pth, refSketch):\n",
    "    compResPerRead = {}\n",
    "    \n",
    "    for l in open(pth, 'r'):\n",
    "        l = l.strip()\n",
    "\n",
    "        if l.startswith('s'):\n",
    "            lastRdId = int(l.split('_')[1])\n",
    "            compResPerRead[lastRdId] = []\n",
    "        else:\n",
    "            cols = l.split(' ')\n",
    "            s = refSketch[int(cols[1])][0] + 1 - K\n",
    "            e = refSketch[int(cols[3])][0]\n",
    "            compResPerRead[lastRdId].append({\"id\": lastRdId, \"start\": s, \"end\": e, \"score\": float(cols[5])})\n",
    "        \n",
    "    return compResPerRead\n",
    "\n",
    "#This function returns a dictionary with integer-based read ids as keys and a list of start stop tupels of homology \n",
    "#intervals reported by FindThoms as values. Stop coordinates may be corrected in the way that the k-1 overlap of the\n",
    "#last k-mer part of the alpha-homology is not included anymore.\n",
    "#Parameter: \n",
    "#res: Result of FindThoms as returned from function loadFindThomsRes\n",
    "#corrEnd: Flag indicating if the end coordinate shall be adaptated\n",
    "def parseCoordsFromFindThomsRes(res, corrEnd):\n",
    "    ints = {}\n",
    "        \n",
    "    for r in res:\n",
    "        if corrEnd:\n",
    "            #Start position was corrected inside loadFindThomsRes already\n",
    "            ints[r] = [(h[\"start\"], h[\"end\"] + 1 - 15) for h in res[r]]\n",
    "        else:\n",
    "            ints[r] = [(h[\"start\"], h[\"end\"]) for h in res[r]]\n",
    "        \n",
    "    return ints\n",
    "\n",
    "#This function filters FindThom's results for intervals above a threshold score relative to read lengths and \n",
    "#interpolated using given decent and intercept\n",
    "#Parameters:\n",
    "#    res: A dictionary of results as returned by loadFindThomsRes(pth, refSketch)\n",
    "#    dec: Float describing the decent of the line used for threshold interpolation\n",
    "#    inter: Float describing the intercept of the line used for threshold interpolation\n",
    "def filterByScoreThreshold(res, dec, inter):\n",
    "    filteredRes = {}\n",
    "    \n",
    "    for r in res:\n",
    "        filteredRes[r] = []\n",
    "        absThres = dec * len(rdSeqs[f\"s_{r}\"].seq) + inter\n",
    "        \n",
    "        for rs in res[r]:\n",
    "            if rs[\"score\"] >= absThres:\n",
    "                filteredRes[r].append(rs)\n",
    "                \n",
    "    return filteredRes\n",
    "\n",
    "def parseCompressedPAF(filepath):\n",
    "    resPerRead = {}\n",
    "    \n",
    "    for l in gzip.open(filepath, 'rt'):\n",
    "        elems = l.strip().split('\\t')\n",
    "        rdid = int(elems[0].split('_')[1])\n",
    "        \n",
    "        if rdid in resPerRead:\n",
    "            resPerRead[rdid].append((int(elems[7]), int(elems[8]) - 1))\n",
    "        else:\n",
    "            resPerRead[rdid] = [(int(elems[7]), int(elems[8]) - 1)]\n",
    "            \n",
    "    return resPerRead\n",
    "\n",
    "#This function reads BLAST results given a list of BLAST result files and adds it to a dictionary\n",
    "#Parameters:\n",
    "#    resFilePaths: List of result file paths\n",
    "#    blastRes: A dictionary assumed to be shaped as described above\n",
    "def loadBLASTresults(resFilePaths, blastRes):\n",
    "    for f in resFilePaths:\n",
    "        for l in open(f, 'r'):\n",
    "            #Header line\n",
    "            if not l.startswith(\"s_\"):\n",
    "                #Find out read id\n",
    "                idPattern = search(\"s_[0-9]+\", l)\n",
    "                #Find out reference substring range\n",
    "                rangePattern = search(\"ref[0-9]+-[0-9]+\", l)\n",
    "                #Parse start and end of reference substring\n",
    "                subStart, subEnd = [int(c) for c in rangePattern.group(0).split(\"ref\")[1].split('-')]\n",
    "                #Parse read id\n",
    "                readId = int(idPattern.group(0).split(\"s_\")[1])\n",
    "\n",
    "                #Check if this read has been seen before and add it if necessary\n",
    "                if not readId in blastRes:\n",
    "                    blastRes[readId] = {}\n",
    "                    \n",
    "                #Check if even this substring has been seen before for this read\n",
    "                if not (subStart, subEnd) in blastRes[readId]:\n",
    "                    blastRes[readId][(subStart, subEnd)] = []\n",
    "            else:\n",
    "                elems = l.split('\\t')\n",
    "                \n",
    "                #Sanity check\n",
    "                #This won't work for FindThoms results since I messed up FASTA headers there...\n",
    "                if len(elems[3].split('-')) < 2:\n",
    "                    hitRefStart = subStart\n",
    "                    hitRefEnd = subEnd\n",
    "                else:\n",
    "                    #Check if reference substring in the header matches that of the current hit\n",
    "                    hitRefStart, hitRefEnd = [int(c) for c in elems[3].split(\"ref\")[1].split('-')]\n",
    "                \n",
    "                if hitRefStart != subStart or hitRefEnd != subEnd:\n",
    "                    print(f\"Some result for read s_{readId} in file {f} has different substring coordinates as \" + \\\n",
    "                          \"expected from the header\")\n",
    "                    print(\"Hit coordinates:\", hitRefStart, hitRefEnd)\n",
    "                    print(\"Coordinates in header:\", subStart, subEnd)\n",
    "                    \n",
    "#                     p0 = f\"../simulations/mappedAreas1/otherTools/sub_s_{readId}_ref{subStart}-{subEnd}.fasta\"\n",
    "#                     p1 = f\"../simulations/mappedAreas1/otherTools/sub_s_{readId}_ref{hitRefStart}-{hitRefEnd}\" + \\\n",
    "#     \".fasta\"\n",
    "        \n",
    "#                     if not exists(p0):\n",
    "#                         SeqIO.write(SeqRecord(Seq(refSeq[subStart:subEnd+1]), id=f\"ref{subStart}-{subEnd}\"), \\\n",
    "#                          0           open(p0, 'w'), \"fasta\")\n",
    "                    \n",
    "#                     if not exists(p1):\n",
    "#                         SeqIO.write(SeqRecord(Seq(refSeq[hitRefStart:hitRefEnd+1]), id=f\"ref{hitRefStart}-\" + \\\n",
    "#                                               f\"{hitRefEnd}\"), open(p1, 'w'), \"fasta\")\n",
    "\n",
    "                    continue\n",
    "                \n",
    "                #Filter for evalue just to be sure\n",
    "                if float(elems[6]) <= 0.01:\n",
    "                    blastRes[readId][(subStart, subEnd)].append({\"rdStart\": int(elems[1]), \"rdEnd\": int(elems[2]), \\\n",
    "                                                                 \"refStart\": int(elems[4]), \"refEnd\": int(elems[5])\\\n",
    "                                                                 , \"strand\": elems[13], \"qcov\": int(elems[14]), \\\n",
    "                                                                 \"eval\": float(elems[6])})\n",
    "                    \n",
    "#Generate BLAST result subsets according to e-value thresholds\n",
    "def filterByEval(res, ethr):\n",
    "    filteredRes = {}\n",
    "    \n",
    "    for r in res:\n",
    "        filteredRes[r] = {}\n",
    "        \n",
    "        for s in res[r]:\n",
    "            for f in res[r][s]:\n",
    "                if f[\"eval\"] <= ethr:\n",
    "                    if s in filteredRes[r]:\n",
    "                        filteredRes[r][s].append(f)\n",
    "                    else:\n",
    "                        filteredRes[r][s] = [f]\n",
    "                        \n",
    "    return filteredRes\n",
    "\n",
    "#Filter BLAST results by e-value\n",
    "def filterByEvalRng(res, evals):\n",
    "    fltRes = {}\n",
    "    \n",
    "    for t in evals:\n",
    "        fltRes[t] = filterByEval(res, t)\n",
    "        \n",
    "    return fltRes\n",
    "\n",
    "#Form BLAST result subsets for tool result subsets we want to evaluate\n",
    "def filterByRes(totRes, filtRes):\n",
    "    blastSubsetRes = {}\n",
    "    \n",
    "    #For each tool\n",
    "    for t in filtRes:\n",
    "        for r in filtRes[t]:\n",
    "            for f in filtRes[t][r]:\n",
    "                if f in totRes[r]:\n",
    "                    if r in blastSubsetRes:\n",
    "                        if not f in blastSubsetRes[r]:\n",
    "                            blastSubsetRes[r][f] = totRes[r][f]\n",
    "                    else:\n",
    "                        blastSubsetRes[r] = {f: totRes[r][f]}\n",
    "            \n",
    "    return blastSubsetRes\n",
    "\n",
    "#Generate BLAST result subsets based on tool thresholds\n",
    "def genBlResSubsFrTlRes(tlThrBlRsDict, res2Flt, tlResCrds2Flt, tlNm):\n",
    "    thresholds = sorted(tlThrBlRsDict.keys(), reverse=True)\n",
    "    \n",
    "    for t in thresholds[1:]:\n",
    "        res2Flt[tlNm] = tlResCrds2Flt[t]\n",
    "        tlThrBlRsDict[t] = filterByRes(tlThrBlRsDict[thresholds[0]], res2Flt)\n",
    "        \n",
    "#Generate ground truth data set for result subsets\n",
    "def detCorrMappings(bRes):\n",
    "    corrMappings = {}\n",
    "\n",
    "    for r in bRes:\n",
    "        corrMappings[r] = []\n",
    "        rdLen = len(rdSeqs[f\"s_{r}\"].seq)\n",
    "\n",
    "        for i in bRes[r]:\n",
    "            subStrLen = i[1] - i[0] + 1\n",
    "            partRes = {\"minus\": [], \"plus\": []}\n",
    "            noRes = True\n",
    "\n",
    "            for br in bRes[r][i]:\n",
    "                covSubStr = abs(br[\"refEnd\"] - br[\"refStart\"] + 1)\n",
    "            \n",
    "                if covSubStr >= 0.9 * subStrLen or br[\"qcov\"] >= 90:\n",
    "                    corrMappings[r].append(i)\n",
    "                    noRes = False\n",
    "                    break\n",
    "                else:\n",
    "                    #TODO: This can be done better\n",
    "                    if br[\"strand\"] == \"plus\":\n",
    "                        start = br[\"refStart\"]\n",
    "                        end = br[\"refEnd\"]\n",
    "                    else:\n",
    "                        start = br[\"refEnd\"]\n",
    "                        end = br[\"refStart\"]\n",
    "\n",
    "                    ovl = False\n",
    "\n",
    "                    for intv in partRes[br[\"strand\"]]:\n",
    "                        if (start <= intv[0] and end >= intv[0]) or (start <= i[1] and end >= i[1]):\n",
    "                            ovl = True\n",
    "                            break\n",
    "\n",
    "                    if not ovl:\n",
    "                        partRes[br[\"strand\"]].append((start, end))\n",
    "                    \n",
    "            if noRes:\n",
    "                lenResParts = [0, 0]\n",
    "\n",
    "                for j, s in enumerate(partRes):\n",
    "                    for ints in partRes[s]:\n",
    "                        lenResParts[j] += ints[1] - ints[0] + 1\n",
    "\n",
    "                if max(lenResParts) >= 0.9 * subStrLen or max(lenResParts) >= 0.9 * rdLen:\n",
    "                    corrMappings[r].append(i)\n",
    "        \n",
    "    return corrMappings\n",
    "\n",
    "#This function returns a dictionary of k-mers from the given reference sketch that are covered by the given results \n",
    "#of some tool per each read from the given set of read ids.\n",
    "#refSk: A list of start positions of k-mers in a sequence that are chosen as being part of the sketch\n",
    "#res: A dictionary with read ids as keys and lists of start and stop coordinate tuples as values\n",
    "def getCovKmers(refSk, res):\n",
    "    covKmers = {}\n",
    "    #Sort start positions in sketch\n",
    "    sortedRefSk = sorted(refSk)\n",
    "    \n",
    "    #Iterate over reads:\n",
    "    for r in res:\n",
    "        #Initialize dictionary for this read\n",
    "        covKmers[r] = {}\n",
    "        \n",
    "        #If there are no results we are done for this read\n",
    "        if len(res[r]) == 0:\n",
    "            continue\n",
    "            \n",
    "        #Initialize a flag marking that we have reached the end of our result list\n",
    "        endReached = False\n",
    "        #Sort result coordinates for this read increasingly by start position\n",
    "        coords = sorted(res[r], key=lambda e: e[0])\n",
    "        #Initialize an index to know which result positions we consider currently (we start with the one that has \n",
    "        #the smallest start coordinate)\n",
    "        i = 0\n",
    "        \n",
    "        #Iterate over sorted list of start positions\n",
    "        for p in sortedRefSk:\n",
    "            #We need to switch to a result which has a larger start position as long as the position of the refer-\n",
    "            #ence sketch k-mer is larger than the end position of our current result\n",
    "            while p > coords[i][1]:\n",
    "                #Move to the next result\n",
    "                i += 1\n",
    "                \n",
    "                #Check if there is another result at all\n",
    "                if i >= len(coords):\n",
    "                    endReached = True\n",
    "                    break\n",
    "                    \n",
    "            #Check if we set the flag\n",
    "            if endReached:\n",
    "                break\n",
    "                \n",
    "            #Check if the k-mer starts before the current result\n",
    "            if p < coords[i][0]:\n",
    "                continue\n",
    "                \n",
    "            #If we have reached this point the k-mer must fall inside the result\n",
    "            covKmers[r][p] = None\n",
    "        \n",
    "    return covKmers\n",
    "\n",
    "#This function creates a dictionary with integer-based read ids as keys and lists of alignment coordinates as values\n",
    "#from a given dictionary of edlib alignment results\n",
    "#Parameter: res: Dictionary of edlib alignment results with integer-based read ids as keys and lists of alignment\n",
    "#                infos parsed from file\n",
    "def getEdlibResAlnCrds(res):\n",
    "    return {r: [(int(rs[0]), int(rs[1])) for rs in res[r]] for r in res}\n",
    "\n",
    "#This function calculates the accurracy of a tool for a given reference sketch and result files of edlib and \n",
    "#the tool using the \"ref sketch method\"\n",
    "def calcFindThomsAccRefSkMeth(refSk, edlibRes, compRes):\n",
    "    tpFpTnFnDict = {\"TPs\": 0, \"FPs\": 0, \"TNs\": 0, \"FNs\": 0}\n",
    "    #Get number of all k-mers in reference sketch\n",
    "    nbAllKmers = len(refSk)\n",
    "    \n",
    "    #Testing\n",
    "    tps = 0\n",
    "    \n",
    "    #Iterate over read ids\n",
    "    for r in edlibRes:\n",
    "        #Iterate over k-mer positions\n",
    "        for p in edlibRes[r]:\n",
    "            #We have a TP if the k-mer is covered by both results\n",
    "            if r in compRes and p in compRes[r]:\n",
    "                tpFpTnFnDict[\"TPs\"] += 1\n",
    "            else:\n",
    "                tpFpTnFnDict[\"FNs\"] += 1\n",
    "                \n",
    "        if r in compRes:\n",
    "            #Iterate over k-mer positions\n",
    "            for p in compRes[r]:\n",
    "                #Count TPs again as a sanity check\n",
    "                if p in edlibRes[r]:\n",
    "                    tps += 1\n",
    "                else:\n",
    "                    tpFpTnFnDict[\"FPs\"] += 1\n",
    "                \n",
    "    #Testing\n",
    "    if not tps == tpFpTnFnDict[\"TPs\"]:\n",
    "        print(\"Sanity check failed\\nTP counts are different\")\n",
    "        \n",
    "    tpFpTnFnDict[\"TNs\"] = (len(edlibRes) * nbAllKmers) - sum([tpFpTnFnDict[c] for c in tpFpTnFnDict])\n",
    "        \n",
    "    print(tpFpTnFnDict)\n",
    "    \n",
    "    return tpFpTnFnDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8085f9bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load reads\n",
    "rarelyMappedReadsFilePath = \"simulations/reads/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0\" + \\\n",
    "\".0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20.fasta\"\n",
    "rdSeqs = readFasta(rarelyMappedReadsFilePath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7fc5079",
   "metadata": {},
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import colors\n",
    "import numpy as np\n",
    "\n",
    "#Load edlib results\n",
    "edlibResultPath = \"simulations/edlibMappings/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0\" + \\\n",
    "\"009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_ri0-69400.er\"\n",
    "edlibRes = {}\n",
    "\n",
    "for l in open(edlibResultPath, 'r'):\n",
    "    if l.startswith(\"simulations/edlibMappings/\"):\n",
    "        readId = int(l.split(\"_ri\")[1].split(\".er\")[0])\n",
    "        edlibRes[readId] = []\n",
    "    else:\n",
    "        edlibRes[readId].append(l.strip().split(' '))\n",
    "        \n",
    "#Create \"rarely mapped 20\" subset of edlib results\n",
    "edlibResMaxRes = {}\n",
    "edlibResMaxRes[20] = createEdlibResSubsetBasedOnMaxResNbPerRead(edlibRes, 20)\n",
    "\n",
    "#Filter results based on distance thresholds\n",
    "edlibThres = [0.01, 0.02, 0.03]\n",
    "filteredEdlibRes = {0.03: edlibResMaxRes[20]}\n",
    "\n",
    "for t in edlibThres[:2]:\n",
    "    filteredEdlibRes[t] = filterByDistanceThreshold(edlibResMaxRes[20], t)\n",
    "    \n",
    "#Load blacklisted k-mers\n",
    "K = 15\n",
    "blKmers = {int(l): None for l in open(\"highAbundKmersMiniK15w10Lrgr100BtStrnds.txt\", 'r')}\n",
    "#Load reference sequence\n",
    "refSeq = readFasta(\"simulations/genomes/t2thumanChrY.fasta\")[0].seq\n",
    "#Calculate minimizer sketch of reference without blacklisted k-mers\n",
    "fltRfSkMiK15W10 = [k for k in calcMiniSketch(str(refSeq), K, 10) if k[2] not in blKmers]\n",
    "fltRfSkMiKmrStPosK15W10 = [k[0] + 1 - K for k in fltRfSkMiK15W10 if k[2] not in blKmers]\n",
    "#Load eskemap results\n",
    "p = \"simulations/homologies/homologies_t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i0.0009090\" + \\\n",
    "\"909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_rm20_k15_w10_c1_u1_de0.08758516_in-2\" + \\\n",
    "\"31.1585158515809.txt\"\n",
    "FindThomsNestResApprxMppngRr20 = loadFindThomsRes(p, fltRfSkMiK15W10)\n",
    "#Generate FindThoms result subsets according to thresholds\n",
    "FindThomsDecs = {0.9: 0.08790719, 0.8: 0.08817002, 0.7: 0.08841644}\n",
    "FindThomsInters = {0.9: -193.79071907191974, 0.8: -110.81700170017575, 0.7: -60.04164416442654}\n",
    "dpFltHomCrdsRr20 = parseCoordsFromFindThomsRes(FindThomsNestResApprxMppngRr20, False)\n",
    "dpFltCrds = {0.95: dpFltHomCrdsRr20}\n",
    "\n",
    "for i in FindThomsDecs:\n",
    "    fltHoms = filterByScoreThreshold(FindThomsNestResApprxMppngRr20, FindThomsDecs[i], FindThomsInters[i])\n",
    "    dpFltCrds[i] = parseCoordsFromFindThomsRes(fltHoms, False)\n",
    "    \n",
    "#Load competitors' results\n",
    "pafRes = {}\n",
    "relPafRes = {}\n",
    "\n",
    "for t in [\"minimap2\", \"Winnowmap2\"]:\n",
    "    p = f\"simulations/{t}Res/t2thumanChrY_sr0.00010909090909090909_dr0.0009818181818181818_i\" + \\\n",
    "    \"0.0009090909090909091_sd7361077429744071834_lmn100_lmx1000000_lavg9000_ls7000_dp10_k15.paf.gz\"\n",
    "    pafRes[t] = parseCompressedPAF(p)\n",
    "    #Only consider results of relevant reads to speed this up\n",
    "    relPafRes[t] = {r: pafRes[t][r] for r in pafRes[t] if r in edlibResMaxRes[20]}\n",
    "    \n",
    "#Load BLAST results\n",
    "blastResReload = {}\n",
    "loadBLASTresults(glob(\"simulations/blastRes/*_e0.01.tsv\"), blastResReload)\n",
    "#Generate BLAST result subsets\n",
    "evalthrs = [0.01, 0.005, 0.001]\n",
    "filteredBlastRes = filterByEvalRng(blastResReload, evalthrs)\n",
    "\n",
    "#Do I need this?\n",
    "# edlibThrsAllHomBlstRs = {0.03: blastResReload, 0.02: None, 0.01: None}\n",
    "# allHomRes2Flt =  {\"ESKEMAP\": dpFltHomCrdsRr20, \"minimap2\": relPafRes[\"minimap2\"], \"Winnowmap2\": relPafRes\\\n",
    "#                   [\"Winnowmap2\"]}\n",
    "# fltEdlibResCrds = {}\n",
    "\n",
    "# for t in edlibThres:\n",
    "#     fltEdlibResCrds[t] = getEdlibResAlnCrds(filteredEdlibRes[t])\n",
    "    \n",
    "# genBlResSubsFrTlRes(edlibThrsAllHomBlstRs, allHomRes2Flt, fltEdlibResCrds, \"edlib\")\n",
    "# scrThrsAllHomBlstRs = {0.95: allHomBLASTres, 0.9: None, 0.8: None, 0.7: None}\n",
    "# genBlResSubsFrTlRes(scrThrsAllHomBlstRs, allHomRes2Flt, dpFltCrds, \"ESKEMAP\")\n",
    "\n",
    "#Generate ground truth sets\n",
    "corrResVarEval = {}\n",
    "\n",
    "for t in evalthrs:\n",
    "    corrResVarEval[t] = detCorrMappings(filteredBlastRes[t])\n",
    "    \n",
    "# corrResAllHomVarDstThr = {}\n",
    "\n",
    "#Do I need this?\n",
    "# for t in edlibThres:\n",
    "#     corrResAllHomVarDstThr[t] = detCorrMappings(edlibThrsAllHomBlstRs[t])\n",
    "    \n",
    "# corrResAllHomVarScrThr = {}\n",
    " \n",
    "# for t in [0.95, 0.9, 0.8, 0.7]:\n",
    "#     corrResAllHomVarScrThr[t] = detCorrMappings(scrThrsAllHomBlstRs[t])\n",
    "\n",
    "#Calculate covered k-mers\n",
    "covKmersCorrAllHomVarEval = {}\n",
    "\n",
    "for t in evalthrs:\n",
    "    covKmersCorrAllHomVarEval[t] = getCovKmers(fltRfSkMiKmrStPosK15W10, corrResVarEval[t])\n",
    "    \n",
    "covKmersEskemap = {}\n",
    "\n",
    "for t in FindThomsDecs:\n",
    "    covKmersEskemap[t] = getCovKmers(fltRfSkMiKmrStPosK15W10, dpFltCrds[i])\n",
    "    \n",
    "for t in edlibThres:\n",
    "    covKmersEdlib[t] = getCovKmers(fltRfSkMiKmrStPosK15W10, getEdlibResAlnCrds(filteredEdlibRes[t]))\n",
    "\n",
    "covKmersCompetitors = {}\n",
    "\n",
    "for t in [\"minimap2\", \"Winnowmap2\"]:\n",
    "    covKmersCompetitors[t] = getCovKmers(fltRfSkMiKmrStPosK15W10, relPafRes[t])\n",
    "\n",
    "#Calculate mapping accuracies\n",
    "tpfptnfnVarEval = {}\n",
    "\n",
    "for t in evalthrs:\n",
    "    tpfptnfnVarEval[t] = {}\n",
    "    tpfptnfnVarEval[t][\"edlib\"] = calcFindThomsAccRefSkMeth(fltRfSkMiK15W10, corrResVarEval[t], covKmersEdlib[0.03])\n",
    "    tpfptnfnVarEval[t][\"eskemap\"] = calcFindThomsAccRefSkMeth(fltRfSkMiK15W10, corrResVarEval[t], \\\n",
    "                                                              covKmersEskemap[0.95])\n",
    "\n",
    "    for c in [\"minimap2\", \"Winnowmap2\"]:\n",
    "        tpfptnfnVarEval[t][c] = calcFindThomsAccRefSkMeth(fltRfSkMiK15W10, corrResVarEval[t], \\\n",
    "                                                          covKmersCompetitors[c])\n",
    "        \n",
    "tpfptnfnVarScrThr = {}\n",
    "\n",
    "for t in FindThomsDecs:\n",
    "    tpfptnfnVarScrThr[t] = calcFindThomsAccRefSkMeth(fltRfSkMiK15W10, corrResVarEval[0.01], covKmersEskemap[t])\n",
    "    \n",
    "tpfptnfnVarDstThr = {}\n",
    "\n",
    "for t in edlibThres:\n",
    "    tpfptnfnVarDstThr[t] = calcFindThomsAccRefSkMeth(fltRfSkMiK15W10, corrResVarEval[0.01], covKmersEdlib[t])\n",
    "        \n",
    "#Make a plot for varying e-values\n",
    "dotshapes = {\"ESKEMAP\": \"o\", \"minimap2\": \"^\", \"Winnowmap2\": \"v\", \"edlib\": \"P\"}\n",
    "dotColors = ['b', 'r', 'g']\n",
    "cnt = 0\n",
    "plt.figure()\n",
    "\n",
    "for et in evalthrs:\n",
    "    for t in tpfptnfnVarEval[et]:    \n",
    "        prec = tpfptnfnVarEval[et][t][\"TPs\"] / (tpfptnfnVarEval[et][t][\"TPs\"] + tpfptnfnVarEval[et][t][\"FPs\"])\n",
    "        recall = tpfptnfnVarEval[et][t][\"TPs\"] / (tpfptnfnVarEval[et][t][\"TPs\"] + tpfptnfnVarEval[et][t][\"FNs\"])\n",
    "        \n",
    "        print(\"E-value:\", et, \"Tool:\", t)\n",
    "        print(\"Precision:\", prec, \"Recall:\", recall)\n",
    "        \n",
    "        if cnt == 0:\n",
    "            plt.plot(prec, recall, dotshapes[t] + 'k', label=t)\n",
    "            \n",
    "        plt.plot(prec, recall, dotshapes[t] + dotColors[cnt])\n",
    "    \n",
    "    cnt += 1\n",
    "    \n",
    "plt.xlabel(\"Precision\")\n",
    "plt.ylabel(\"Recall\")\n",
    "plt.legend()\n",
    "plt.savefig(\"precRecVarEval.pdf\", format=\"pdf\")\n",
    "plt.show()\n",
    "#Make a plot for varying thresholds\n",
    "fig, ax = plt.subplots()\n",
    "scrThrs = sorted([0.95, 0.9, 0.8, 0.7])\n",
    "cmapDsts = plt.get_cmap('winter', 3)\n",
    "cmapScrs = plt.get_cmap('cool', len(scrThrs))\n",
    "\n",
    "for t in tpfptnfnVarEval[0.01]:        \n",
    "    prec = tpfptnfnVarEval[0.01][t][\"TPs\"] / (tpfptnfnVarEval[0.01][t][\"TPs\"] + tpfptnfnVarEval[0.01][t][\"FPs\"])\n",
    "    recall = tpfptnfnVarEval[0.01][t][\"TPs\"] / (tpfptnfnVarEval[0.01][t][\"TPs\"] + tpfptnfnVarEval[0.01][t][\"FNs\"])\n",
    "    ax.scatter(prec, recall, marker=dotshapes[t], c='k', label=t)\n",
    "    \n",
    "precs = []\n",
    "recalls = []\n",
    "shift = (max(scrThrs) - min(scrThrs)) / (len(scrThrs) - 1)\n",
    "scrThrsColVals = [min(scrThrs) + i * shift for i in range(len(scrThrs))]\n",
    "padding = (scrThrsColVals[1] - scrThrsColVals[0]) / 2\n",
    "\n",
    "for th in scrThrs:\n",
    "    precs.append(tpfptnfnVarScrThr[th][\"TPs\"] / (tpfptnfnVarScrThr[th][\"TPs\"] + tpfptnfnVarScrThr[th][\"FPs\"]))\n",
    "    recalls.append(tpfptnfnVarScrThr[th][\"TPs\"] / (tpfptnfnVarScrThr[th][\"TPs\"] + tpfptnfnVarScrThr[th][\"FNs\"]))\n",
    "        \n",
    "cax = ax.scatter(precs, recalls, marker=dotshapes[\"ESKEMAP\"], c=scrThrs, cmap=cmapScrs, vmin=scrThrsColVals[0] - \\\n",
    "                 padding, vmax=scrThrsColVals[-1] + padding)\n",
    "cbar = fig.colorbar(cax, ticks=scrThrsColVals, label=\"ESKEMAP target confidence interval\")\n",
    "cbar.ax.set_yticklabels(scrThrs)\n",
    "precs = []\n",
    "recalls = []\n",
    "dstThrs = sorted([0.03, 0.02, 0.01])\n",
    "padding = (dstThrs[1] - dstThrs[0]) / 2\n",
    "        \n",
    "for th in dstThrs:\n",
    "    precs.append(tpfptnfnVarDstThr[th][\"TPs\"] / (tpfptnfnVarDstThr[th][\"TPs\"] + tpfptnfnVarDstThr[th][\"FPs\"]))\n",
    "    recalls.append(tpfptnfnVarDstThr[th][\"TPs\"] / (tpfptnfnVarDstThr[th][\"TPs\"] + tpfptnfnVarDstThr[th][\"FNs\"]))\n",
    "        \n",
    "cax = ax.scatter(precs, recalls, marker=dotshapes[\"edlib\"], c=dstThrs, cmap=cmapDsts, vmin=dstThrs[0] - padding,\\\n",
    "                 vmax=dstThrs[2] + padding)\n",
    "fig.colorbar(cax, ticks=dstThrs, label=\"Edlib edit distance treshold\")\n",
    "ax.set_xlabel(\"Precision\")\n",
    "ax.set_ylabel(\"Recall\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"precRecVarEvalAllPnts.pdf\", format=\"pdf\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}