aldenjava.opticalmapping.data.mappingresult.OptMapResultNode.java Source code

Java tutorial

Introduction

Here is the source code for aldenjava.opticalmapping.data.mappingresult.OptMapResultNode.java

Source

/**************************************************************************
**  OMBlast
**  Software aligning optical maps
**  
**  Version 1.0 -- September 1, 2015
**  
**  Copyright (C) 2015 by Alden Leung, All rights reserved.
**  Contact:  aldenleung@link.cuhk.edu.hk
**  Organization:  Hong Kong Bioinformatics Centre, School of Life Sciences, The
**                 Chinese University of Hong Kong, Shatin, NT,
**                 Hong Kong SAR
**  
**  This file is part of OMBlast.
**  
**  OMBlast is free software; you can redistribute it and/or 
**  modify it under the terms of the GNU General Public License 
**  as published by the Free Software Foundation; either version 
**  3 of the License, or (at your option) any later version.
**  
**  OMBlast is distributed in the hope that it will be useful,
**  but WITHOUT ANY WARRANTY; without even the implied warranty of
**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**  GNU General Public License for more details.
**  
**  You should have received a copy of the GNU General Public 
**  License along with OMBlast; if not, see 
**  <http://www.gnu.org/licenses/>.
**************************************************************************/

package aldenjava.opticalmapping.data.mappingresult;

import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashSet;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Set;

import org.apache.commons.lang.ArrayUtils;
import org.apache.commons.lang.StringUtils;

import aldenjava.common.SimpleLocation;
import aldenjava.common.SimpleLongLocation;
import aldenjava.opticalmapping.Cigar;
import aldenjava.opticalmapping.GenomicPosNode;
import aldenjava.opticalmapping.data.data.DataNode;
import aldenjava.opticalmapping.data.data.SimulationInfo;

public class OptMapResultNode {
    // Parameters
    public DataNode parentFrag;
    public GenomicPosNode mappedRegion;
    public int mappedstrand;
    public int subrefstart;
    public int subrefstop;
    public int subfragstart;
    public int subfragstop;
    public Cigar cigar;
    public double mappedscore;
    public double confidence;

    // refp indicator
    // Due to inefficient initialization, this is implemented but not used in constructor.
    public List<MatchingSignalPair> mspsQ;
    public List<MatchingSignalPair> mspsR;

    // Init
    public OptMapResultNode(DataNode f, GenomicPosNode mappedRegion, int mappedstrand, int subrefstart,
            int subrefstop, int subfragstart, int subfragstop, Cigar cigar, double mappedscore, double confidence) {
        parentFrag = f;
        this.mappedRegion = mappedRegion;
        this.mappedstrand = mappedstrand;
        this.subrefstart = subrefstart;
        this.subrefstop = subrefstop;
        this.subfragstart = subfragstart;
        this.subfragstop = subfragstop;
        this.cigar = cigar;
        this.mappedscore = mappedscore;
        this.confidence = confidence;
        // this.updateMSP();
    }

    public OptMapResultNode(OptMapResultNode r) {
        parentFrag = r.parentFrag;
        this.mappedRegion = r.mappedRegion;
        this.mappedstrand = r.mappedstrand;
        this.subrefstart = r.subrefstart;
        this.subrefstop = r.subrefstop;
        this.subfragstart = r.subfragstart;
        this.subfragstop = r.subfragstop;
        if (r.cigar != null)
            this.cigar = new Cigar(r.cigar);
        else
            this.cigar = null;
        this.mappedscore = r.mappedscore;
        this.confidence = r.confidence;
        // this.updateMSP();
    }

    public void setCigar(Cigar cigar) {
        this.cigar = cigar;
        // this.updateMSP();
    }

    @Deprecated
    public int[] getRelativeRefPos() {
        if (cigar == null)
            return null;
        if (parentFrag == null)
            return null;
        int[] relativeRefPos = new int[this.parentFrag.getTotalSegment()];
        for (int i = 0; i < relativeRefPos.length; i++)
            relativeRefPos[i] = -1;
        boolean lastMatched = false;
        int fragpos = subfragstart - mappedstrand;
        int refpos = subrefstart - 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                if (lastMatched) {
                    relativeRefPos[fragpos] = refpos;
                }
                refpos++;
                fragpos += mappedstrand;
                lastMatched = true;
                break;
            case 'I':
                // relative RefPos is still -1
                fragpos += mappedstrand;
                lastMatched = false;
                break;
            case 'D':
                refpos++;
                lastMatched = false;
                break;
            default:
                System.err.println("Warning!");
                break;
            }
        }
        return relativeRefPos;

    }

    public List<OptMapResultNode> getBreakResult(LinkedHashMap<String, DataNode> optrefmap, int meas, double ear) {
        if (parentFrag == null || optrefmap == null)
            return null;
        List<OptMapResultNode> resultList = new ArrayList<OptMapResultNode>();
        if (!isUsed()) {
            resultList.add(new OptMapResultNode(this));
            return resultList;
        }

        DataNode ref = optrefmap.get(mappedRegion.ref);
        String precigar = this.cigar.getPrecigar();
        int lastrefpos = -1;
        int lastfragpos = -1;
        int currentrefpos = subrefstart;
        int currentfragpos = subfragstart;
        int lastrefstart = -1;
        int lastfragstart = -1;
        for (char c : precigar.toCharArray())
            switch (c) {
            case 'M':
                if (lastrefpos != -1) {
                    if (lastrefstart == -1) {
                        lastrefstart = lastrefpos;
                        lastfragstart = lastfragpos;
                    }
                    // start comparison
                    long reflen = ref.length(lastrefpos, currentrefpos - 1);
                    long fraglen;
                    if (mappedstrand == 1)
                        fraglen = parentFrag.length(lastfragpos, currentfragpos - mappedstrand);
                    else if (mappedstrand == -1)
                        fraglen = parentFrag.length(currentfragpos - mappedstrand, lastfragpos);
                    else
                        fraglen = 0;

                    boolean pass = (reflen * (1 - ear) - meas < fraglen && fraglen < reflen * (1 + ear) + meas);
                    if (!pass) {
                        if (lastfragstart != lastfragpos) {
                            resultList.add(this.getSubResult(ref, lastfragstart, lastfragpos - mappedstrand)); // the current one should not be included

                        }
                        lastrefstart = -1;
                        lastfragstart = -1;
                    } else {
                        if (pass && currentrefpos == subrefstop + 1) {
                            resultList.add(this.getSubResult(ref, lastfragstart, subfragstop)); // the current one should not be included
                            lastrefstart = -1;
                            lastfragstart = -1;
                        }
                    }
                }
                lastrefpos = currentrefpos;
                lastfragpos = currentfragpos;
                currentfragpos += mappedstrand;
                currentrefpos += 1;
                break;
            case 'I':
                currentfragpos += mappedstrand;
                break;
            case 'D':
                currentrefpos += 1;
                break;
            default:
                ;
            }
        return resultList;
    }

    public List<OptMapResultNode> getBreakResult(LinkedHashMap<String, DataNode> optrefmap,
            GenomicPosNode unwantedRegion) {
        if (parentFrag == null || optrefmap == null)
            return null;
        List<OptMapResultNode> resultList = new ArrayList<OptMapResultNode>();
        if (!isUsed()) {
            resultList.add(new OptMapResultNode(this));
            return resultList;
        }

        DataNode ref = optrefmap.get(mappedRegion.ref);
        String precigar = this.cigar.getPrecigar();
        int lastrefpos = -1;
        int lastfragpos = -1;
        int currentrefpos = subrefstart;
        int currentfragpos = subfragstart;
        int lastrefstart = -1;
        int lastfragstart = -1;
        for (char c : precigar.toCharArray()) {
            switch (c) {
            case 'M':
                if (lastrefpos != -1) {
                    if (lastrefstart == -1) {
                        lastrefstart = lastrefpos;
                        lastfragstart = lastfragpos;
                    }
                    GenomicPosNode currentRegion = new GenomicPosNode(ref.name, ref.refp[lastrefstart - 1],
                            ref.refp[currentrefpos - 1]);

                    boolean pass = !(unwantedRegion.overlapSize(currentRegion) > 0);
                    if (!pass) {
                        if (lastfragstart != lastfragpos)
                            resultList.add(this.getSubResult(ref, lastfragstart, lastfragpos - mappedstrand));

                        lastrefstart = -1;
                        lastfragstart = -1;
                    } else {
                        if (pass && currentrefpos == subrefstop + 1) {
                            resultList.add(this.getSubResult(ref, lastfragstart, subfragstop));
                            lastrefstart = -1;
                            lastfragstart = -1;
                        }
                    }
                }
                lastrefpos = currentrefpos;
                lastfragpos = currentfragpos;
                currentfragpos += mappedstrand;
                currentrefpos += 1;
                break;
            case 'I':
                currentfragpos += mappedstrand;
                break;
            case 'D':
                currentrefpos += 1;
                break;
            default:
                ;
            }
        }
        return resultList;

    }

    public int getFP() {
        if (cigar != null)
            return cigar.getFP();
        else
            return -1;
    }

    public int getFN() {
        if (cigar != null)
            return cigar.getFN();
        else
            return -1;
    }

    public double getFPRate() {
        if (cigar == null)
            return -1;
        return this.getFP() / (double) (Math.abs(length(subfragstart, subfragstop)) + 1);
    }

    public double getFNRate() {
        if (cigar == null)
            return -1;
        return this.getFN() / (double) (Math.abs(subrefstop - subrefstart) + 1);
    }

    public long getMapLength() {
        if (subfragstart == -1 || subfragstop == -1)
            return -1;
        if (mappedstrand == 1)
            return parentFrag.length(subfragstart, subfragstop);
        else if (mappedstrand == -1)
            return parentFrag.length(subfragstop, subfragstart);
        else
            return 0;
    }

    public double getMapScale() {
        if (this.isUsed())
            return (getMapLength() / (double) (mappedRegion.length()));
        else
            return -1;
    }

    public double getMapSigRatio() {
        return (cigar.getMatch() / (double) (this.parentFrag.getTotalSegment() - 1));
    }

    public int[] getRefMapSignal() {
        if (cigar == null)
            return null;
        if (parentFrag == null)
            return null;

        int[] matchSignals = new int[cigar.getMatch()];
        int index = 0;
        int fragpos = subfragstart - mappedstrand;
        int refpos = subrefstart - 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                matchSignals[index++] = refpos; // need to be checked
                refpos++;
                fragpos += mappedstrand;
                break;
            case 'I':
                // relative RefPos is still -1
                fragpos += mappedstrand;
                break;
            case 'D':
                refpos++;
                break;
            default:
                throw new RuntimeException(cigar + " character");
                // break;
            }
        }
        return matchSignals;
    }

    public int[] getMapSignal() {
        if (cigar == null)
            return null;
        if (parentFrag == null)
            return null;

        int[] matchSignals = new int[cigar.getMatch()];
        int index = 0;
        int fragpos = subfragstart - mappedstrand;
        int refpos = subrefstart - 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                matchSignals[index++] = fragpos; // need to be checked
                refpos++;
                fragpos += mappedstrand;
                break;
            case 'I':
                // relative RefPos is still -1
                fragpos += mappedstrand;
                break;
            case 'D':
                refpos++;
                break;
            default:
                throw new RuntimeException(cigar + " character");
                // break;
            }
        }
        return matchSignals;

    }

    public int getMatch() {
        if (cigar != null)
            return cigar.getMatch();
        else
            return -1;
    }

    public GenomicPosNode getMoleMappedRegion() {
        // Not-reversed
        if (mappedstrand == 1)
            return new GenomicPosNode(parentFrag.name, parentFrag.refp[subfragstart - 1],
                    parentFrag.refp[subfragstop]);
        else
            return new GenomicPosNode(parentFrag.name, parentFrag.refp[subfragstop - 1],
                    parentFrag.refp[subfragstart]);
    }

    public int getSignal(String refname, int targetSig) {
        if (!mappedRegion.ref.equals(refname))
            return -1;
        int currRefSig = subrefstart - 1;
        int currFragSig = subfragstart;

        if (mappedstrand == 1)
            currFragSig -= 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                if (targetSig == currRefSig) {
                    return currFragSig;
                }
                currRefSig++;
                currFragSig += mappedstrand;
                break;
            case 'I':
                currFragSig += mappedstrand;
                break;
            case 'D':
                currRefSig++;
                break;
            default:
                break;
            }
        }
        return -1;
    }

    public long getSignalDisplace(String refname, int targetSig) {
        int currRefSig = subrefstart - 1;
        int currFragSig = subfragstart;

        if (mappedstrand == 1)
            currFragSig -= 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                if (targetSig == currRefSig) {
                    if (mappedstrand == 1)
                        if (subfragstart > currFragSig)
                            return 0;
                        else
                            return length(subfragstart, currFragSig);
                    else if (currFragSig + 1 > subfragstart)
                        return 0;
                    else
                        return length(currFragSig + 1, subfragstart);
                }
                currRefSig++;
                currFragSig += mappedstrand;
                break;
            case 'I':
                currFragSig += mappedstrand;
                break;
            case 'D':
                currRefSig++;
                break;
            default:
                break;
            }
        }
        return 0;
    }

    public double getSubFragRatio() {
        if (getMapLength() == -1)
            return -1;
        return (getMapLength() / (double) (length(1, this.getTotalSegment() - 2)));
    }

    public OptMapResultNode getSubResult(DataNode ref, int fragstart, int fragstop) {
        if ((mappedstrand == 1 && subfragstart <= fragstart && fragstart <= subfragstop && subfragstart <= fragstop
                && fragstop <= subfragstop)
                || (mappedstrand == -1 && subfragstop <= fragstart && fragstart <= subfragstart
                        && subfragstop <= fragstop && fragstop <= subfragstart)) {
            String precigar = this.cigar.getPrecigar();
            int currentrefpos = subrefstart;
            int currentfragpos = subfragstart;

            int newfragstart = -1;
            int newfragstop = -1;
            int newrefstart = -1;
            int newrefstop = -1;

            StringBuilder storedPrecigar = new StringBuilder();
            StringBuilder recentPrecigar = new StringBuilder();
            boolean start = false;
            boolean stop = false;
            for (char c : precigar.toCharArray()) {
                recentPrecigar.append(c);
                switch (c) {

                case 'M':
                    if (fragstart * mappedstrand <= currentfragpos * mappedstrand) // we need direction
                    {
                        start = true;
                        if (newrefstart == -1) {
                            newrefstart = currentrefpos;
                            newfragstart = currentfragpos;
                        }
                    }
                    if ((fragstop + mappedstrand) * mappedstrand <= currentfragpos * mappedstrand) {
                        stop = true;
                        if (fragstop + mappedstrand == currentfragpos) {
                            // if fragstop + mappedstrand < currentfragpos, some extra unmatches occur
                            storedPrecigar.append(c);
                            newrefstop = currentrefpos - 1;
                            newfragstop = currentfragpos - mappedstrand;
                        } else
                            // Error occurs if you are using exact Match Position to draw sub result
                            ;

                        break;
                    } else {
                        // continue to update stop position
                        newrefstop = currentrefpos - 1; // correct?
                        newfragstop = currentfragpos - mappedstrand;
                    }
                    currentfragpos += mappedstrand;
                    currentrefpos += 1;
                    break;
                case 'I':
                    currentfragpos += mappedstrand;
                    break;
                case 'D':
                    currentrefpos += 1;
                    break;
                default:
                    ;
                }
                if (stop)
                    break;
                if (start)
                    storedPrecigar.append(c);
            }
            Cigar newcigar;
            newcigar = new Cigar(storedPrecigar.toString());
            return (new OptMapResultNode(parentFrag, ref.getGenomicPos(newrefstart, newrefstop, true), mappedstrand,
                    newrefstart, newrefstop, newfragstart, newfragstop, newcigar, -1, -1));
        } else {
            throw new IllegalArgumentException(
                    "Sub-result argument fails boundary checking. " + parentFrag == null ? ""
                            : ("Result ID: " + parentFrag.name));

        }

    }

    public int getSubRefSigStart() {
        if (subrefstart < 0)
            return -1;
        return subrefstart - 1;
    }

    public int getSubRefSigStop() {
        if (subrefstop < 0)
            return -1;
        return subrefstop;
    }

    public int getSubMoleSigStart() {
        if (subfragstart < 0)
            return -1;
        if (mappedstrand == 1)
            return subfragstart - 1;
        else
            return subfragstart;

    }

    public int getSubMoleSigStop() {
        if (subfragstop < 0)
            return -1;
        if (mappedstrand == 1)
            return subfragstop;
        else
            return subfragstop - 1;
    }

    public void updateMSP() {
        if (cigar == null || subfragstart < 0 || subfragstop < 0 || subrefstart < 0 || subrefstop < 0)
            return;
        List<MatchingSignalPair> msps = new ArrayList<MatchingSignalPair>();

        int currRefSig = subrefstart - 1;
        int currQuerySig = subfragstart;

        if (mappedstrand == 1)
            currQuerySig -= 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                msps.add(new MatchingSignalPair(currRefSig, currQuerySig));
                currRefSig++;
                currQuerySig += mappedstrand;
                break;
            case 'I':
                currQuerySig += mappedstrand;
                break;
            case 'D':
                currRefSig++;
                break;
            default:
                // Other characters are not supported
                break;
            }
        }
        List<MatchingSignalPair> mspsQ = new ArrayList<MatchingSignalPair>(msps);
        Collections.sort(msps, MatchingSignalPair.qposComparator);
        List<MatchingSignalPair> mspsR = new ArrayList<MatchingSignalPair>(msps);
        Collections.sort(msps, MatchingSignalPair.rposComparator);

        this.mspsQ = mspsQ;
        this.mspsR = mspsR;
    }

    public Set<MatchingSignalPair> getMSPOverlapped(OptMapResultNode result) {
        Set<MatchingSignalPair> mspSet = new HashSet<MatchingSignalPair>();

        int i1 = 0;
        int i2 = 0;
        while (i1 < this.mspsQ.size() && i2 < result.mspsQ.size())
            if (this.mspsQ.get(i1).qpos == result.mspsQ.get(i2).qpos) {
                mspSet.add(this.mspsQ.get(i1));
                i1++;
                i2++;
            } else if (this.mspsQ.get(i1).qpos > result.mspsQ.get(i2).qpos)
                i2++;
            else
                i1++;

        int j1 = 0;
        int j2 = 0;
        while (j1 < this.mspsR.size() && j2 < result.mspsR.size())
            if (this.mspsR.get(j1).rpos == result.mspsR.get(j2).rpos) {
                mspSet.add(this.mspsR.get(j1));
                j1++;
                j2++;
            } else if (this.mspsR.get(j1).rpos > result.mspsR.get(j2).rpos)
                j2++;
            else
                j1++;
        return mspSet;
    }

    public OptMapResultNode getSubResultbyRefPos(DataNode ref, int refstart, int refstop) {
        // input is signal start and stop
        // pter points at every signal
        int[] pter = new int[parentFrag.getTotalSegment() - 1];

        String precigar = this.cigar.getPrecigar();
        int currentrefpos = subrefstart;
        int currentfragpos = subfragstart;

        if (mappedstrand == 1)
            for (int i = 0; i < subfragstart - 1; i++)
                pter[i] = -1;
        else
            for (int i = subfragstart; i < parentFrag.getTotalSegment() - 1; i++)
                pter[i] = -1;
        for (char c : precigar.toCharArray()) {
            switch (c) {
            case 'M':
                if (mappedstrand == 1)
                    pter[currentfragpos - 1] = currentrefpos - 1;
                else
                    pter[currentfragpos] = currentrefpos - 1;
                currentfragpos += mappedstrand;
                currentrefpos += 1;
                break;
            case 'I':
                if (mappedstrand == 1)
                    pter[currentfragpos - 1] = -1;
                else
                    pter[currentfragpos] = -1;
                currentfragpos += mappedstrand;
                break;
            case 'D':
                currentrefpos += 1;
                break;
            default:
                ;
            }
        }
        if (mappedstrand == 1)
            for (int i = subfragstop + 1; i < parentFrag.getTotalSegment() - 1; i++)
                pter[i] = -1;
        else
            for (int i = 0; i < subfragstop - 1; i++)
                pter[i] = -1;

        int targetfragstart = -1;
        int targetfragstop = -1;
        for (int i = 0; i < parentFrag.getTotalSegment() - 1; i++) {
            if (pter[i] == refstart)
                targetfragstart = i;
            if (pter[i] == refstop)
                targetfragstop = i;
        }
        if (targetfragstart == -1 || targetfragstop == -1)
            return null;
        else {
            if (mappedstrand == 1)
                return getSubResult(ref, targetfragstart + 1, targetfragstop);
            else
                return getSubResult(ref, targetfragstart, targetfragstop + 1);
        }

    }

    public int getTotalSegment() {
        return parentFrag.getTotalSegment();
    }

    public long length() {
        return parentFrag.length();
    }

    public long length(int start, int stop) {
        return parentFrag.length(start, stop);
    }

    @Override
    public String toString() {
        // Name, Aligned Region, Cigar
        if (this.isUsed())
            return String.format("%s\t%s\t%s", parentFrag.toString(), mappedRegion.toString(), cigar);
        else
            return parentFrag.toString();
    }

    // Modify
    public void reverse() {
        parentFrag = parentFrag.getReverse();
        this.mappedstrand *= -1;
        this.subfragstart = parentFrag.getTotalSegment() - subfragstart - 1;
        this.subfragstop = parentFrag.getTotalSegment() - subfragstop - 1;
    }

    public void scale(double ratio) {
        parentFrag = new DataNode(parentFrag);
        parentFrag.scale(ratio);
        // super.scale(ratio);
    }

    private void trimResult(int step) {
        String precigar = cigar.getPrecigar();
        if (step == 0) {
            // updateMSP();
            return;
        } else if (step < 0) // right to left trim, according to ref
            precigar = StringUtils.reverse(precigar);
        if (precigar.length() == 1) // will trim to nothing, cannot trim anymore
            return;
        String removed = precigar.substring(0, precigar.indexOf('M', 1));
        ;
        int m = 0;
        int i = 0;
        int d = 0;
        for (char c : removed.toCharArray())
            switch (c) {
            case 'M':
                m++;
                break;
            case 'I':
                i++;
                break;
            case 'D':
                d++;
                break;
            default:
                break;
            }
        if (step > 0) {
            subrefstart += m + d;
            if (mappedstrand == -1)
                subfragstart -= m + i; // + or -
            else
                subfragstart += m + i;
        } else {
            subrefstop -= m + d; // + or -
            if (mappedstrand == -1)
                subfragstop += m + i;
            else
                subfragstop -= m + i; // + or -
        }
        String newprecigar = precigar.substring(precigar.indexOf('M', 1), precigar.length());
        if (step < 0)
            newprecigar = StringUtils.reverse(newprecigar);
        this.cigar = new Cigar(newprecigar);
        trimResult(step - Math.abs(step) / step);
    }

    public void trimResult(int step, LinkedHashMap<String, DataNode> optrefmap) {
        this.trimResult(step, optrefmap.get(mappedRegion.ref));
    }

    public void trimResult(int step, DataNode ref) {
        this.trimResult(step);
        this.updateMappedRegion(ref);
    }

    public void updateScore(LinkedHashMap<String, DataNode> optrefmap, int match, int fpp, int fnp) {
        this.mappedscore = (cigar.calcScore(match, fpp, fnp)) * (1 - Math.abs(1 - getMapScale()));
    }

    public void updateMappedRegion(DataNode ref) {
        long mappedstart;
        long mappedstop;
        if (subrefstart >= 1)
            mappedstart = ref.refp[subrefstart - 1];
        else
            mappedstart = 0;
        if (subrefstop < ref.refp.length)
            mappedstop = ref.refp[subrefstop];
        else
            mappedstop = ref.size;
        mappedRegion = new GenomicPosNode(ref.name, mappedstart, mappedstop);
    }

    // Test
    public boolean correctlyMapped() {
        // double defaultMaxMappedFactor = 0.50;
        // return correctlyMapped(defaultMaxMappedFactor);
        if (parentFrag == null)
            return false;
        if (!parentFrag.hasSimulationInfo())
            return false;
        return correctlyMapped(0, parentFrag.simuInfo);
    }

    public boolean correctlyMapped(SimulationInfo simuInfo) {
        return correctlyMapped(0, simuInfo);
    }

    public boolean correctlyMapped(long tolerance, SimulationInfo simuInfo) {
        return mappedRegion.isClose(simuInfo.simuRegion, tolerance) && mappedstrand == simuInfo.simuStrand;
    }

    public boolean isRefClose(OptMapResultNode map, long gapAllowed) {
        return this.mappedRegion.isClose(map.mappedRegion, gapAllowed);
    }

    public boolean isFragClose(OptMapResultNode map, long gapAllowed) {
        SimpleLongLocation s1;
        if (mappedstrand == 1)
            s1 = new SimpleLongLocation(length(0, subfragstart - 1), length(0, subfragstop));
        else if (mappedstrand == -1)
            s1 = new SimpleLongLocation(length(0, subfragstop - 1), length(0, subfragstart));
        else
            s1 = null;
        SimpleLongLocation s2;
        if (map.mappedstrand == 1)
            s2 = new SimpleLongLocation(length(0, map.subfragstart - 1), length(0, map.subfragstop));
        else if (map.mappedstrand == -1)
            s2 = new SimpleLongLocation(length(0, map.subfragstop - 1), length(0, map.subfragstart));
        else
            s2 = null;

        return (s1.overlap(s2, gapAllowed));

    }

    public boolean isClose(OptMapResultNode map, long gapAllowed) {
        return (this.isFragClose(map, gapAllowed) && this.isRefClose(map, gapAllowed));
    }

    public boolean isClose(GenomicPosNode region, long gapAllowed) {
        return this.mappedRegion.isClose(region, gapAllowed);
    }

    public boolean isSimilar(OptMapResultNode map, long gapAllowed) {
        return isClose(map, gapAllowed) && (this.mappedstrand == map.mappedstrand);
    }

    public boolean isUsed() {
        if (mappedRegion == null)
            return false;
        else if (mappedRegion.ref.equalsIgnoreCase("Discarded"))
            return false;
        else if (mappedRegion.ref.equalsIgnoreCase("Unmapped"))
            return false;
        return true;
    }

    public Boolean isSubFragInfoValid() {
        if (parentFrag.refp == null) {
            System.err.println("No fragment information for checking.");
            return null;
        }
        int max = parentFrag.getTotalSegment() - 1;
        int min = 0;
        return (subfragstart <= max && subfragstart >= min && subfragstop <= max && subfragstop >= min);
    }

    public boolean isSubRefInfoValid(DataNode ref) {
        int max = ref.refp.length;
        int min = 0;
        return (subrefstart <= max && subrefstart >= min && subrefstop <= max && subrefstop >= min)
                && (subrefstart <= subrefstop);

    }

    public boolean isSubRefInfoValid(LinkedHashMap<String, DataNode> optrefmap) {
        if (optrefmap.containsKey(mappedRegion.ref))
            return isSubRefInfoValid(optrefmap.get(mappedRegion.ref));
        else
            return false;
    }

    public boolean mapSignal(String refname, int targetSig) {
        if (!mappedRegion.ref.equalsIgnoreCase(refname))
            return false;
        int currRefSig = subrefstart - 1;
        for (char c : cigar.getPrecigar().toCharArray()) {
            switch (c) {
            case 'M':
                if (targetSig == currRefSig)
                    return true;
                currRefSig++;
                break;
            case 'I':
                break;
            case 'D':
                currRefSig++;
                break;
            default:
                break;
            }
            if (currRefSig > targetSig)
                return false;
        }
        return false;
    }

    public boolean overlap(OptMapResultNode result) {
        return (overlapQuery(result) || overlapRef(result));
    }

    public boolean overlapQuery(OptMapResultNode result) {
        SimpleLocation thisQueryLoc = new SimpleLocation(this.subfragstart, this.subfragstop);
        SimpleLocation resultQueryloc = new SimpleLocation(result.subfragstart, result.subfragstop);

        return (thisQueryLoc.overlapsize(resultQueryloc) > -1); // even subfragstart not overlap is not enough; the signal is still shared.
    }

    public boolean overlapRef(OptMapResultNode map) {
        SimpleLocation thisrefloc = new SimpleLocation(this.subrefstart, this.subrefstop);
        SimpleLocation maprefloc = new SimpleLocation(map.subrefstart, map.subrefstop);
        return ((this.mappedRegion.ref.equals(map.mappedRegion.ref) && thisrefloc.overlapsize(maprefloc) > -1)); // even subrefstart not overlap is not enough; the signal is still shared.
    }

    // Copy
    public List<OptMapResultNode> getRealMap() {
        List<OptMapResultNode> mapList = new ArrayList<OptMapResultNode>();
        mapList.add(this);
        return mapList;
    }

    public OptMapResultNode newReverseFragment() {
        OptMapResultNode map = new OptMapResultNode(this);
        map.reverse();
        return map;
    }

    // Static
    public static OptMapResultNode newBlankMapNode(DataNode f) {
        return new OptMapResultNode(f, null, 0, -1, -1, -1, -1, null, Double.NEGATIVE_INFINITY, -1);
    }

    public static List<GenomicPosNode> getMappedRegion(List<OptMapResultNode> resultList) {
        if (resultList == null)
            throw new NullPointerException("resultList");
        List<GenomicPosNode> posList = new ArrayList<>();
        for (OptMapResultNode result : resultList)
            posList.add(result.mappedRegion);
        return posList;
    }

    public static List<GenomicPosNode> getMoleMappedRegion(List<OptMapResultNode> resultList) {
        if (resultList == null)
            throw new NullPointerException("resultList");
        List<GenomicPosNode> posList = new ArrayList<>();
        for (OptMapResultNode result : resultList)
            posList.add(result.getMoleMappedRegion());
        return posList;
    }

    public static List<GenomicPosNode> getPotentiallyMappedRegion(LinkedHashMap<String, DataNode> optrefmap,
            List<OptMapResultNode> resultList) {
        List<GenomicPosNode> targetRegionList = new ArrayList<GenomicPosNode>();
        if (resultList.isEmpty())
            return targetRegionList;
        if (!resultList.get(0).isUsed())
            return targetRegionList;

        for (OptMapResultNode result : resultList) {
            String ref = result.mappedRegion.ref;
            long start = result.mappedRegion.start - result.parentFrag.size;
            long stop = result.mappedRegion.stop + result.parentFrag.size;
            if (start < 1)
                start = 1;
            if (stop > optrefmap.get(ref).size)
                stop = optrefmap.get(ref).size;
            targetRegionList.add(new GenomicPosNode(ref, start, stop));
        }
        return GenomicPosNode.merge(targetRegionList);
    }

    public static LinkedHashMap<String, List<GenomicPosNode>> getPotentiallyMappedRegion(
            LinkedHashMap<String, DataNode> optrefmap,
            LinkedHashMap<String, List<OptMapResultNode>> resultListMap) {

        LinkedHashMap<String, List<GenomicPosNode>> targetRegionMap = new LinkedHashMap<>();
        for (List<OptMapResultNode> resultList : resultListMap.values()) {
            List<GenomicPosNode> targetRegionList = getPotentiallyMappedRegion(optrefmap, resultList);
            targetRegionMap.put(resultList.get(0).parentFrag.name, targetRegionList);
        }
        return targetRegionMap;
    }

    public static int[] getMapSignal(List<OptMapResultNode> resultList) {
        Set<Integer> set = new HashSet<Integer>();
        for (OptMapResultNode result : resultList) {
            for (int i : result.getMapSignal())
                set.add(i);
        }
        List<Integer> list = new ArrayList<Integer>(set);
        Collections.sort(list);
        return ArrayUtils.toPrimitive(list.toArray(new Integer[list.size()]));
    }

    public static List<SimpleLocation> getUnmapPortion(List<OptMapResultNode> fragmentmaplist) {
        if (fragmentmaplist == null || fragmentmaplist.isEmpty())
            return null;
        DataNode f = fragmentmaplist.get(0).parentFrag;
        // boolean[] usedSegment = new boolean[f.getTotalSegment()];
        List<SimpleLocation> subfraglist = new ArrayList<SimpleLocation>();
        for (OptMapResultNode map : fragmentmaplist)
            subfraglist.add(new SimpleLocation(map.subfragstart, map.subfragstop));
        SimpleLocation.mergeLocationList(subfraglist);
        List<SimpleLocation> remainfraglist = new ArrayList<SimpleLocation>();
        remainfraglist.add(new SimpleLocation(0, f.getTotalSegment() - 1));
        SimpleLocation.removeSimpleLocationList(remainfraglist, subfraglist);
        return remainfraglist;
    }

    public static List<OptMapResultNode> reconstruct(LinkedHashMap<String, DataNode> optrefmap,
            OptMapResultNode map) {
        // Only reconstruct the subrefstart,.... and subfragratio, and cigar string
        int direction = map.mappedstrand;
        List<OptMapResultNode> fragmentmaplist = new ArrayList<OptMapResultNode>();

        // String precigar = Cigar.convertpreCIGAR(map.cigar);
        String precigar = map.cigar.getPrecigar();
        // if (direction == -1)
        // precigar = StringUtils.reverse(precigar);
        String[] cigarlist = precigar.split("S");
        for (int i = 1; i <= cigarlist.length; i += 2) {
            int subrefstart = -1;
            int subrefstop = -1;
            int subfragstart = -1;
            int subfragstop = -1;
            // double scale = 0;
            if (i - 2 < 0) {
                subrefstart = map.subrefstart;
                subfragstart = map.subfragstart;
            } else {
                StringBuilder concat = new StringBuilder();
                for (int j = 0; j < i - 1; j++)
                    concat.append(cigarlist[j]);
                String previousprecigar = concat.toString();
                int match = Cigar.getCertainNumberFromPrecigar(previousprecigar, 'M');
                int insert = Cigar.getCertainNumberFromPrecigar(previousprecigar, 'I');
                int delete = Cigar.getCertainNumberFromPrecigar(previousprecigar, 'D');
                subfragstart = map.subfragstart + (match + insert) * direction;
                subrefstart = map.subrefstart + match + delete;
            }
            String recentprecigar = cigarlist[i - 1];
            if (i == cigarlist.length) {
                subrefstop = map.subrefstop;
                subfragstop = map.subfragstop;
            } else {
                int match = Cigar.getCertainNumberFromPrecigar(recentprecigar, 'M');
                int insert = Cigar.getCertainNumberFromPrecigar(recentprecigar, 'I');
                int delete = Cigar.getCertainNumberFromPrecigar(recentprecigar, 'D');
                subfragstop = subfragstart + (match + insert - 1) * direction;
                subrefstop = subrefstart + match + delete - 1;
            }
            DataNode ref = optrefmap.get(map.mappedRegion.ref);
            long estimatestartpos = ref.refp[subrefstart - 1];
            long estimatestoppos = ref.refp[subrefstop];

            Cigar newcigar = new Cigar(recentprecigar);
            fragmentmaplist.add(new OptMapResultNode(map.parentFrag,
                    new GenomicPosNode(map.mappedRegion.ref, estimatestartpos, estimatestoppos), map.mappedstrand,
                    subrefstart, subrefstop, subfragstart, subfragstop, newcigar, map.mappedscore, -1));
        }
        return fragmentmaplist;
    }

    public static OptMapResultNode reverseRefAndFrag(OptMapResultNode result, DataNode originalRef) {
        DataNode newRef = new DataNode(result.parentFrag);
        OptMapResultNode newResult;
        int newsubfragstart;
        int newsubfragstop;
        int newsubrefstart;
        int newsubrefstop;
        Cigar newcigar;
        if (result.mappedstrand == 1) {
            newsubfragstart = result.subrefstart;
            newsubfragstop = result.subrefstop;
            newsubrefstart = result.subfragstart;
            newsubrefstop = result.subfragstop;
            newcigar = new Cigar(result.cigar);
        } else {
            newsubfragstart = result.subrefstop;
            newsubfragstop = result.subrefstart;
            newsubrefstart = result.subfragstop;
            newsubrefstop = result.subfragstart;
            newcigar = result.cigar.getReverseCigar();
        }
        newResult = new OptMapResultNode(new DataNode(originalRef), null, result.mappedstrand, newsubrefstart,
                newsubrefstop, newsubfragstart, newsubfragstop, newcigar, result.mappedscore, result.confidence);
        newResult.updateMappedRegion(newRef);
        return newResult;
    }

    public static OptMapResultNode reverseRefAndFrag(OptMapResultNode result,
            LinkedHashMap<String, DataNode> originalOptRefMap) {
        return OptMapResultNode.reverseRefAndFrag(result, originalOptRefMap.get(result.mappedRegion.ref));
    }

    public static List<OptMapResultNode> reverseRefAndFrag(List<OptMapResultNode> resultList,
            DataNode originalRef) {
        List<OptMapResultNode> newResultList = new ArrayList<OptMapResultNode>();
        for (OptMapResultNode result : resultList)
            newResultList.add(OptMapResultNode.reverseRefAndFrag(result, originalRef));
        return newResultList;
    }

    public static List<OptMapResultNode> reverseRefAndFrag(List<OptMapResultNode> resultList,
            LinkedHashMap<String, DataNode> originalOptRefMap) {
        List<OptMapResultNode> newResultList = new ArrayList<OptMapResultNode>();
        for (OptMapResultNode result : resultList)
            newResultList.add(OptMapResultNode.reverseRefAndFrag(result, originalOptRefMap));
        return newResultList;
    }

    public static Comparator<OptMapResultNode> mappedstartcomparator = new Comparator<OptMapResultNode>() {
        @Override
        public int compare(OptMapResultNode f1, OptMapResultNode f2) {
            return f1.mappedRegion.compareTo(f2.mappedRegion);
        }
    };
    public static Comparator<OptMapResultNode> subfragstartcomparator = new Comparator<OptMapResultNode>() {
        @Override
        public int compare(OptMapResultNode f1, OptMapResultNode f2) {
            int i = Integer.valueOf(f1.subfragstart).compareTo(Integer.valueOf(f2.subfragstart));
            if (i == 0)
                i = Integer.valueOf(f1.subfragstop).compareTo(Integer.valueOf(f2.subfragstop));
            return i;
        }
    };
    public static Comparator<OptMapResultNode> subfragstartstopcomparator = new Comparator<OptMapResultNode>() {
        // consider the comparison when reverse and forward strand is together
        @Override
        public int compare(OptMapResultNode f1, OptMapResultNode f2) {
            int i = Integer.valueOf(f1.mappedstrand == 1 ? f1.subfragstart : f1.subfragstop)
                    .compareTo(Integer.valueOf(f2.mappedstrand == 1 ? f2.subfragstart : f2.subfragstop));
            if (i == 0)
                i = Integer.valueOf(f1.mappedstrand == 1 ? f1.subfragstop : f1.subfragstart)
                        .compareTo(Integer.valueOf(f2.mappedstrand == 1 ? f2.subfragstop : f2.subfragstart));
            return i;
        }
    };
    public static Comparator<OptMapResultNode> mappedscorecomparator = new Comparator<OptMapResultNode>() {
        @Override
        public int compare(OptMapResultNode f1, OptMapResultNode f2) {
            return Double.valueOf(f1.mappedscore).compareTo(Double.valueOf(f2.mappedscore));
        }
    };

}

class MatchingSignalPair {
    public final int rpos;
    public final int qpos;

    public MatchingSignalPair(int rpos, int qpos) {
        this.rpos = rpos;
        this.qpos = qpos;
    }

    public boolean equals(MatchingSignalPair msp) {
        return this.rpos == msp.rpos && this.qpos == msp.qpos;
    }

    public boolean overlap(MatchingSignalPair msp) {
        return this.rpos == msp.rpos || this.qpos == msp.qpos;
    }

    // In a list of MSPs in an alignment, qpos and rpos should be unique
    public static Comparator<MatchingSignalPair> rposComparator = new Comparator<MatchingSignalPair>() {
        @Override
        public int compare(MatchingSignalPair msp1, MatchingSignalPair msp2) {
            return Integer.compare(msp1.rpos, msp2.rpos);
        }
    };
    public static Comparator<MatchingSignalPair> qposComparator = new Comparator<MatchingSignalPair>() {
        @Override
        public int compare(MatchingSignalPair msp1, MatchingSignalPair msp2) {
            return Integer.compare(msp1.qpos, msp2.qpos);
        }
    };
}