[med-svn] [sniffles] 01/02: Imported Upstream version 0.0.1+ds
Afif Elghraoui
afif at moszumanska.debian.org
Sat Jun 18 23:28:36 UTC 2016
This is an automated email from the git hooks/post-receive script.
afif pushed a commit to branch master
in repository sniffles.
commit 33452fbffafdd0410ebe01697b8d3397506540a2
Author: Afif Elghraoui <afif at ghraoui.name>
Date: Sat Jun 18 13:43:38 2016 -0700
Imported Upstream version 0.0.1+ds
---
CMakeLists.txt | 36 +++
LICENSE | 22 ++
README.md | 92 ++++++
cmake/Sniffles | Bin 0 -> 261820 bytes
cmake/makefile | 62 ++++
cmake/objects.mk | 8 +
cmake/sources.mk | 30 ++
cmake/src/Alignment.d | 54 ++++
cmake/src/Alignment.o | Bin 0 -> 628428 bytes
cmake/src/BamParser.d | 58 ++++
cmake/src/BamParser.o | Bin 0 -> 298480 bytes
cmake/src/Parser.d | 56 ++++
cmake/src/Parser.o | Bin 0 -> 8376 bytes
cmake/src/Sniffles.d | 148 +++++++++
cmake/src/Sniffles.o | Bin 0 -> 812932 bytes
cmake/src/plane-sweep/Main.d | 1 +
cmake/src/plane-sweep/Main.o | Bin 0 -> 1380 bytes
cmake/src/plane-sweep/MyHeap.d | 62 ++++
cmake/src/plane-sweep/MyHeap.o | Bin 0 -> 179100 bytes
cmake/src/plane-sweep/MyList.d | 62 ++++
cmake/src/plane-sweep/MyList.o | Bin 0 -> 91880 bytes
cmake/src/plane-sweep/PlaneSweep.d | 67 ++++
cmake/src/plane-sweep/PlaneSweep.o | Bin 0 -> 125776 bytes
cmake/src/plane-sweep/subdir.mk | 33 ++
cmake/src/sub/Detect_Breapoints.d | 82 +++++
cmake/src/sub/Detect_Breapoints.o | Bin 0 -> 877496 bytes
cmake/src/sub/subdir.mk | 24 ++
cmake/src/subdir.mk | 33 ++
cmake/src/tree/IntervallTree.d | 67 ++++
cmake/src/tree/IntervallTree.o | Bin 0 -> 297288 bytes
cmake/src/tree/Scapegoat.d | 66 ++++
cmake/src/tree/Scapegoat.o | Bin 0 -> 8688 bytes
cmake/src/tree/subdir.mk | 27 ++
src/Alignment.cpp | 607 +++++++++++++++++++++++++++++++++++++
src/Alignment.h | 110 +++++++
src/BamParser.cpp | 56 ++++
src/BamParser.h | 42 +++
src/CMakeLists.txt | 60 ++++
src/Ignore_Regions.cpp | 112 +++++++
src/Ignore_Regions.h | 15 +
src/Paramer.h | 93 ++++++
src/Parser.cpp | 8 +
src/Parser.h | 21 ++
src/Sniffles.cpp | 178 +++++++++++
src/Version.h | 9 +
src/Version.h.in | 9 +
src/phasing/PhaserSV.cpp | 13 +
src/phasing/PhaserSV.h | 24 ++
src/plane-sweep/IContainer.h | 24 ++
src/plane-sweep/Main.cpp | 81 +++++
src/plane-sweep/MyHeap.cpp | 93 ++++++
src/plane-sweep/MyHeap.h | 31 ++
src/plane-sweep/MyList.cpp | 68 +++++
src/plane-sweep/MyList.h | 51 ++++
src/plane-sweep/Node.h | 79 +++++
src/plane-sweep/Plane-sweep.h | 58 ++++
src/plane-sweep/PlaneSweep.cpp | 43 +++
src/print/BedePrinter.cpp | 59 ++++
src/print/BedePrinter.h | 26 ++
src/print/IPrinter.cpp | 51 ++++
src/print/IPrinter.h | 65 ++++
src/print/MariaPrinter.cpp | 58 ++++
src/print/MariaPrinter.h | 24 ++
src/print/NGMPrinter.cpp | 70 +++++
src/print/NGMPrinter.h | 30 ++
src/print/VCFPrinter.cpp | 84 +++++
src/print/VCFPrinter.h | 29 ++
src/realign/IAlignment.h | 63 ++++
src/realign/Realign.cpp | 172 +++++++++++
src/realign/Realign.h | 40 +++
src/realign/SWCPU.cpp | 563 ++++++++++++++++++++++++++++++++++
src/realign/SWCPU.h | 119 ++++++++
src/sub/Breakpoint.cpp | 412 +++++++++++++++++++++++++
src/sub/Breakpoint.h | 149 +++++++++
src/sub/Container.h | 15 +
src/sub/Detect_Breakpoints.cpp | 465 ++++++++++++++++++++++++++++
src/sub/Detect_Breakpoints.h | 29 ++
src/sub/IRegion.h | 56 ++++
src/tree/BinTree.cpp | 293 ++++++++++++++++++
src/tree/BinTree.h | 42 +++
src/tree/IntervallTree.cpp | 269 ++++++++++++++++
src/tree/IntervallTree.h | 40 +++
src/tree/Intervall_bed.cpp | 245 +++++++++++++++
src/tree/Intervall_bed.h | 39 +++
src/tree/Leaf.h | 64 ++++
src/tree/Scapegoat.cpp | 10 +
src/tree/Scapegoat.h | 234 ++++++++++++++
src/tree/TNode.h | 68 +++++
88 files changed, 6858 insertions(+)
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..65d74f9
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 2.8)
+project(Sniffles)
+
+set( NGM_VERSION_MAJOR 0 )
+set( NGM_VERSION_MINOR 1 )
+IF(CMAKE_BUILD_TYPE STREQUAL "Debug")
+ message(STATUS "Building in debug mode!")
+ set( NGM_VERSION_BUILD 0-debug )
+else()
+ set( NGM_VERSION_BUILD 0 )
+ENDIF()
+
+
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin/sniffles-${NGM_VERSION_MAJOR}.${NGM_VERSION_MINOR}.${NGM_VERSION_BUILD}/)
+file(MAKE_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+
+
+# Set a default build type for single-configuration
+# CMake generators if no build type is set.
+IF(NOT CMAKE_CONFIGURATION_TYPES AND NOT CMAKE_BUILD_TYPE)
+ message(STATUS "No build type specified. Using 'release'")
+ SET(CMAKE_BUILD_TYPE Release)
+ENDIF(NOT CMAKE_CONFIGURATION_TYPES AND NOT CMAKE_BUILD_TYPE)
+
+FIND_PACKAGE(OpenMP REQUIRED)
+if(OPENMP_FOUND)
+ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+ set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+endif()
+
+
+add_subdirectory(lib/zlib-1.2.7)
+add_subdirectory(lib/bamtools-2.3.0)
+
+add_subdirectory(src)
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..ac24c9a
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2015 Fritz Sedlazeck
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..2c6d9e0
--- /dev/null
+++ b/README.md
@@ -0,0 +1,92 @@
+# Sniffles
+Sniffles is a structural variation caller using third generation sequencing (PacBio or Oxford Nanopore). It detects all types of SVs using evidence from split-read alignments, high-mismatch regions, and coverage analysis. Please note the current version of Sniffles requires output from BWA-MEM with the optional SAM attributes enabled! If you experience problems or have suggestions please contact: fritz.sedlazeck at gmail.com
+
+**************************************
+
+INSTALL:
+
+Download Sniffles:
+```
+git clone https://github.com/fritzsedlazeck/Sniffles
+```
+
+ cd Sniffles
+
+ mkdir build
+
+ cd build
+
+ cmake ..
+
+ make
+
+ cd ../bin/
+
+
+**************************************
+
+USAGE:
+```
+./sniffles -m reads.bam -v calls.vcf
+```
+
+Options:
+
+ ./sniffles -m <string> [-s <int>] [--max_num_splits <int>] [-q <int>]
+ [-l <int>] [-v <string>] [--bede <string>] [-c <int>] [-t
+ <int>] [-d <int>] [-n <int>] [--use_MD_Cigar] [--]
+ [--version] [-h]
+
+
+Where:
+
+ -m <string>, --mapped_reads <string>
+ (required) Bam File
+
+ -s <int>, --min_support <int>
+ Minimum number of reads that support a SV. Default: 10
+
+ --max_num_splits <int>
+ Maximum number of splits per read to be still taken into account.
+ Default: 4
+
+ -q <int>, --minmapping_qual <int>
+ Minimum Mapping Quality. Default: 20
+
+ -l <int>, --min_length <int>
+ Minimum length of SV to be reported. Default:0
+
+ -v <string>, --vcf <string>
+ VCF output file name
+
+ --bede <string>
+ Simplified format of bede Format.
+
+ -c <int>, --min_cigar_event <int>
+ Minimum Cigar Event (e.g. Insertion, deletion) to take into account.
+ Default:50
+
+ -t <int>, --threads <int>
+ Number of threads to use. Default: 3
+
+ -d <int>, --max_distance <int>
+ Maximum distance to group SV together. Default: 1kb
+
+ -n <int>, --num_reads_report <int>
+ Report up to N reads that support the SV. Default: 0
+
+ --use_MD_Cigar
+ Enables Sniffles to use the alignment information to screen for
+ suspicious regions.
+
+ --, --ignore_rest
+ Ignores the rest of the labeled arguments following this flag.
+
+ --version
+ Displays version information and exits.
+
+ -h, --help
+ Displays usage information and exits.
+
+
+ Sniffles version 0.0.1
diff --git a/cmake/Sniffles b/cmake/Sniffles
new file mode 100755
index 0000000..50306ac
Binary files /dev/null and b/cmake/Sniffles differ
diff --git a/cmake/makefile b/cmake/makefile
new file mode 100644
index 0000000..73965cd
--- /dev/null
+++ b/cmake/makefile
@@ -0,0 +1,62 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+-include ../makefile.init
+
+RM := rm -rf
+
+# All of the sources participating in the build are defined here
+-include sources.mk
+-include src/tree/subdir.mk
+-include src/sub/subdir.mk
+-include src/plane-sweep/subdir.mk
+-include src/subdir.mk
+-include subdir.mk
+-include objects.mk
+
+ifneq ($(MAKECMDGOALS),clean)
+ifneq ($(strip $(CC_DEPS)),)
+-include $(CC_DEPS)
+endif
+ifneq ($(strip $(C++_DEPS)),)
+-include $(C++_DEPS)
+endif
+ifneq ($(strip $(C_UPPER_DEPS)),)
+-include $(C_UPPER_DEPS)
+endif
+ifneq ($(strip $(CXX_DEPS)),)
+-include $(CXX_DEPS)
+endif
+ifneq ($(strip $(C_DEPS)),)
+-include $(C_DEPS)
+endif
+ifneq ($(strip $(CPP_DEPS)),)
+-include $(CPP_DEPS)
+endif
+endif
+
+-include ../makefile.defs
+
+# Add inputs and outputs from these tool invocations to the build variables
+
+# All Target
+all: Sniffles
+
+# Tool invocations
+Sniffles: $(OBJS) $(USER_OBJS)
+ @echo 'Building target: $@'
+ @echo 'Invoking: Cross G++ Linker'
+ g++ -L/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/lib -o "Sniffles" $(OBJS) $(USER_OBJS) $(LIBS)
+ @echo 'Finished building target: $@'
+ @echo ' '
+
+# Other Targets
+clean:
+ -$(RM) $(CC_DEPS)$(C++_DEPS)$(EXECUTABLES)$(OBJS)$(C_UPPER_DEPS)$(CXX_DEPS)$(C_DEPS)$(CPP_DEPS) Sniffles
+ - at echo ' '
+
+.PHONY: all clean dependents
+.SECONDARY:
+
+-include ../makefile.targets
diff --git a/cmake/objects.mk b/cmake/objects.mk
new file mode 100644
index 0000000..deb2291
--- /dev/null
+++ b/cmake/objects.mk
@@ -0,0 +1,8 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+USER_OBJS :=
+
+LIBS := -lbamtools
+
diff --git a/cmake/sources.mk b/cmake/sources.mk
new file mode 100644
index 0000000..2818222
--- /dev/null
+++ b/cmake/sources.mk
@@ -0,0 +1,30 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+C_UPPER_SRCS :=
+CXX_SRCS :=
+C++_SRCS :=
+OBJ_SRCS :=
+CC_SRCS :=
+ASM_SRCS :=
+C_SRCS :=
+CPP_SRCS :=
+O_SRCS :=
+S_UPPER_SRCS :=
+CC_DEPS :=
+C++_DEPS :=
+EXECUTABLES :=
+OBJS :=
+C_UPPER_DEPS :=
+CXX_DEPS :=
+C_DEPS :=
+CPP_DEPS :=
+
+# Every subdirectory with source files must be described here
+SUBDIRS := \
+src/tree \
+src/sub \
+src/plane-sweep \
+src \
+
diff --git a/cmake/src/Alignment.d b/cmake/src/Alignment.d
new file mode 100644
index 0000000..05f069c
--- /dev/null
+++ b/cmake/src/Alignment.d
@@ -0,0 +1,54 @@
+src/Alignment.d: ../src/Alignment.cpp ../src/Alignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/Paramer.h
+
+../src/Alignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/Paramer.h:
diff --git a/cmake/src/Alignment.o b/cmake/src/Alignment.o
new file mode 100644
index 0000000..c796031
Binary files /dev/null and b/cmake/src/Alignment.o differ
diff --git a/cmake/src/BamParser.d b/cmake/src/BamParser.d
new file mode 100644
index 0000000..f6c63d6
--- /dev/null
+++ b/cmake/src/BamParser.d
@@ -0,0 +1,58 @@
+src/BamParser.d: ../src/BamParser.cpp ../src/BamParser.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/Alignment.h ../src/Paramer.h ../src/Parser.h
+
+../src/BamParser.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/Alignment.h:
+
+../src/Paramer.h:
+
+../src/Parser.h:
diff --git a/cmake/src/BamParser.o b/cmake/src/BamParser.o
new file mode 100644
index 0000000..0721929
Binary files /dev/null and b/cmake/src/BamParser.o differ
diff --git a/cmake/src/Parser.d b/cmake/src/Parser.d
new file mode 100644
index 0000000..cc1ff22
--- /dev/null
+++ b/cmake/src/Parser.d
@@ -0,0 +1,56 @@
+src/Parser.d: ../src/Parser.cpp ../src/Parser.h ../src/Alignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/Paramer.h
+
+../src/Parser.h:
+
+../src/Alignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/Paramer.h:
diff --git a/cmake/src/Parser.o b/cmake/src/Parser.o
new file mode 100644
index 0000000..011d63e
Binary files /dev/null and b/cmake/src/Parser.o differ
diff --git a/cmake/src/Sniffles.d b/cmake/src/Sniffles.d
new file mode 100644
index 0000000..c2d7b20
--- /dev/null
+++ b/cmake/src/Sniffles.d
@@ -0,0 +1,148 @@
+src/Sniffles.d: ../src/Sniffles.cpp ../src/Paramer.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLine.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/SwitchArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Arg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ArgException.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Visitor.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLineInterface.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ArgTraits.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/StandardTraits.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/MultiSwitchArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/UnlabeledValueArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ValueArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Constraint.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/OptionalUnlabeledTracker.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/UnlabeledMultiArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/MultiArg.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/XorHandler.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/HelpVisitor.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLineOutput.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/VersionVisitor.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/IgnoreRestVisitor.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/StdOutput.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ValuesConstraint.h \
+ ../src/sub/Detect_Breakpoints.h ../src/sub/../BamParser.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/Alignment.h ../src/Parser.h \
+ ../src/sub/../plane-sweep/Plane-sweep.h \
+ ../src/sub/../plane-sweep/IContainer.h \
+ ../src/sub/../plane-sweep/Node.h ../src/sub/../plane-sweep/MyList.h \
+ ../src/sub/../plane-sweep/MyHeap.h ../src/sub/../tree/IntervallTree.h \
+ ../src/sub/../tree/TNode.h ../src/sub/../tree/../sub/Breakpoint.h
+
+../src/Paramer.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLine.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/SwitchArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Arg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ArgException.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Visitor.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLineInterface.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ArgTraits.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/StandardTraits.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/MultiSwitchArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/UnlabeledValueArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ValueArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/Constraint.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/OptionalUnlabeledTracker.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/UnlabeledMultiArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/MultiArg.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/XorHandler.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/HelpVisitor.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/CmdLineOutput.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/VersionVisitor.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/IgnoreRestVisitor.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/StdOutput.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include/tclap/ValuesConstraint.h:
+
+../src/sub/Detect_Breakpoints.h:
+
+../src/sub/../BamParser.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/Alignment.h:
+
+../src/Parser.h:
+
+../src/sub/../plane-sweep/Plane-sweep.h:
+
+../src/sub/../plane-sweep/IContainer.h:
+
+../src/sub/../plane-sweep/Node.h:
+
+../src/sub/../plane-sweep/MyList.h:
+
+../src/sub/../plane-sweep/MyHeap.h:
+
+../src/sub/../tree/IntervallTree.h:
+
+../src/sub/../tree/TNode.h:
+
+../src/sub/../tree/../sub/Breakpoint.h:
diff --git a/cmake/src/Sniffles.o b/cmake/src/Sniffles.o
new file mode 100644
index 0000000..44b4603
Binary files /dev/null and b/cmake/src/Sniffles.o differ
diff --git a/cmake/src/plane-sweep/Main.d b/cmake/src/plane-sweep/Main.d
new file mode 100644
index 0000000..ab8853b
--- /dev/null
+++ b/cmake/src/plane-sweep/Main.d
@@ -0,0 +1 @@
+src/plane-sweep/Main.d: ../src/plane-sweep/Main.cpp
diff --git a/cmake/src/plane-sweep/Main.o b/cmake/src/plane-sweep/Main.o
new file mode 100644
index 0000000..65ded99
Binary files /dev/null and b/cmake/src/plane-sweep/Main.o differ
diff --git a/cmake/src/plane-sweep/MyHeap.d b/cmake/src/plane-sweep/MyHeap.d
new file mode 100644
index 0000000..5cbba9f
--- /dev/null
+++ b/cmake/src/plane-sweep/MyHeap.d
@@ -0,0 +1,62 @@
+src/plane-sweep/MyHeap.d: ../src/plane-sweep/MyHeap.cpp \
+ ../src/plane-sweep/MyHeap.h ../src/plane-sweep/IContainer.h \
+ ../src/plane-sweep/Node.h ../src/plane-sweep/../Alignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/plane-sweep/../Paramer.h
+
+../src/plane-sweep/MyHeap.h:
+
+../src/plane-sweep/IContainer.h:
+
+../src/plane-sweep/Node.h:
+
+../src/plane-sweep/../Alignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/plane-sweep/../Paramer.h:
diff --git a/cmake/src/plane-sweep/MyHeap.o b/cmake/src/plane-sweep/MyHeap.o
new file mode 100644
index 0000000..01bff2a
Binary files /dev/null and b/cmake/src/plane-sweep/MyHeap.o differ
diff --git a/cmake/src/plane-sweep/MyList.d b/cmake/src/plane-sweep/MyList.d
new file mode 100644
index 0000000..e4df628
--- /dev/null
+++ b/cmake/src/plane-sweep/MyList.d
@@ -0,0 +1,62 @@
+src/plane-sweep/MyList.d: ../src/plane-sweep/MyList.cpp \
+ ../src/plane-sweep/MyList.h ../src/plane-sweep/Node.h \
+ ../src/plane-sweep/../Alignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/plane-sweep/../Paramer.h ../src/plane-sweep/IContainer.h
+
+../src/plane-sweep/MyList.h:
+
+../src/plane-sweep/Node.h:
+
+../src/plane-sweep/../Alignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/plane-sweep/../Paramer.h:
+
+../src/plane-sweep/IContainer.h:
diff --git a/cmake/src/plane-sweep/MyList.o b/cmake/src/plane-sweep/MyList.o
new file mode 100644
index 0000000..250cda9
Binary files /dev/null and b/cmake/src/plane-sweep/MyList.o differ
diff --git a/cmake/src/plane-sweep/PlaneSweep.d b/cmake/src/plane-sweep/PlaneSweep.d
new file mode 100644
index 0000000..4dc4f32
--- /dev/null
+++ b/cmake/src/plane-sweep/PlaneSweep.d
@@ -0,0 +1,67 @@
+src/plane-sweep/PlaneSweep.d: ../src/plane-sweep/PlaneSweep.cpp \
+ ../src/plane-sweep/Plane-sweep.h ../src/plane-sweep/IContainer.h \
+ ../src/plane-sweep/Node.h ../src/plane-sweep/../Alignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/plane-sweep/../Paramer.h ../src/plane-sweep/MyList.h \
+ ../src/plane-sweep/MyHeap.h
+
+../src/plane-sweep/Plane-sweep.h:
+
+../src/plane-sweep/IContainer.h:
+
+../src/plane-sweep/Node.h:
+
+../src/plane-sweep/../Alignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/plane-sweep/../Paramer.h:
+
+../src/plane-sweep/MyList.h:
+
+../src/plane-sweep/MyHeap.h:
diff --git a/cmake/src/plane-sweep/PlaneSweep.o b/cmake/src/plane-sweep/PlaneSweep.o
new file mode 100644
index 0000000..3214bfb
Binary files /dev/null and b/cmake/src/plane-sweep/PlaneSweep.o differ
diff --git a/cmake/src/plane-sweep/subdir.mk b/cmake/src/plane-sweep/subdir.mk
new file mode 100644
index 0000000..09def0d
--- /dev/null
+++ b/cmake/src/plane-sweep/subdir.mk
@@ -0,0 +1,33 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+# Add inputs and outputs from these tool invocations to the build variables
+CPP_SRCS += \
+../src/plane-sweep/Main.cpp \
+../src/plane-sweep/MyHeap.cpp \
+../src/plane-sweep/MyList.cpp \
+../src/plane-sweep/PlaneSweep.cpp
+
+OBJS += \
+./src/plane-sweep/Main.o \
+./src/plane-sweep/MyHeap.o \
+./src/plane-sweep/MyList.o \
+./src/plane-sweep/PlaneSweep.o
+
+CPP_DEPS += \
+./src/plane-sweep/Main.d \
+./src/plane-sweep/MyHeap.d \
+./src/plane-sweep/MyList.d \
+./src/plane-sweep/PlaneSweep.d
+
+
+# Each subdirectory must supply rules for building sources it contributes
+src/plane-sweep/%.o: ../src/plane-sweep/%.cpp
+ @echo 'Building file: $<'
+ @echo 'Invoking: Cross G++ Compiler'
+ g++ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include -O3 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<"
+ @echo 'Finished building: $<'
+ @echo ' '
+
+
diff --git a/cmake/src/sub/Detect_Breapoints.d b/cmake/src/sub/Detect_Breapoints.d
new file mode 100644
index 0000000..642c0fe
--- /dev/null
+++ b/cmake/src/sub/Detect_Breapoints.d
@@ -0,0 +1,82 @@
+src/sub/Detect_Breapoints.d: ../src/sub/Detect_Breapoints.cpp \
+ ../src/sub/Detect_Breakpoints.h ../src/sub/../BamParser.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/sub/../Alignment.h ../src/sub/../Paramer.h \
+ ../src/sub/../Parser.h ../src/sub/../plane-sweep/Plane-sweep.h \
+ ../src/sub/../plane-sweep/IContainer.h \
+ ../src/sub/../plane-sweep/Node.h ../src/sub/../plane-sweep/MyList.h \
+ ../src/sub/../plane-sweep/MyHeap.h ../src/sub/../tree/IntervallTree.h \
+ ../src/sub/../tree/TNode.h ../src/sub/../tree/../sub/Breakpoint.h
+
+../src/sub/Detect_Breakpoints.h:
+
+../src/sub/../BamParser.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/sub/../Alignment.h:
+
+../src/sub/../Paramer.h:
+
+../src/sub/../Parser.h:
+
+../src/sub/../plane-sweep/Plane-sweep.h:
+
+../src/sub/../plane-sweep/IContainer.h:
+
+../src/sub/../plane-sweep/Node.h:
+
+../src/sub/../plane-sweep/MyList.h:
+
+../src/sub/../plane-sweep/MyHeap.h:
+
+../src/sub/../tree/IntervallTree.h:
+
+../src/sub/../tree/TNode.h:
+
+../src/sub/../tree/../sub/Breakpoint.h:
diff --git a/cmake/src/sub/Detect_Breapoints.o b/cmake/src/sub/Detect_Breapoints.o
new file mode 100644
index 0000000..90caaa5
Binary files /dev/null and b/cmake/src/sub/Detect_Breapoints.o differ
diff --git a/cmake/src/sub/subdir.mk b/cmake/src/sub/subdir.mk
new file mode 100644
index 0000000..52e7d27
--- /dev/null
+++ b/cmake/src/sub/subdir.mk
@@ -0,0 +1,24 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+# Add inputs and outputs from these tool invocations to the build variables
+CPP_SRCS += \
+../src/sub/Detect_Breapoints.cpp
+
+OBJS += \
+./src/sub/Detect_Breapoints.o
+
+CPP_DEPS += \
+./src/sub/Detect_Breapoints.d
+
+
+# Each subdirectory must supply rules for building sources it contributes
+src/sub/%.o: ../src/sub/%.cpp
+ @echo 'Building file: $<'
+ @echo 'Invoking: Cross G++ Compiler'
+ g++ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include -O3 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<"
+ @echo 'Finished building: $<'
+ @echo ' '
+
+
diff --git a/cmake/src/subdir.mk b/cmake/src/subdir.mk
new file mode 100644
index 0000000..628d7e8
--- /dev/null
+++ b/cmake/src/subdir.mk
@@ -0,0 +1,33 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+# Add inputs and outputs from these tool invocations to the build variables
+CPP_SRCS += \
+../src/Alignment.cpp \
+../src/BamParser.cpp \
+../src/Parser.cpp \
+../src/Sniffles.cpp
+
+OBJS += \
+./src/Alignment.o \
+./src/BamParser.o \
+./src/Parser.o \
+./src/Sniffles.o
+
+CPP_DEPS += \
+./src/Alignment.d \
+./src/BamParser.d \
+./src/Parser.d \
+./src/Sniffles.d
+
+
+# Each subdirectory must supply rules for building sources it contributes
+src/%.o: ../src/%.cpp
+ @echo 'Building file: $<'
+ @echo 'Invoking: Cross G++ Compiler'
+ g++ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include -O3 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<"
+ @echo 'Finished building: $<'
+ @echo ' '
+
+
diff --git a/cmake/src/tree/IntervallTree.d b/cmake/src/tree/IntervallTree.d
new file mode 100644
index 0000000..81c1d79
--- /dev/null
+++ b/cmake/src/tree/IntervallTree.d
@@ -0,0 +1,67 @@
+src/tree/IntervallTree.d: ../src/tree/IntervallTree.cpp \
+ ../src/tree/IntervallTree.h ../src/tree/TNode.h \
+ ../src/tree/../sub/Breakpoint.h ../src/tree/../sub/../Paramer.h \
+ ../src/tree/../sub/../BamParser.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/tree/../sub/../Alignment.h ../src/tree/../sub/../Parser.h
+
+../src/tree/IntervallTree.h:
+
+../src/tree/TNode.h:
+
+../src/tree/../sub/Breakpoint.h:
+
+../src/tree/../sub/../Paramer.h:
+
+../src/tree/../sub/../BamParser.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/tree/../sub/../Alignment.h:
+
+../src/tree/../sub/../Parser.h:
diff --git a/cmake/src/tree/IntervallTree.o b/cmake/src/tree/IntervallTree.o
new file mode 100644
index 0000000..6d7955d
Binary files /dev/null and b/cmake/src/tree/IntervallTree.o differ
diff --git a/cmake/src/tree/Scapegoat.d b/cmake/src/tree/Scapegoat.d
new file mode 100644
index 0000000..aa6ba2d
--- /dev/null
+++ b/cmake/src/tree/Scapegoat.d
@@ -0,0 +1,66 @@
+src/tree/Scapegoat.d: ../src/tree/Scapegoat.cpp ../src/tree/Scapegoat.h \
+ ../src/tree/TNode.h ../src/tree/../sub/Breakpoint.h \
+ ../src/tree/../sub/../Paramer.h ../src/tree/../sub/../BamParser.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h \
+ /Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h \
+ ../src/tree/../sub/../Alignment.h ../src/tree/../sub/../Parser.h
+
+../src/tree/Scapegoat.h:
+
+../src/tree/TNode.h:
+
+../src/tree/../sub/Breakpoint.h:
+
+../src/tree/../sub/../Paramer.h:
+
+../src/tree/../sub/../BamParser.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamMultiReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/api_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/shared/bamtools_global.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamReader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAlignment.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamAux.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamConstants.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamIndex.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamHeader.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgramChain.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamProgram.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroupDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamReadGroup.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequenceDictionary.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/SamSequence.h:
+
+/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/api/BamWriter.h:
+
+../src/tree/../sub/../Alignment.h:
+
+../src/tree/../sub/../Parser.h:
diff --git a/cmake/src/tree/Scapegoat.o b/cmake/src/tree/Scapegoat.o
new file mode 100644
index 0000000..2d0980a
Binary files /dev/null and b/cmake/src/tree/Scapegoat.o differ
diff --git a/cmake/src/tree/subdir.mk b/cmake/src/tree/subdir.mk
new file mode 100644
index 0000000..9e331a8
--- /dev/null
+++ b/cmake/src/tree/subdir.mk
@@ -0,0 +1,27 @@
+################################################################################
+# Automatically-generated file. Do not edit!
+################################################################################
+
+# Add inputs and outputs from these tool invocations to the build variables
+CPP_SRCS += \
+../src/tree/IntervallTree.cpp \
+../src/tree/Scapegoat.cpp
+
+OBJS += \
+./src/tree/IntervallTree.o \
+./src/tree/Scapegoat.o
+
+CPP_DEPS += \
+./src/tree/IntervallTree.d \
+./src/tree/Scapegoat.d
+
+
+# Each subdirectory must supply rules for building sources it contributes
+src/tree/%.o: ../src/tree/%.cpp
+ @echo 'Building file: $<'
+ @echo 'Invoking: Cross G++ Compiler'
+ g++ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/bamtools/include/ -I/Users/fsedlaze/Documents/workspace/Sniffles/lib/tclap-1.2.1/include -O3 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<"
+ @echo 'Finished building: $<'
+ @echo ' '
+
+
diff --git a/src/Alignment.cpp b/src/Alignment.cpp
new file mode 100644
index 0000000..d7b28d2
--- /dev/null
+++ b/src/Alignment.cpp
@@ -0,0 +1,607 @@
+/*
+ * Alignments.cpp
+ *
+ * Created on: May 25, 2012
+ * Author: fritz
+ */
+
+#include "Alignment.h"
+
+void Alignment::setRef(string sequence) {
+ alignment.second = sequence;
+}
+void Alignment::initAlignment() {
+ al = new BamAlignment();
+}
+void Alignment::setAlignment(BamAlignment * align) {
+ al = align;
+ alignment.first.clear();
+ alignment.second.clear();
+ is_computed = false;
+
+ orig_length = al->QueryBases.size();
+ for (size_t i = 0; i < al->QueryBases.size(); i++) {
+ alignment.first += toupper(al->QueryBases[i]);
+ }
+ stop = this->getPosition() + this->getRefLength();
+}
+void Alignment::computeAlignment() {
+ int pos = 0;
+ for (size_t i = 0; i < al->CigarData.size(); i++) {
+ if (al->CigarData[i].Type == 'I') {
+ for (uint32_t t = 0; t < al->CigarData[i].Length; t++) {
+ alignment.second.insert(pos, "-");
+ alignment.second.erase(alignment.second.size() - 1, 1);
+ pos++;
+ }
+ } else if (al->CigarData[i].Type == 'D') {
+ for (uint32_t t = 0; t < al->CigarData[i].Length; t++) {
+ alignment.first.insert(pos, "-");
+ pos++;
+ }
+ } else if (al->CigarData[i].Type == 'S') {
+
+ if (pos == 0) { //front side
+ alignment.second.erase(((int) alignment.second.size()) - al->CigarData[i].Length, al->CigarData[i].Length);
+ } else { //backside
+ alignment.second.erase(pos, al->CigarData[i].Length);
+ }
+ alignment.first.erase(pos, al->CigarData[i].Length);
+
+ } else if (al->CigarData[i].Type == 'M') {
+ pos += al->CigarData[i].Length;
+ } else if (al->CigarData[i].Type == 'H') {
+
+ } else if (al->CigarData[i].Type == 'N') {
+ alignment.second.erase(pos, al->CigarData[i].Length);
+ }
+ }
+ for (size_t i = 0; i < alignment.first.size(); i++) {
+ if (alignment.first[i] == '=') {
+ alignment.first[i] = alignment.second[i];
+ }
+ }
+
+ is_computed = true;
+
+ if (alignment.first.size() != alignment.second.size()) {
+ cerr << "Error alignment has different length" << endl;
+ cerr << " ignoring alignment " << al->Name << endl;
+ cerr << al->Position << endl;
+
+ cerr << endl;
+ cerr << "read: " << alignment.first << endl;
+ cerr << endl;
+ cerr << " ref: " << alignment.second << endl;
+ cerr << endl;
+ cerr << orig_length << endl;
+ vector<CigarOp> cig = getCigar();
+
+ for (size_t i = 0; i < cig.size(); i++) {
+ cerr << cig[i].Length << cig[i].Type << " ";
+ }
+ cerr << endl;
+ exit(0);
+ return;
+ }
+}
+int32_t Alignment::getPosition() {
+ return al->Position;
+}
+int32_t Alignment::getRefID() {
+ return al->RefID;
+}
+bool Alignment::getStrand() {
+ return !al->IsReverseStrand();
+}
+vector<CigarOp> Alignment::getCigar() {
+ return al->CigarData;
+}
+string Alignment::getQualitValues() {
+ return al->Qualities;
+}
+size_t Alignment::get_length(std::vector<CigarOp> CigarData) {
+ size_t len = 0; //orig_length;
+ for (size_t i = 0; i < CigarData.size(); i++) {
+ if (CigarData[i].Type == 'D' || CigarData[i].Type == 'M' || CigarData[i].Type == 'N') {
+
+ len += CigarData[i].Length;
+ }
+ }
+ return len;
+}
+size_t Alignment::getRefLength() {
+ return get_length(this->al->CigarData);
+}
+size_t Alignment::getOrigLen() {
+ return orig_length;
+}
+pair<string, string> Alignment::getSequence() {
+ return alignment;
+}
+BamAlignment * Alignment::getAlignment() {
+ return al;
+}
+string Alignment::getName() {
+ return al->Name;
+}
+uint16_t Alignment::getMappingQual() {
+ return al->MapQuality;
+}
+float Alignment::getIdentity() {
+ if (is_computed) {
+ float match = 0;
+ for (size_t i = 0; i < alignment.first.size(); i++) {
+ if (alignment.first[i] == alignment.second[i]) {
+ match++;
+ }
+ }
+ return match / (float) alignment.first.size();
+ }
+ return -1;
+}
+int Alignment::getAlignmentFlag() {
+ return al->AlignmentFlag;
+}
+string Alignment::getQueryBases() {
+ return al->QueryBases;
+}
+string Alignment::getQualities() {
+ return al->Qualities;
+}
+string convertInt(int number) {
+ stringstream ss; //create a stringstream
+ ss << number; //add number to the stream
+ return ss.str(); //return a string with the contents of the stream
+}
+string Alignment::getTagData() {
+ vector<string> tags;
+
+ uint32_t i = 0;
+ if (al->GetTag("AS", i)) {
+ string tmp = "AS:i:";
+ tmp += convertInt(i);
+ tags.push_back(tmp);
+
+ }
+ i = 0;
+ if (al->GetTag("NM", i)) {
+ string tmp = "NM:i:";
+ tmp += convertInt(i);
+ tags.push_back(tmp);
+ }
+
+ string md;
+ if (al->GetTag("MD", md)) {
+ string tmp = "MD:Z:";
+ tmp += md;
+ tags.push_back(tmp);
+ }
+
+ i = 0;
+ if (al->GetTag("UQ", i)) {
+ string tmp = "UQ:i:";
+ tmp += convertInt(i);
+ tags.push_back(tmp);
+ }
+ string sa;
+ if (al->GetTag("SA", sa)) {
+ string tmp = "SA:Z:";
+ tmp += sa;
+ tags.push_back(tmp);
+ }
+
+ string res;
+ for (size_t i = 0; i < tags.size(); i++) {
+ res += tags[i];
+ if (i + 1 < tags.size()) {
+ res += '\t';
+ }
+ }
+ return res;
+}
+void Alignment::initSequence() {
+ this->alignment.first.clear();
+ this->alignment.second.clear();
+}
+
+int Alignment::get_id(RefVector ref, std::string chr) {
+ for (size_t i = 0; i < ref.size(); i++) {
+ if (strcmp(ref[i].RefName.c_str(), chr.c_str()) == 0) {
+ return i;
+ }
+ }
+ return -1; //should not happen!
+}
+void Alignment::get_coords(aln_str tmp, int & start, int &stop) {
+
+ if (tmp.cigar[0].Type == 'S' || tmp.cigar[0].Type == 'H') {
+ start = tmp.cigar[0].Length;
+ } else {
+ start = 0;
+ }
+
+ size_t index = tmp.cigar.size() - 1;
+ if (tmp.cigar[index].Type == 'S' || tmp.cigar[index].Type == 'H') {
+ stop = tmp.cigar[index].Length;
+ } else {
+ stop = 0;
+ }
+
+ if (!tmp.strand) {
+ int h = start;
+ start = stop;
+ stop = h;
+ }
+
+ /*start = 0;
+ stop = 0;
+ for (size_t i = 0; i < cigar.size(); i++) {
+ if (start == 0 && (cigar[i].Type == 'H' || cigar[i].Type == 'S')) {
+ start += cigar[i].Length;
+ stop += cigar[i].Length;
+ }
+ if (cigar[i].Type == 'I' || cigar[i].Type == 'M') {
+ stop += cigar[i].Length;
+ }
+ }*/
+}
+void sort_insert(aln_str tmp, vector<aln_str> &entries) {
+
+ for (vector<aln_str>::iterator i = entries.begin(); i != entries.end(); i++) {
+ if (tmp.read_pos_start < (*i).read_pos_start) {
+ entries.insert(i, tmp);
+ return;
+ }
+ }
+ entries.push_back(tmp);
+}
+vector<aln_str> Alignment::getSA(RefVector ref) {
+ string sa;
+
+ vector<aln_str> entries;
+ if (al->GetTag("SA", sa) && !sa.empty()) {
+
+ //store the main aln:
+ aln_str tmp;
+ tmp.RefID = this->getRefID();
+ tmp.cigar = this->getCigar();
+ tmp.length = this->getRefLength();
+ tmp.mq = this->getMappingQual();
+ tmp.pos = this->getPosition(); //+get_ref_lengths(tmp.RefID, ref);
+ tmp.strand = getStrand();
+ bool flag = strcmp(getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0;
+
+ get_coords(tmp, tmp.read_pos_start, tmp.read_pos_stop);
+ entries.push_back(tmp);
+ if (flag) {
+ std::cout << "Main Read: " << tmp.read_pos_start << " REF: " << tmp.pos << " " << tmp.RefID << std::endl;
+ }
+
+ //parse the rest:
+ size_t i = 0;
+ int count = 0;
+
+ std::string cigar;
+ std::string chr;
+ while (i < sa.size()) {
+ if (count == 0 && sa[i] != ',') {
+ chr += sa[i];
+ }
+ if (count == 1 && sa[i - 1] == ',') {
+ tmp.pos = atoi(&sa[i]);
+ }
+ if (count == 2 && sa[i - 1] == ',') {
+ tmp.strand = (bool) (sa[i] == '+');
+ }
+ if (count == 3 && sa[i] != ',') {
+ cigar += sa[i];
+ }
+ if (count == 4 && sa[i - 1] == ',') {
+ tmp.mq = atoi(&sa[i]);
+ }
+ if (count == 5 && sa[i] != ';') {
+ tmp.nm = atoi(&sa[i]);
+ }
+
+ if (sa[i] == ',') {
+ count++;
+ }
+ if (sa[i] == ';') {
+ if (tmp.mq > Parameter::Instance()->min_mq && entries.size() <= Parameter::Instance()->max_splits) {
+ //check this!
+ tmp.cigar = translate_cigar(cigar); //translates the cigar (string) to a type vector
+ get_coords(tmp, tmp.read_pos_start, tmp.read_pos_stop); //get the coords on the read.
+ tmp.length = get_length(tmp.cigar); //gives the length on the reference.
+ tmp.RefID = get_id(ref, chr); //translates back the chr to the id of the chr;
+ //TODO: should we do something about the MD string?
+ if (flag) {
+ std::cout << "Read: " << tmp.read_pos_start << " REF: " << tmp.pos << " " << tmp.RefID;
+ if (tmp.strand) {
+ std::cout << "+" << std::endl;
+ } else {
+ std::cout << "+" << std::endl;
+ }
+ }
+ //tmp.pos+=get_ref_lengths(tmp.RefID, ref);
+ //insert sorted:
+ includes_SV = true;
+ sort_insert(tmp, entries);
+ }
+ chr.clear();
+ cigar.clear();
+ tmp.cigar.clear();
+ count = 0;
+ tmp.mq = 0;
+ }
+ i++;
+ }
+ }
+ return entries;
+}
+
+//returns -1 if flags are not set!
+double Alignment::get_scrore_ratio() {
+ uint score = 0;
+ uint subscore = 0;
+ if (al->GetTag("AS", score) && al->GetTag("XS", subscore)) {
+ if(subscore==0){
+ return 40;// -1;
+ }
+ if (strcmp(getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0) {
+ std::cout << getName()<<" score: "<<(double) score / (double) subscore << std::endl;
+ }
+ return (double) score / (double) subscore;
+ }
+ return 100; //TODO: -1
+}
+bool Alignment::get_is_save() {
+ string sa;
+
+ double score = get_scrore_ratio();
+
+ /*if((al->GetTag("XA", sa) && !sa.empty())){
+ std::cout<<this->getName()<<"XA"<<std::endl;
+ }
+ if( (al->GetTag("XT", sa) && !sa.empty()) ){
+ std::cout<<this->getName()<<"XT"<<std::endl;
+ }*/
+
+ return !((al->GetTag("XA", sa) && !sa.empty()) || (al->GetTag("XT", sa) && !sa.empty()) || (score == -1 || score< Parameter::Instance()->score_treshold)); //TODO: 7.5
+}
+
+std::vector<CigarOp> Alignment::translate_cigar(std::string cigar) {
+
+ std::vector<CigarOp> new_cigar;
+
+ size_t i = 0;
+ bool first = true;
+ CigarOp tmp;
+ tmp.Length = -1;
+ while (i < cigar.size()) {
+ if (tmp.Length == -1) {
+ tmp.Length = atoi(&cigar[i]);
+ } else if (tmp.Length != -1 && atoi(&cigar[i]) == 0 && cigar[i] != '0') {
+ tmp.Type = cigar[i];
+ new_cigar.push_back(tmp);
+
+ tmp.Length = -1;
+ first = false;
+ }
+ i++;
+ }
+ return new_cigar;
+}
+
+double Alignment::get_avg_indel_length_Cigar() {
+ double len=0;
+ double num=0;
+ for (size_t i = 0; i < al->CigarData.size(); i++) {
+ if ((al->CigarData[i].Type == 'I'||al->CigarData[i].Type == 'D') && al->CigarData[i].Length > 1) {
+ len+= al->CigarData[i].Length;
+ num++;
+ }
+ }
+
+ return len/num;
+}
+
+vector<str_event> Alignment::get_events_CIGAR() {
+
+ size_t read_pos = 0;
+ size_t pos = this->getPosition(); //orig_length;
+ vector<str_event> events;
+ for (size_t i = 0; i < al->CigarData.size(); i++) {
+ if (al->CigarData[i].Type == 'H' || (al->CigarData[i].Type == 'S' || al->CigarData[i].Type == 'M')) {
+ read_pos += al->CigarData[i].Length;
+ }
+ if (al->CigarData[i].Type == 'D' && al->CigarData[i].Length > Parameter::Instance()->min_cigar_event) {
+ str_event ev;
+ ev.read_pos = read_pos;
+ ev.length = al->CigarData[i].Length; //deletion
+ ev.pos = pos;
+ includes_SV = true;
+ events.push_back(ev);
+ }
+ if (al->CigarData[i].Type == 'I' && al->CigarData[i].Length > Parameter::Instance()->min_cigar_event) {
+ // std::cout<<"CIGAR: "<<al->CigarData[i].Length<<" "<<this->getName()<<std::endl;
+ str_event ev;
+ ev.length = al->CigarData[i].Length * -1; //insertion;
+ ev.pos = pos;
+ ev.read_pos = read_pos;
+ includes_SV = true;
+ events.push_back(ev);
+ read_pos += al->CigarData[i].Length;
+ }
+ if (al->CigarData[i].Type == 'D' || al->CigarData[i].Type == 'M' || al->CigarData[i].Type == 'N') {
+ pos += al->CigarData[i].Length;
+ }
+ }
+
+ return events;
+}
+
+double Alignment::get_num_mismatches(std::string md) {
+ bool deletion = false;
+ bool match = false;
+ vector<int> helper;
+ double mis = 0;
+ double len = 0;
+ double maxim = 0;
+ for (size_t i = 0; i < md.size(); i += 20) {
+ mis = 0;
+ len = 0;
+ for (size_t j = 0; len < 100 && j + i < md.size(); j++) {
+ if (match && atoi(&md[i + j]) == 0 && md[i + j] != '0') { //is not a number:
+ if (md[i] == '^') {
+ deletion = true;
+ } else {
+ len++;
+ }
+ if (!deletion) {
+ //mistmatch!!
+ mis++;
+ match = false;
+ }
+ } else {
+ len += atoi(&md[i + j]);
+
+ match = true;
+ deletion = false;
+ }
+ }
+
+ if (strcmp(getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0) {
+ std::cout << (mis / len) << std::endl;
+ }
+ if ((mis / len) > maxim) {
+ maxim = (mis / len);
+ }
+ }
+ return maxim; // 0.03);
+}
+std::string Alignment::get_md() {
+ std::string md;
+ if (al->GetTag("MD", md)) {
+ return md;
+ }
+ return md;
+}
+vector<str_event> Alignment::get_events_MD(int min_mis) {
+ vector<str_event> events;
+ std::string md;
+ if (al->GetTag("MD", md)) {
+ //TODO: remove:
+ bool flag = strcmp(getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0;
+
+ if (flag) {
+ std::cout << "found!" << std::endl;
+ }
+ //TODO think of a good threshold!
+ if (get_num_mismatches(md) > Parameter::Instance()->min_num_mismatches) {
+ if (flag) {
+ std::cout << "is_strange!" << std::endl;
+ }
+ //generate a vector that holds the positions of the read
+ std::vector<int> aln;
+ int pos = getPosition();
+
+ for (size_t i = 0; i < al->CigarData.size(); i++) {
+ if (al->CigarData[i].Type == 'I') { //TODO check
+ }
+ if (al->CigarData[i].Type == 'D') {
+ pos += al->CigarData[i].Length;
+ }
+ if (al->CigarData[i].Type == 'M') {
+ for (size_t j = 0; j < al->CigarData[i].Length; j++) {
+ aln.push_back(pos);
+ pos++;
+ //aln += "=";
+ }
+ }
+ }
+ //fill in the mismatches:
+ bool deletion = false;
+ bool match = false;
+ double mis = 0;
+ double len = 0;
+ //MD:Z:106G48^C33C41T0G5^G15G0C52^T15^C5^GC16G0A48^T3^C9^G3^G17^G23^G3^C2^G40C0G6
+ for (size_t i = 0; i < md.size(); i++) {
+ if ((atoi(&md[i]) == 0 && md[i] != '0')) { //is not a number:
+ if (md[i] == '^') {
+ deletion = true;
+ }
+ if (!deletion) {
+ //mistmatch!!
+ mis++;
+ aln[len] = aln[len] * -1;
+ len++;
+ }
+ match = false;
+ } else if (!match) {
+ len += atoi(&md[i]);
+ match = true;
+ deletion = false;
+ }
+ }
+
+ /* if (flag) {
+ for (size_t i = 0; i < aln.size(); i++) {
+ std::cout << aln[i] << " ";
+ }
+ std::cout << endl;
+ }
+*/
+ int runlength = 100;
+ str_event ev;
+ ev.pos = -1;
+ ev.length = -1;
+ ev.read_pos = 0;
+ int start = 0;
+ int last = 0;
+ for (size_t i = 0; i < aln.size(); i += 50) { //+=runlength/2 ??
+ //std::cout<<aln[i]<<";";
+ int mis = 0;
+ int first = 0;
+
+ for (size_t j = 0; (j + i) < aln.size() && j < runlength; j++) {
+ if (aln[(i + j)] < 0) {
+ if (first == 0) {
+ first = abs(aln[(i + j)]);
+ }
+ last = abs(aln[(i + j)]);
+ mis++;
+ }
+ }
+ if (mis > min_mis) { //TOOD ratio?
+ if (ev.pos == -1) {
+ start = i;
+ ev.pos = first;
+ ev.read_pos = ev.pos - getPosition();
+ }
+ } else {
+ if ((start > 20 && abs((int) (i + runlength) - (int) aln.size()) > 20) && ev.pos != -1) {
+ if (flag) {
+ std::cout << i << " " << (i + runlength) << " " << aln.size() << std::endl;
+ std::cout << ev.pos << " " << last << " " << std::endl;
+ }
+ includes_SV = true;
+ ev.length = last - ev.pos;
+ if (flag) {
+ std::cout << ev.pos << " " << ev.length << std::endl;
+ }
+ if (ev.length > runlength) {
+ events.push_back(ev);
+ }
+ last = 0;
+ ev.pos = -1;
+ } else {
+ ev.pos = -1;
+ }
+ }
+ }
+ }
+
+ }
+ return events;
+}
diff --git a/src/Alignment.h b/src/Alignment.h
new file mode 100644
index 0000000..1894519
--- /dev/null
+++ b/src/Alignment.h
@@ -0,0 +1,110 @@
+/*
+ * Alignments.h
+ *
+ * Created on: May 25, 2012
+ * Author: fritz
+ */
+
+#ifndef ALIGNMENTS_H_
+#define ALIGNMENTS_H_
+
+#include <string.h>
+#include <vector>
+#include <math.h>
+#include "api/BamAux.h"
+#include "api/BamMultiReader.h"
+#include "api/BamWriter.h"
+#include "Paramer.h"
+using namespace BamTools;
+using namespace std;
+
+typedef unsigned short ushort;
+typedef unsigned int uint;
+
+struct str_event{
+ short length;
+ int pos;
+ int read_pos;
+};
+struct aln_str{
+ int RefID;
+ int pos;
+ bool strand;
+ std::vector<CigarOp> cigar;
+ ushort mq;
+ ushort nm;
+ ushort length;
+ int read_pos_start;
+ int read_pos_stop;
+};
+
+class Alignment {
+
+private:
+ BamAlignment * al;
+ bool includes_SV;
+ pair<string,string> alignment;
+ bool is_computed;
+ int32_t orig_length;
+ int stop;
+ std::vector<CigarOp> translate_cigar(std::string cigar);
+ size_t get_length(std::vector<CigarOp> CigarData);
+
+ int get_id(RefVector ref, std::string chr);
+public:
+ Alignment(){
+ stop=0;
+ orig_length=0;
+ al=NULL;
+ is_computed=false;
+ includes_SV=false;
+ }
+ ~Alignment(){
+ alignment.first.clear();
+ alignment.second.clear();
+ delete al;
+ }
+ void setAlignment(BamAlignment * al);
+ void setRef(string sequence);
+ void computeAlignment();
+
+ pair<string,string> getSequence();
+ int32_t getPosition();
+ int32_t getRefID();
+ bool getStrand();
+ uint16_t getMappingQual();
+ string getName();
+ vector<CigarOp> getCigar();
+ string getQualitValues();
+ size_t getRefLength();
+ size_t getOrigLen();
+ BamAlignment * getAlignment();
+ float getIdentity();
+ void initAlignment();
+ int getAlignmentFlag();
+ string getQueryBases();
+ string getQualities();
+ string getTagData();
+ vector<aln_str> getSA(RefVector ref);
+ void initSequence();
+ int get_stop(){
+ return stop;
+ }
+ vector<str_event> get_events_CIGAR();
+ vector<str_event> get_events_MD(int min_length);
+ void get_coords(aln_str tmp,int & start,int &stop);
+ bool supports_SV(){
+ return this->includes_SV;
+ }
+ void set_supports_SV(bool flag){
+ this->includes_SV=flag;
+ }
+ bool get_is_save();
+ double get_num_mismatches(std::string md);
+ double get_scrore_ratio();
+ std::string get_md();
+ double get_avg_indel_length_Cigar();
+};
+
+
+#endif /* ALIGNMENTS_H_ */
diff --git a/src/BamParser.cpp b/src/BamParser.cpp
new file mode 100644
index 0000000..83fb9b3
--- /dev/null
+++ b/src/BamParser.cpp
@@ -0,0 +1,56 @@
+/*
+ * parser.cpp
+ *
+ * Created on: Apr 17, 2012
+ * Author: fritz
+ */
+
+#include "BamParser.h"
+
+BamParser::BamParser(string file){
+ vector<string > tmps;
+ tmps.push_back(file);
+
+ if(!reader.Open(tmps)){
+ cerr<<"BAM Parser: could not open file: "<<file<<endl;
+ exit(0);
+ }
+
+}
+
+Alignment* BamParser::parseRead(uint16_t mappingQv){
+
+ Alignment *align = new Alignment();
+ BamAlignment* al = new BamAlignment();
+ while(reader.GetNextAlignmentCore(al[0])){
+ if( al->IsMapped() && al->MapQuality > mappingQv){
+ al->BuildCharData();
+ align->setAlignment(al);
+ return align;
+ }
+ }
+ return align;
+
+}
+void BamParser::parseReadFast(uint16_t mappingQv,Alignment*& align){
+
+// Alignment *align = new Alignment();
+ BamAlignment* al = align->getAlignment();
+// getSequence().first
+ align->initSequence();
+
+ while(reader.GetNextAlignmentCore(al[0])){
+ if( al->IsMapped() && al->MapQuality > mappingQv){
+ al->BuildCharData();
+ align->setAlignment(al);
+ return;
+ }
+ }
+}
+RefVector BamParser::get_refInfo(){
+ return reader.GetReferenceData();
+}
+
+string BamParser::get_header(){
+ return reader.GetHeaderText();
+}
diff --git a/src/BamParser.h b/src/BamParser.h
new file mode 100644
index 0000000..6344ddc
--- /dev/null
+++ b/src/BamParser.h
@@ -0,0 +1,42 @@
+/*
+ * parser.h
+ *
+ * Created on: Apr 17, 2012
+ * Author: fritz
+ */
+
+#ifndef BAMPARSER_H_
+#define BAMPARSER_H_
+
+#include "api/BamMultiReader.h"
+#include "api/BamWriter.h"
+#include "Alignment.h"
+#include "Parser.h"
+
+#include <iostream>
+#include <fstream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+
+
+using namespace BamTools;
+using namespace std;
+
+class BamParser: public Parser {
+private:
+ BamMultiReader reader;
+
+public:
+ BamParser(string file);
+ ~BamParser(){
+ reader.Close();
+ }
+ Alignment * parseRead(uint16_t mappingQv);
+ void parseReadFast(uint16_t mappingQv,Alignment*& aln);
+ string get_header();
+ RefVector get_refInfo();
+};
+
+#endif /* PARSER_H_ */
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000..4a524c2
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,60 @@
+cmake_minimum_required(VERSION 2.8)
+project(Sniffles)
+
+include_directories (../lib/bamtools-2.3.0/src)
+include_directories(../lib/tclap-1.2.1/include)
+
+configure_file( Version.h.in ${CMAKE_SOURCE_DIR}/src/Version.h )
+
+add_executable(sniffles
+ Alignment.cpp
+ BamParser.cpp
+ Parser.cpp
+ Sniffles.cpp
+ Ignore_Regions.cpp
+ tree/Intervall_bed.cpp
+ sub/Detect_Breakpoints.cpp
+ sub/Breakpoint.cpp
+ tree/IntervallTree.cpp
+ realign/SWCPU.cpp
+ realign/Realign.cpp
+ print/VCFPrinter.cpp
+ print/BedePrinter.cpp
+ print/MariaPrinter.cpp
+ print/IPrinter.cpp
+ phasing/PhaserSV.cpp
+ tree/BinTree.cpp
+ print/NGMPrinter.cpp
+ )
+
+#target_link_libraries(ngm-core pthread)
+TARGET_LINK_LIBRARIES(sniffles BamTools-static)
+TARGET_LINK_LIBRARIES(sniffles zlibstatic)
+
+add_executable(sniffles-debug
+ Alignment.cpp
+ BamParser.cpp
+ Parser.cpp
+ Sniffles.cpp
+ Ignore_Regions.cpp
+ tree/Intervall_bed.cpp
+ sub/Detect_Breakpoints.cpp
+ sub/Breakpoint.cpp
+ tree/IntervallTree.cpp
+ realign/SWCPU.cpp
+ realign/Realign.cpp
+ print/VCFPrinter.cpp
+ print/BedePrinter.cpp
+ print/MariaPrinter.cpp
+ print/IPrinter.cpp
+ phasing/PhaserSV.cpp
+ tree/BinTree.cpp
+ print/NGMPrinter.cpp
+ )
+
+SET_TARGET_PROPERTIES(sniffles-debug PROPERTIES COMPILE_FLAGS "-g3 -O0")
+
+#target_link_libraries(sniffles-debug pthread)
+TARGET_LINK_LIBRARIES(sniffles-debug BamTools-static)
+TARGET_LINK_LIBRARIES(sniffles-debug zlibstatic)
+
diff --git a/src/Ignore_Regions.cpp b/src/Ignore_Regions.cpp
new file mode 100644
index 0000000..bac0881
--- /dev/null
+++ b/src/Ignore_Regions.cpp
@@ -0,0 +1,112 @@
+/*
+ * Ignore_Regions.cpp
+ *
+ * Created on: Feb 4, 2016
+ * Author: fsedlaze
+ */
+
+#include "Ignore_Regions.h"
+#include "sub/Detect_Breakpoints.h"
+
+long get_ref_coords(std::string chr, RefVector ref) {
+ long length = 0;
+ for (size_t i = 0; i < ref.size(); i++) {
+ if (strcmp(ref[i].RefName.c_str(), chr.c_str()) == 0) {
+ return length;
+ }
+ length += ref[i].RefLength + Parameter::Instance()->max_dist;
+ }
+ return -1; //should not happen
+
+}
+
+
+long get_ref_lengths2(int id, RefVector ref) {
+ long length = 0;
+
+ for (size_t i = 0; i < (size_t) id && i < ref.size(); i++) {
+ length += ref[i].RefLength + Parameter::Instance()->max_dist;
+ }
+ return length;
+}
+
+int get_id2(RefVector ref, std::string chr) {
+ for (size_t i = 0; i < ref.size(); i++) {
+ if (strcmp(ref[i].RefName.c_str(), chr.c_str()) == 0) {
+ return i;
+ }
+ }
+ return -1; //should not happen!
+}
+
+
+void initialize_bed(IntervallTree_bed &bed_tree, Leaf *&root,RefVector ref) {
+ //bst.insert(point, root);
+ size_t buffer_size = 2000000;
+ char*buffer = new char[buffer_size];
+ std::ifstream myfile;
+
+ myfile.open(Parameter::Instance()->ignore_regions_bed.c_str(), std::ifstream::in);
+ if (!myfile.good()) {
+ std::cout << "SAM Parser: could not open file: " << Parameter::Instance()->ignore_regions_bed.c_str() << std::endl;
+ exit(0);
+ }
+
+ myfile.getline(buffer, buffer_size);
+
+ while (!myfile.eof()) {
+ int count = 0;
+ string chr;
+ int p1;
+ int p2;
+ for (size_t i = 0; i < buffer_size && buffer[i] != '\0' && buffer[i] != '\n'; i++) {
+ if (count == 0 && buffer[i] != '\t') {
+ chr += buffer[i];
+ }
+ if (count == 1 && buffer[i - 1] == '\t') {
+ p1 = atoi(&buffer[i]);
+ }
+ if (count == 2 && buffer[i - 1] == '\t') {
+ p2 = atoi(&buffer[i]);
+ break;
+ }
+ if (buffer[i] == '\t') {
+ count++;
+ }
+ }
+ //transfer coordinates:
+ long ref_dist = get_ref_coords(chr, ref);
+ std::cout << (long) p1 + ref_dist << " " << (long) p2 + ref_dist << std::endl;
+ bed_tree.insert((long) p1 + ref_dist, (long) p2 + ref_dist, root);
+
+ myfile.getline(buffer, buffer_size);
+ }
+
+}
+
+void ignore_regions(std::vector<Breakpoint *> & final_SV,RefVector ref) {
+
+ IntervallTree_bed bed_tree;
+ Leaf *root = NULL;
+ initialize_bed(bed_tree, root,ref);
+
+ bed_tree.postorder(root);
+ size_t i = 0;
+ while (i < final_SV.size()) {
+ if (final_SV[i]->get_SVtype() & DUP) {
+ std::cout << final_SV[i]->get_coordinates().start.most_support << " ";
+ }
+ if (bed_tree.is_in(final_SV[i]->get_coordinates().start.most_support, root) || bed_tree.is_in(final_SV[i]->get_coordinates().stop.most_support, root)) {
+ if (final_SV[i]->get_SVtype() & DUP) {
+ std::cout << "erase" << endl;
+ }
+ final_SV.erase(final_SV.begin() + i);
+ } else {
+ if (final_SV[i]->get_SVtype() & DUP) {
+ std::cout << "keep" << endl;
+ }
+ i++;
+ }
+ }
+}
+
diff --git a/src/Ignore_Regions.h b/src/Ignore_Regions.h
new file mode 100644
index 0000000..414a139
--- /dev/null
+++ b/src/Ignore_Regions.h
@@ -0,0 +1,15 @@
+/*
+ * Ignore_Regions.h
+ *
+ * Created on: Feb 4, 2016
+ * Author: fsedlaze
+ */
+
+#ifndef IGNORE_REGIONS_H_
+#define IGNORE_REGIONS_H_
+#include "sub/Breakpoint.h"
+#include "tree/Intervall_bed.h"
+
+void ignore_regions(std::vector<Breakpoint *> & final_SV);
+void initialize_bed(IntervallTree_bed &bed_tree, Leaf *&root,RefVector ref);
+#endif /* IGNORE_REGIONS_H_ */
diff --git a/src/Paramer.h b/src/Paramer.h
new file mode 100644
index 0000000..6c4e4f4
--- /dev/null
+++ b/src/Paramer.h
@@ -0,0 +1,93 @@
+/*
+ * Paramer.h
+ *
+ * Created on: Aug 20, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PARAMER_H_
+#define PARAMER_H_
+
+#include <string.h>
+#include <string>
+#include <vector>
+#include <stdlib.h>
+#include <iostream>
+struct region_str {
+ std::string chr;
+ int start;
+ int stop;
+};
+
+class Parameter {
+private:
+ Parameter() {
+ score_treshold = 0;
+ min_num_mismatches = 0.3;
+ }
+ ~Parameter() {
+
+ }
+ static Parameter* m_pInstance;
+ std::vector<region_str> regions;
+public:
+ std::string output_vcf;
+ std::string output_bede;
+ std::string ref_seq;
+ std::string read_name;
+ std::string output_maria;
+ std::string ignore_regions_bed;
+
+ std::vector<std::string> bam_files;
+ short min_mq;
+ short min_cigar_event;
+ short report_n_reads;
+ short corridor;
+
+ double error_rate;
+ double score_treshold;
+ double min_num_mismatches;
+
+ int min_support;
+ int max_splits;
+ int max_dist;
+ int min_length;
+ int min_reads_phase;
+ int num_threads;
+
+ bool realign;
+ bool useMD_CIGAR;
+
+ void set_regions(std::string reg) {
+ size_t i = 0;
+ region_str tmp;
+ short sep;
+ while (i < reg.size()) {
+ tmp.chr = reg.substr(i, reg.find_first_of(':'));
+ i += tmp.chr.size() + 1;
+ tmp.start = atoi(®[i]);
+ i += reg.find_first_of('-') + 1;
+ tmp.stop = atoi(®[i]);
+ i += reg.find_first_of(';') + 1;
+ regions.push_back(tmp);
+ }
+ std::cout << "found regions: " << regions.size() << std::endl;
+ }
+ bool overlaps(std::string chr, int start, int stop) {
+ for (size_t i = 0; i < regions.size(); i++) {
+ if (strcmp(chr.c_str(), regions[i].chr.c_str()) == 0 && (abs(start - regions[i].start) < max_dist && abs(stop - regions[i].stop) < max_dist)) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ static Parameter* Instance() {
+ if (!m_pInstance) {
+ m_pInstance = new Parameter;
+ }
+ return m_pInstance;
+ }
+};
+
+#endif /* PARAMER_H_ */
diff --git a/src/Parser.cpp b/src/Parser.cpp
new file mode 100644
index 0000000..82c373b
--- /dev/null
+++ b/src/Parser.cpp
@@ -0,0 +1,8 @@
+/*
+ * Parser.cpp
+ *
+ * Created on: May 29, 2012
+ * Author: fritz
+ */
+
+#include "Parser.h"
diff --git a/src/Parser.h b/src/Parser.h
new file mode 100644
index 0000000..a15ff42
--- /dev/null
+++ b/src/Parser.h
@@ -0,0 +1,21 @@
+/*
+ * Parser.h
+ *
+ * Created on: May 29, 2012
+ * Author: fritz
+ */
+
+#ifndef PARSER_H_
+#define PARSER_H_
+
+#include "Alignment.h"
+
+class Parser {
+
+public:
+ virtual ~Parser(){};
+ virtual Alignment * parseRead(uint16_t mappingQv) = 0;
+ virtual void parseReadFast(uint16_t mappingQv,Alignment *& align) = 0;
+};
+
+#endif /* PARSER_H_ */
diff --git a/src/Sniffles.cpp b/src/Sniffles.cpp
new file mode 100644
index 0000000..e14f25d
--- /dev/null
+++ b/src/Sniffles.cpp
@@ -0,0 +1,178 @@
+//============================================================================
+// Name : Sniffles.cpp
+// Author : Fritz Sedlazeck
+// Version :
+// Copyright : Your copyright notice
+// Description : Hello World in C++, Ansi-style
+//============================================================================
+// phil: cd ~/hetero/philipp/pacbio/example-svs/reads
+//For mac: cmake -D CMAKE_C_COMPILER=/opt/local/bin/gcc-mp-4.7 -D CMAKE_CXX_COMPILER=/opt/local/bin/g++-mp-4.7 ..
+#include <iostream>
+#include "Paramer.h"
+#include <tclap/CmdLine.h>
+#include <omp.h>
+#include "realign/Realign.h"
+#include "sub/Detect_Breakpoints.h"
+#include "print/IPrinter.h"
+#include "print/BedePrinter.h"
+#include "print/VCFPrinter.h"
+#include "print/MariaPrinter.h"
+#include "phasing/PhaserSV.h"
+#include "print/NGMPrinter.h"
+#include "Ignore_Regions.h"
+
+Parameter* Parameter::m_pInstance = NULL;
+//cmake -D CMAKE_C_COMPILER=/opt/local/bin/gcc-mp-4.7 -D CMAKE_CXX_COMPILER=/opt/local/bin/g++-mp-4.7 ..
+
+//TODO: Friday: redo the positioning. Min,Max as representation, but collect the position and give the most likely!
+
+// TODO: Finish virtual Region class, Inherit Breakpoint and strange regions; Implement the region into the intervall tree.
+
+// Think of method to filter out strange SV.
+// flushing the interval tree after a certain amount of SV/ bp? Would speed up things.
+// Allow for Illumina split read data?
+// Think again of computing the coverage?
+// VCF and bede output
+// Think of multiple bam files -> setting genotypes
+// Think about overlapping SV, maybe flag to report if they share the same read -> phasing info?
+// Regular scan through the SV and move those where the end point lies far behind the current pos or reads. Eg. 1MB?
+
+//TODO: look into the gene fusions.
+
+void read_parameters(int argc, char *argv[]) {
+ TCLAP::CmdLine cmd("Sniffles version 0.0.1", ' ', "0.0.1");
+
+ TCLAP::ValueArg<std::string> arg_bamfile("m", "mapped_reads", "Bam File", true, "", "string");
+ TCLAP::ValueArg<std::string> arg_vcf("v", "vcf", "VCF output file name", false, "", "string");
+ //TCLAP::ValueArg<std::string> arg_bede("b", "bede", "Bede output file name", false, "", "string");
+ TCLAP::ValueArg<std::string> arg_maria("b", "bede", "Simplified format of bede Format.", false, "", "string");
+ //TCLAP::ValueArg<std::string> arg_noregions("", "bed", "Ignore SV overlapping with those regions.", false, "", "string");
+// TCLAP::ValueArg<std::string> arg_ref("r", "reference", "Reference fasta sequence. Activates realignment step", false, "", "string");
+ //TCLAP::ValueArg<std::string> arg_region("", "regions", "List of regions CHR:start-stop; to check", false, "", "string");
+ TCLAP::ValueArg<int> arg_support("s", "min_support", "Minimum number of reads that support a SV. Default: 10", false, 10, "int");
+ TCLAP::ValueArg<int> arg_splits("", "max_num_splits", "Maximum number of splits per read to be still taken into account. Default: 4", false, 4, "int");
+ TCLAP::ValueArg<int> arg_dist("d", "max_distance", "Maximum distance to group SV together. Default: 1kb", false, 1000, "int");
+ TCLAP::ValueArg<int> arg_threads("t", "threads", "Number of threads to use. Default: 3", false, 3, "int");
+ //TCLAP::ValueArg<int> arg_corridor("", "corridor", "Maximum size of corridor for realignment. Default: 2000", false, 2000, "int");
+ TCLAP::ValueArg<int> arg_minlength("l", "min_length", "Minimum length of SV to be reported. Default:0", false, 0, "int");
+ TCLAP::ValueArg<int> arg_mq("q", "minmapping_qual", "Minimum Mapping Quality. Default: 20", false, 20, "int");
+ TCLAP::ValueArg<int> arg_cigar("c", "min_cigar_event", "Minimum Cigar Event (e.g. Insertion, deletion) to take into account. Default:50 ", false, 50, "int");
+ TCLAP::ValueArg<int> arg_numreads("n", "num_reads_report", "Report up to N reads that support the SV. Default: 0", false, 0, "int");
+// TCLAP::ValueArg<int> arg_phase_minreads("", "min_reads_phase", "Minimum reads overlapping two SV to phase them together. Default: 1", false, 1, "int");
+
+// TCLAP::SwitchArg arg_realign("", "re-align", "Enables the realignment of reads at predicted SV sites. Leads to more accurate breakpoint predictions.", cmd, false);
+ TCLAP::SwitchArg arg_MD_cigar("", "use_MD_Cigar", "Enables Sniffles to use the alignment information to screen for suspicious regions.", cmd, false);
+
+ cmd.add(arg_numreads);
+ cmd.add(arg_dist);
+ //cmd.add(arg_bede);
+// cmd.add(arg_noregions);
+ cmd.add(arg_threads);
+ cmd.add(arg_cigar);
+ cmd.add(arg_maria);
+ cmd.add(arg_vcf);
+// cmd.add(arg_ref);
+ //cmd.add(arg_corridor);
+ //cmd.add(arg_region);
+ cmd.add(arg_minlength);
+ //cmd.add(arg_phase_minreads);
+ cmd.add(arg_mq);
+ cmd.add(arg_splits);
+ cmd.add(arg_support);
+ cmd.add(arg_bamfile);
+ //parse cmd:
+ cmd.parse(argc, argv);
+
+ Parameter::Instance()->read_name = " "; //m150118_093731_00118_c100767352550000001823169407221590_s1_p0/27613/0_15682" ;//TODO: just for debuging reasons!
+ Parameter::Instance()->bam_files.push_back(arg_bamfile.getValue());
+ Parameter::Instance()->min_mq = arg_mq.getValue();
+ Parameter::Instance()->output_vcf = arg_vcf.getValue();
+ Parameter::Instance()->min_cigar_event = arg_cigar.getValue();
+ Parameter::Instance()->report_n_reads = arg_numreads.getValue();
+ Parameter::Instance()->min_support = arg_support.getValue();
+ Parameter::Instance()->max_splits = arg_splits.getValue();
+ Parameter::Instance()->max_dist = arg_dist.getValue();
+ //Parameter::Instance()->ref_seq = arg_ref.getValue();
+ Parameter::Instance()->realign = false; //arg_realign.getValue();
+ //Parameter::Instance()->corridor = arg_corridor.getValue();
+// Parameter::Instance()->output_bede = arg_bede.getValue();
+ Parameter::Instance()->min_length = arg_minlength.getValue();
+ //Parameter::Instance()->min_reads_phase = arg_phase_minreads.getValue();
+ Parameter::Instance()->useMD_CIGAR = arg_MD_cigar.getValue();
+ Parameter::Instance()->num_threads = arg_threads.getValue();
+ Parameter::Instance()->output_maria = arg_maria.getValue();
+ //Parameter::Instance()->ignore_regions_bed = arg_noregions.getValue();
+}
+
+int main(int argc, char *argv[]) {
+
+ try {
+ //init parameter and reads user defined parameter from command line.
+ read_parameters(argc, argv);
+
+ //init openmp:
+ omp_set_dynamic(0);
+ omp_set_num_threads(Parameter::Instance()->num_threads);
+
+ //init printer:
+ IPrinter * printer;
+ if (!Parameter::Instance()->ref_seq.empty()) {
+ printer = new NGMPrinter();
+ } else if (!Parameter::Instance()->output_vcf.empty()) {
+ printer = new VCFPrinter();
+ std::cout << "VCF parser" << std::endl;
+ } else if (!Parameter::Instance()->output_bede.empty()) {
+ printer = new BedePrinter();
+ } else if (!Parameter::Instance()->output_maria.empty()) {
+ printer = new MariaPrinter();
+ } else {
+ std::cerr << "Please specify an output file using -v or -b" << std::endl;
+ return -1;
+ }
+
+ printer->init();
+
+ //TODO add support of multiple files!
+ for (size_t i = 0; i < Parameter::Instance()->bam_files.size(); i++) {
+ detect_breakpoints(Parameter::Instance()->bam_files[i], printer);
+ }
+
+ //realignment: Using NGM???
+ if (!Parameter::Instance()->ref_seq.empty()) {
+ std::cout<<"Realignment step activated:"<<std::endl;
+ //1. parse in output file(s)
+ //2. extract reads from bam files (read groups?) //shellscript: picard/samtools??
+ //3. realign
+ //4. call Sniffles
+ }
+
+ //grouping SV together:
+ //PhaserSV * phaser= new PhaserSV();
+ //phaser->phase(final_SV);
+ //delete phaser;
+
+ //printing results:
+
+ } catch (TCLAP::ArgException &e) // catch any exceptions
+ {
+ std::cerr << "error: " << e.error() << " for arg " << e.argId() << std::endl;
+ }
+ return 0;
+}
+/*
+ * Input:
+ * Regions to check
+ *
+ * Parameters
+ *
+ *
+ */
+
+/*
+ * 1. Detect strange regions
+ * Using MD, Cigar, Split reads
+ *
+ * 2. Extract reads from region (?) The main read holds the whole sequence for BWA-MEM
+ *
+ * 3. Realign regions (?)
+ */
diff --git a/src/Version.h b/src/Version.h
new file mode 100644
index 0000000..09ef3c8
--- /dev/null
+++ b/src/Version.h
@@ -0,0 +1,9 @@
+#ifndef VERSION_H
+#define VERSION_H
+
+#define VERSION_MAJOR ""
+#define VERSION_MINOR ""
+#define VERSION_BUILD ""
+
+#endif // VERSION_H
+
diff --git a/src/Version.h.in b/src/Version.h.in
new file mode 100644
index 0000000..8e7f2bd
--- /dev/null
+++ b/src/Version.h.in
@@ -0,0 +1,9 @@
+#ifndef VERSION_H
+#define VERSION_H
+
+#define VERSION_MAJOR "@NGM_VERSION_MAJOR@"
+#define VERSION_MINOR "@NGM_VERSION_MINOR@"
+#define VERSION_BUILD "@NGM_VERSION_BUILD@"
+
+#endif // VERSION_H
+
diff --git a/src/phasing/PhaserSV.cpp b/src/phasing/PhaserSV.cpp
new file mode 100644
index 0000000..0e61541
--- /dev/null
+++ b/src/phasing/PhaserSV.cpp
@@ -0,0 +1,13 @@
+/*
+ * PhaserSV.cpp
+ *
+ * Created on: Sep 2, 2015
+ * Author: fsedlaze
+ */
+
+#include "PhaserSV.h"
+
+void PhaserSV::phase(std::vector<Breakpoint *> &svs) {
+ //run all against all to check read names.. Is there an easier method?
+
+}
diff --git a/src/phasing/PhaserSV.h b/src/phasing/PhaserSV.h
new file mode 100644
index 0000000..8acafcc
--- /dev/null
+++ b/src/phasing/PhaserSV.h
@@ -0,0 +1,24 @@
+/*
+ * PhaserSV.h
+ *
+ * Created on: Sep 2, 2015
+ * Author: fsedlaze
+ */
+
+#include "api/BamReader.h"
+#include "../sub/Breakpoint.h"
+#include "../tree/BinTree.h"
+
+class PhaserSV{
+private:
+
+public:
+ PhaserSV(){
+
+ }
+ ~PhaserSV(){
+
+ }
+ void phase(std::vector<Breakpoint *> &svs);
+
+};
diff --git a/src/plane-sweep/IContainer.h b/src/plane-sweep/IContainer.h
new file mode 100644
index 0000000..da90405
--- /dev/null
+++ b/src/plane-sweep/IContainer.h
@@ -0,0 +1,24 @@
+/*
+ * IContainer.h
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef ICONTAINER_H_
+#define ICONTAINER_H_
+#include "Node.h"
+class IContainer{
+protected:
+
+public:
+ IContainer(){}
+ virtual ~IContainer(){};
+ virtual void push(Alignment * read)=0;
+ virtual int size()=0;
+ virtual Node * get_entries(int new_start)=0;
+};
+
+
+
+#endif /* ICONTAINER_H_ */
diff --git a/src/plane-sweep/Main.cpp b/src/plane-sweep/Main.cpp
new file mode 100644
index 0000000..0f7731e
--- /dev/null
+++ b/src/plane-sweep/Main.cpp
@@ -0,0 +1,81 @@
+//============================================================================
+// Name : plane_sweep.cpp
+// Author : Fritz Sedlazeck
+// Version :
+// Copyright : Your copyright notice
+// Description : Hello World in C++, Ansi-style
+//============================================================================
+/*
+#include <iostream>
+#include <algorithm>
+#include <chrono>
+#include "SamParser.h"
+#include "Plane-sweep.h"
+
+using namespace std::chrono;
+int main(int argc, char *argv[]) {
+ if (argc > 1) {
+ IParser *parser;
+ std::string filename = std::string(argv[2]);
+ std::cout << filename << std::endl;
+ std::string samfile = "sam";
+ std::size_t found = filename.find(samfile);
+ if (found != std::string::npos) {
+ parser = new SamParser();
+ } else {
+ parser = new TableParser();
+ }
+ std::vector<read_str> reads = parser->parse_reads(filename);
+ std::cout << "Read in reads: " << reads.size() << std::endl;
+
+ switch (atoi(argv[1])) {
+ case 1:
+ if (argc == 3) {
+ //exchangeable by implementing IParser:
+ PlaneSweep * sweep = new PlaneSweep();
+ std::cout << "start storing:" << std::endl;
+ high_resolution_clock::time_point t1 = high_resolution_clock::now();
+ for (size_t i = 0; i < reads.size(); i++) {
+ sweep->add_read(reads[i]);
+ }
+ sweep->finalyze();
+ high_resolution_clock::time_point t2 = high_resolution_clock::now();
+ auto duration = std::chrono::duration_cast<std::chrono::milliseconds>( t2 - t1 ).count();
+ std::cerr<<"time elapsed: "<<duration<<" milliseconds"<<std::endl;
+ }
+ break;
+ case 2:
+ high_resolution_clock::time_point t1 = high_resolution_clock::now();
+ int maxim = 0;
+ for (size_t i = 0; i < reads.size(); i++) {
+ maxim = std::max(maxim, (int) reads[i].stop);
+ }
+ std::cerr << maxim << std::endl;
+ std::vector<int> genome;
+ genome.resize(maxim + 1, 0);
+ for (size_t i = 0; i < reads.size(); i++) {
+ size_t j = reads[i].start;
+ while (j < reads[i].stop) {
+ genome[j]++;
+ j++;
+ }
+ }
+ for (size_t i = 0; i < genome.size(); i++) {
+ if (i > 0 && genome[i - 1] != genome[i]) {
+ std::cout << i << '\t' << genome[i] << std::endl;
+ }
+ }
+ genome.clear();
+ high_resolution_clock::time_point t2 = high_resolution_clock::now();
+ auto duration = std::chrono::duration_cast<std::chrono::milliseconds>( t2 - t1 ).count();
+ std::cerr<<"time elapsed: "<<duration<<" milliseconds"<<std::endl;
+ break;
+ }
+ }else{
+ std::cout<<"Choose:"<<std::endl;
+ std::cout<<"1: Plain sweep"<<std::endl;
+ std::cout<<"2: Brute force"<<std::endl;
+ }
+ return 0;
+}
+*/
diff --git a/src/plane-sweep/MyHeap.cpp b/src/plane-sweep/MyHeap.cpp
new file mode 100644
index 0000000..c64892a
--- /dev/null
+++ b/src/plane-sweep/MyHeap.cpp
@@ -0,0 +1,93 @@
+/*
+ * MyHeap.cpp
+ *
+ * Created on: Apr 15, 2015
+ * Author: fsedlaze
+ */
+#include "MyHeap.h"
+
+
+void MyHeap::push(Alignment * read){
+ //std::cout<<"push"<<std::endl;
+ Node * tmp=new Node(read);
+ heap.push_back(tmp);
+ heapifyup(heap.size() - 1);
+ //std::cout<<"push end"<<std::endl;
+}
+
+Node * MyHeap::get_entries(int new_start)
+{
+
+ if(!heap.empty() && heap[0]->get_stop() <= new_start){
+ Node * min = heap[0];
+ heap[0] = heap.at(heap.size() - 1);
+ heap.pop_back();
+ heapifydown(0);
+ std::cout<<"get entries end "<<this->size()<<std::endl;
+ print();
+ return min;
+ }
+ return NULL;
+}
+
+void MyHeap::print()
+{
+ std::vector<Node*>::iterator pos = heap.begin();
+ std::cout << "Heap = ";
+ while ( pos != heap.end() ) {
+ std::cout << (*pos)->get_stop() << " ";
+ ++pos;
+ }
+ std::cout << std::endl;
+}
+
+void MyHeap::heapifyup(int index)
+{
+ while ( ( index > 0 ) && ( parent(index) >= 0 ) &&
+ ( heap[parent(index)]->get_stop() > heap[index]->get_stop() ) )
+ {
+ Node * tmp = heap[parent(index)];
+ heap[parent(index)] = heap[index];
+ heap[index] = tmp;
+ index = parent(index);
+ }
+}
+
+void MyHeap::heapifydown(int index)
+{
+ int child = left(index);
+ if ( ( child > 0 ) && ( right(index) > 0 ) &&
+ ( heap[child]->get_stop() > heap[right(index)]->get_stop() ) )
+ {
+ child = right(index);
+ }
+ if ( child > 0 )
+ {
+ Node* tmp = heap[index];
+ heap[index] = heap[child];
+ heap[child] = tmp;
+ heapifydown(child);
+ }
+}
+
+int MyHeap::left(int parent)
+{
+ int i = ( parent << 1 ) + 1; // 2 * parent + 1
+ return ( i < heap.size() ) ? i : -1;
+}
+
+int MyHeap::right(int parent)
+{
+ int i = ( parent << 1 ) + 2; // 2 * parent + 2
+ return ( i < heap.size() ) ? i : -1;
+}
+
+int MyHeap::parent(int child)
+{
+ if (child != 0)
+ {
+ int i = (child - 1) >> 1;
+ return i;
+ }
+ return -1;
+}
diff --git a/src/plane-sweep/MyHeap.h b/src/plane-sweep/MyHeap.h
new file mode 100644
index 0000000..f69133e
--- /dev/null
+++ b/src/plane-sweep/MyHeap.h
@@ -0,0 +1,31 @@
+/*
+ * MyHeap.h
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef MYHEAP_H_
+#define MYHEAP_H_
+#include <vector>
+#include <iostream>
+#include "IContainer.h"
+#include "Node.h"
+class MyHeap: public IContainer{
+ private:
+ std::vector<Node*> heap;
+ int left(int parent);
+ int right(int parent);
+ int parent(int child);
+ void heapifyup(int index);
+ void heapifydown(int index);
+public:
+ MyHeap(){};
+ ~MyHeap(){};
+ void push(Alignment * read);
+ Node * get_entries(int new_start);
+ int size() { return this->heap.size(); }
+ void print();
+};
+
+#endif /* MYHEAP_H_ */
diff --git a/src/plane-sweep/MyList.cpp b/src/plane-sweep/MyList.cpp
new file mode 100644
index 0000000..45a98c5
--- /dev/null
+++ b/src/plane-sweep/MyList.cpp
@@ -0,0 +1,68 @@
+/*
+ * MyList.cpp
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#include "MyList.h"
+
+void MyList::push(Alignment * read) {
+ Node * new_node = new Node(read);
+ //std::cout << new_node->get_read()->getPosition() << " "
+ // << new_node->get_read()->getPosition() + new_node->get_stop()
+ // << std::endl;
+
+ this->num_nodes++;
+ if (this->start == NULL) { //set first node:
+ start = new_node;
+ last = new_node;
+ } else { //we have already something:
+ if (this->last != NULL
+ && new_node->get_stop() > this->last->get_stop()) { //should be inserted after the last node
+ last->set_next(new_node);
+ last = new_node;
+ } else { //insert within the list:
+ if (start->get_stop() > new_node->get_stop()) {
+ //insert before start;
+ new_node->set_next(start);
+ start = new_node;
+ } else {
+ Node * prev = start;
+ Node * curr = start->get_next();
+
+ while (curr != NULL && curr->get_stop() < new_node->get_stop()) {
+ prev = curr;
+ curr = curr->get_next();
+ }
+ new_node->set_next(curr);
+ prev->set_next(new_node);
+ }
+ }
+ }
+}
+
+Node* MyList::pop() {
+ Node * tmp = start;
+ if (tmp != NULL) {
+ start = start->get_next();
+ if (tmp == this->last) {
+ last = NULL;
+ }
+ this->num_nodes--;
+ tmp->set_next(NULL);
+ }
+ return tmp; //TODO: care about deleting the object!
+}
+
+int MyList::size() {
+ return num_nodes;// + num_nodes_lowMQ;
+}
+
+Node * MyList::get_entries(int new_start) {
+ if (new_start == -1
+ || (this->start != NULL && this->start->get_stop() <= new_start)) {
+ return pop();
+ }
+ return NULL;
+}
diff --git a/src/plane-sweep/MyList.h b/src/plane-sweep/MyList.h
new file mode 100644
index 0000000..3f842c4
--- /dev/null
+++ b/src/plane-sweep/MyList.h
@@ -0,0 +1,51 @@
+/*
+ * MyList.h
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef MYLIST_H_
+#define MYLIST_H_
+#include <iostream>
+#include "Node.h"
+#include "IContainer.h"
+class MyList: public IContainer{
+private:
+ Node * start;
+ Node * last;
+ int num_nodes;
+ int split_read;
+ int num_nodes_lowMQ;
+public:
+ MyList(){
+ split_read=0;
+ num_nodes=0;
+ num_nodes_lowMQ=0;
+ start=NULL;
+ last=NULL;
+ }
+ ~MyList(){
+ delete [] start;
+ num_nodes=0;
+ split_read=0;
+ }
+ Node * pop();
+ void push(Alignment * read);
+ int size();
+ Node * get_entries(int new_start);
+ Node * get_start(){
+ return start;
+ }
+ int get_normal_reads(){
+ return this->num_nodes; //strange value!
+ }
+ int get_split_reads(){
+ return this->split_read;
+ }
+ int get_lowMqcov(){
+ return this->num_nodes_lowMQ;
+ }
+};
+
+#endif /* MYLIST_H_ */
diff --git a/src/plane-sweep/Node.h b/src/plane-sweep/Node.h
new file mode 100644
index 0000000..469032d
--- /dev/null
+++ b/src/plane-sweep/Node.h
@@ -0,0 +1,79 @@
+/*
+ * Node.h
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef NODE_H_
+#define NODE_H_
+#include "../Alignment.h"
+
+
+class Node{
+private:
+ //Alignment * read;
+ int start;
+ int stop;
+ short mq;
+ Node * next;
+ short times; //indicate if two reads stop at the same location
+ bool processed;
+ bool support; //to del
+public:
+ Node(){
+ processed=false;
+ next=NULL;
+ times=0;
+ start=0;
+ stop=0;
+ // read=NULL;
+ mq=0;
+ }
+ Node(Alignment * read){
+ processed=false;
+ next=NULL;
+ times=1;
+ mq=read->getMappingQual();
+ //this->read=read;
+ this->start=read->getPosition();
+ this->stop=this->start+read->getRefLength();
+ support=read->supports_SV();
+ }
+ ~Node(){
+
+ }
+ bool supports_SV(){
+ return support;
+ }
+ Node * get_next(){
+ return this->next;
+ }
+ int get_stop(){
+ return this->stop;
+ }
+ void set_next(Node *next){
+ this->next=next;
+ }
+ /*Alignment * get_read(){
+ return this->read;
+ }
+ void set_read(Alignment * read){
+ this->start=read->getPosition();
+ this->stop=read->getRefLength();
+ this->read=read;
+ }*/
+ void set_processed (bool value){
+ this->processed=value;
+ }
+
+ bool was_processed(){
+ return this->processed;
+ }
+ uint16_t getmq(){
+ return mq;
+ }
+};
+
+
+#endif /* NODE_H_ */
diff --git a/src/plane-sweep/Plane-sweep.h b/src/plane-sweep/Plane-sweep.h
new file mode 100644
index 0000000..5b00b18
--- /dev/null
+++ b/src/plane-sweep/Plane-sweep.h
@@ -0,0 +1,58 @@
+/*
+ * Plane-sweep.h
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PLANESWEEP_H_
+#define PLANESWEEP_H_
+
+#include <iostream>
+#include <string.h>
+#include "IContainer.h"
+#include "../Alignment.h"
+
+#include "MyList.h"
+#include "MyHeap.h"
+//#include "MyHeap.h"
+
+class PlaneSweep {
+private:
+ MyList* current_reads; // will be used as heap.
+// void release_pos(int new_start);
+ int RefID; //the current chr/contig;
+ int curr_pos;
+public:
+ PlaneSweep(){
+ curr_pos=0;
+ RefID=-1;
+ current_reads=new MyList();
+ //current_reads=new MyHeap();
+ }
+ ~PlaneSweep(){
+ delete current_reads;
+ }
+ void release_pos(int new_start);
+ void add_read(Alignment* read);
+ void finalyze();
+ Node * get_current_reads(){
+ return current_reads->get_start();
+ }
+ int get_num_reads(){
+ return current_reads->get_normal_reads();
+ }
+ int get_num_SVreads(){
+ return current_reads->get_split_reads();
+ }
+ int get_num_lowMQ_reads(){
+ return current_reads->get_lowMqcov();
+ }
+ int get_RefID(){
+ return RefID;
+ }
+};
+
+
+
+#endif /* PLANESWEEP_H_ */
diff --git a/src/plane-sweep/PlaneSweep.cpp b/src/plane-sweep/PlaneSweep.cpp
new file mode 100644
index 0000000..dcb7c12
--- /dev/null
+++ b/src/plane-sweep/PlaneSweep.cpp
@@ -0,0 +1,43 @@
+/*
+ * PlaneSweep.cpp
+ *
+ * Created on: Apr 14, 2015
+ * Author: fsedlaze
+ */
+
+#include "Plane-sweep.h"
+
+void PlaneSweep::release_pos(int new_start) {
+
+ Node * curr = this->current_reads->get_entries(new_start);
+ while (curr != NULL) {
+ delete curr;
+ curr = NULL;
+ curr = this->current_reads->get_entries(new_start);
+ }
+
+}
+
+void PlaneSweep::add_read(Alignment* read) {
+ //std::cout<<"\tnew: "<<read->getPosition()<<" "<<read->getPosition()+read->getRefLength()<<std::endl;
+ if (this->current_reads->size() > 0 && read->getRefID() != this->RefID) {
+ finalyze();
+ this->RefID = read->getRefID();
+ } else if (this->current_reads->size() == 0) {
+ this->RefID = read->getRefID();
+ }
+ //first check if we can release already some positions:
+ release_pos(read->getPosition());
+
+ //insert read to our list if it does not support an SV:
+ if(!read->supports_SV() && read->getMappingQual()>20){
+ current_reads->push(read);
+ }
+}
+
+void PlaneSweep::finalyze() {
+ //report the remaining positions/reads;
+ std::cout << "Finalize" << std::endl;
+ release_pos(-1); //flag for releasing all;
+}
+
diff --git a/src/print/BedePrinter.cpp b/src/print/BedePrinter.cpp
new file mode 100644
index 0000000..04c6e67
--- /dev/null
+++ b/src/print/BedePrinter.cpp
@@ -0,0 +1,59 @@
+/*
+ * BedePrinter.cpp
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#include "BedePrinter.h"
+
+void BedePrinter::print_header() {
+//nothing to be done!
+}
+void BedePrinter::print_body(Breakpoint * &SV, RefVector ref) {
+ //chr1 934247 934273 chr1 934690 934692 1 1.36749e-107 + - TYPE:DELETION IDS:2,36 STRANDS:+-,36 MAX:chr1:934248;chr1:934692 95:chr1:934248-934252;chr1:934692-934692
+
+ if (!this->bed_tree.is_in(SV->get_coordinates().start.most_support, this->root) && !this->bed_tree.is_in(SV->get_coordinates().stop.most_support, this->root)) {
+ std::string chr;
+ int pos = IPrinter::calc_pos(SV->get_coordinates().start.max_pos, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos + 1);
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.max_pos, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos + 1);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", id);
+ id++;
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%c", '0'); //TODO: think about eval!
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", SV->get_strand(1).c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", IPrinter::get_type(SV->get_SVtype()).c_str()); //TODO
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", "IDS:??");
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", "STRANDS:??,666");
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", "MAX:");
+ pos = IPrinter::calc_pos(SV->get_coordinates().start.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", ':');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", ',');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", ':');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i",SV->get_length());
+ fprintf(file, "%c", '\n');
+ }
+}
diff --git a/src/print/BedePrinter.h b/src/print/BedePrinter.h
new file mode 100644
index 0000000..5b5b47a
--- /dev/null
+++ b/src/print/BedePrinter.h
@@ -0,0 +1,26 @@
+/*
+ * BedePrinter.h
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PRINT_BEDEPRINTER_H_
+#define PRINT_BEDEPRINTER_H_
+#include "IPrinter.h"
+class BedePrinter:public IPrinter{
+private:
+ void print_header();
+ void print_body(Breakpoint *& SV, RefVector ref);
+public:
+ BedePrinter(){
+
+ }
+ ~BedePrinter(){
+
+ }
+};
+
+
+
+#endif /* PRINT_BEDEPRINTER_H_ */
diff --git a/src/print/IPrinter.cpp b/src/print/IPrinter.cpp
new file mode 100644
index 0000000..1bb9f4a
--- /dev/null
+++ b/src/print/IPrinter.cpp
@@ -0,0 +1,51 @@
+/*
+ * IPrinter.cpp
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#include "IPrinter.h"
+
+std::string IPrinter::get_chr(long pos, RefVector ref) {
+ // std::cout << "pos: " << pos << std::endl;
+ size_t id = 0;
+ while (id < ref.size() && pos >= 0) {
+ pos -= ((long) ref[id].RefLength + (long) Parameter::Instance()->max_dist);
+ // std::cout << id << std::endl;
+ id++;
+ }
+
+ return ref[id - 1].RefName;
+}
+
+long IPrinter::calc_pos(long pos, RefVector ref, std::string &chr) {
+ size_t i = 0;
+ pos -= (ref[i].RefLength + Parameter::Instance()->max_dist);
+
+ while (i < ref.size() && pos >= 0) {
+ i++;
+ pos -= ((long) ref[i].RefLength + (long) Parameter::Instance()->max_dist);
+ }
+ chr = ref[i].RefName;
+ return pos + ref[i].RefLength + (long) Parameter::Instance()->max_dist;
+}
+
+std::string IPrinter::get_type(char type) {
+ if (type & DEL) {
+ return "DEL";
+ }
+ if (type & INV) {
+ return "INV";
+ }
+ if (type & DUP) {
+ return "DUP";
+ }
+ if (type & INS) {
+ return "INS";
+ }
+ if (type & TRA) {
+ return "TRA";
+ }
+ return "WTF"; // should not occur!
+}
diff --git a/src/print/IPrinter.h b/src/print/IPrinter.h
new file mode 100644
index 0000000..be10a18
--- /dev/null
+++ b/src/print/IPrinter.h
@@ -0,0 +1,65 @@
+/*
+ * IPrinter.h
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PRINT_IPRINTER_H_
+#define PRINT_IPRINTER_H_
+#include <vector>
+#include <iostream>
+
+#include "../tree/Intervall_bed.h"
+#include "api/BamReader.h"
+#include "../Ignore_Regions.h"
+#include "../sub/Breakpoint.h"
+
+class IPrinter {
+protected:
+ FILE *file;
+ uint id;
+ RefVector ref;
+ BamParser *mapped_file;
+ IntervallTree_bed bed_tree;
+ Leaf *root;
+
+ virtual void print_header()=0;
+ virtual void print_body(Breakpoint * &SV, RefVector ref)=0;
+ long calc_pos(long pos, RefVector ref, std::string &chr);
+ std::string get_chr(long pos, RefVector ref);
+ std::string get_type(char type);
+
+public:
+ IPrinter() {
+ id = 0;
+ root = NULL;
+ //we just need the ref information:
+ }
+ virtual ~IPrinter() {
+ delete mapped_file;
+ fclose(this->file);
+ }
+ ;
+ void printSV(Breakpoint * SV) {
+ print_body(SV, ref);
+ }
+ void init() {
+ if(!Parameter::Instance()->output_vcf.empty()){
+ file = fopen(Parameter::Instance()->output_vcf.c_str(), "w");
+ }else if(!Parameter::Instance()->output_maria.empty()){
+ file = fopen(Parameter::Instance()->output_maria.c_str(), "w");
+ }
+ print_header();
+ BamParser *mapped_file = new BamParser(Parameter::Instance()->bam_files[0]);
+ this->ref = mapped_file->get_refInfo();
+
+ if (!Parameter::Instance()->ignore_regions_bed.empty()) {
+ std::cout << "Cross checking..." << std::endl;
+ initialize_bed(bed_tree, root, ref);
+ }
+
+ }
+};
+
+#endif /* PRINT_IPRINTER_H_ */
diff --git a/src/print/MariaPrinter.cpp b/src/print/MariaPrinter.cpp
new file mode 100644
index 0000000..b5cc161
--- /dev/null
+++ b/src/print/MariaPrinter.cpp
@@ -0,0 +1,58 @@
+/*
+ * MariaPrinter.cpp
+ *
+ * Created on: Sep 4, 2015
+ * Author: fsedlaze
+ */
+
+#include "MariaPrinter.h"
+
+void MariaPrinter::print_header() {
+ fprintf(file, "%s", "Chrom\tstart\tstop\tchrom2\tstart2\tstop2\tvariant_name/ID\tscore (smaller is better)\tstrand1\tstrand2\ttype\tnumber_of_split_reads\tbest_chrom\tbest_pos\tbest_chrom2\tbest_pos2\tpredicted_length\n");
+}
+void MariaPrinter::print_body(Breakpoint *& SV, RefVector ref) {
+ //"Chrom\tstart\tstop\tchrom2\tstart2\tstop2\tvariant_name/ID\tscore (smaller is better)\tstrand1\tstrand2\ttype\tnumber_of_split_reads\n"
+ if (!this->bed_tree.is_in(SV->get_coordinates().start.most_support, this->root) && !this->bed_tree.is_in(SV->get_coordinates().stop.most_support, this->root)) {
+ std::string chr;
+ std::string strands = SV->get_strand(2);
+ int pos = IPrinter::calc_pos(SV->get_coordinates().start.min_pos, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", IPrinter::calc_pos(SV->get_coordinates().start.max_pos, ref, chr));
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.min_pos, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", IPrinter::calc_pos(SV->get_coordinates().stop.max_pos, ref, chr));
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", id);
+ id++;
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", -1); //TODO: score
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%c", strands[0]);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%c", strands[1]);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", IPrinter::get_type(SV->get_SVtype()).c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", SV->get_support());
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().start.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", SV->get_length());
+ fprintf(file, "%c", '\n');
+ }
+}
diff --git a/src/print/MariaPrinter.h b/src/print/MariaPrinter.h
new file mode 100644
index 0000000..e73c4ae
--- /dev/null
+++ b/src/print/MariaPrinter.h
@@ -0,0 +1,24 @@
+/*
+ * MariaPrinter.h
+ *
+ * Created on: Sep 4, 2015
+ * Author: fsedlaze
+ */
+
+#include "IPrinter.h"
+
+
+
+class MariaPrinter: public IPrinter{
+private:
+ void print_header();
+ void print_body(Breakpoint * &SV, RefVector ref);
+
+public:
+ MariaPrinter(){
+
+ }
+ ~MariaPrinter(){
+
+ }
+};
diff --git a/src/print/NGMPrinter.cpp b/src/print/NGMPrinter.cpp
new file mode 100644
index 0000000..9a5b735
--- /dev/null
+++ b/src/print/NGMPrinter.cpp
@@ -0,0 +1,70 @@
+/*
+ * NGMPrinter.cpp
+ *
+ * Created on: Sep 23, 2015
+ * Author: fsedlaze
+ */
+
+/*
+ * MariaPrinter.cpp
+ *
+ * Created on: Sep 4, 2015
+ * Author: fsedlaze
+ */
+
+#include "NGMPrinter.h"
+
+void NGMPrinter::print_header() {
+
+}
+void NGMPrinter::print_body(Breakpoint *& SV, RefVector ref) {
+ //"Chrom\tstart\tstop\tchrom2\tstart2\tstop2\tvariant_name/ID\tscore (smaller is better)\tstrand1\tstrand2\ttype\tnumber_of_split_reads\n"
+
+ if (!this->bed_tree.is_in(SV->get_coordinates().start.most_support, this->root) && !this->bed_tree.is_in(SV->get_coordinates().stop.most_support, this->root)) {
+ std::string chr;
+ std::string strands = SV->get_strand(2);
+ if ((SV->get_SVtype() & TRA) || SV->get_length() > 1000000) { //1MB??
+ int pos = IPrinter::calc_pos(SV->get_coordinates().start.min_pos, ref, chr) - Parameter::Instance()->max_dist;
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ if (pos > 0) {
+ fprintf(file, "%i", pos);
+ } else {
+ fprintf(file, "%i", 0);
+ }
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().start.max_pos, ref, chr) + Parameter::Instance()->max_dist;
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\n');
+
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.min_pos, ref, chr) - Parameter::Instance()->max_dist;
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ if (pos > 0) {
+ fprintf(file, "%i", pos);
+ } else {
+ fprintf(file, "%i", 0);
+ }
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.max_pos, ref, chr) - Parameter::Instance()->max_dist;
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\n');
+ } else { //smaller SV:
+
+ int pos = IPrinter::calc_pos(SV->get_coordinates().start.min_pos, ref, chr) - Parameter::Instance()->max_dist;
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ if (pos > 0) {
+ fprintf(file, "%i", pos);
+ } else {
+ fprintf(file, "%i", 0);
+ }
+ fprintf(file, "%c", '\t');
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.max_pos, ref, chr) + Parameter::Instance()->max_dist;
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\n');
+
+ }
+
+ }
+}
diff --git a/src/print/NGMPrinter.h b/src/print/NGMPrinter.h
new file mode 100644
index 0000000..4a56ef9
--- /dev/null
+++ b/src/print/NGMPrinter.h
@@ -0,0 +1,30 @@
+/*
+ * NGMPrinter.h
+ *
+ * Created on: Sep 23, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PRINT_NGMPRINTER_H_
+#define PRINT_NGMPRINTER_H_
+
+#include "IPrinter.h"
+class NGMPrinter:public IPrinter{
+
+private:
+ void print_header();
+ void print_body(Breakpoint * &SV, RefVector ref);
+public:
+ NGMPrinter(){
+
+ }
+
+ ~NGMPrinter(){
+ }
+
+};
+
+
+
+
+#endif /* PRINT_NGMPRINTER_H_ */
diff --git a/src/print/VCFPrinter.cpp b/src/print/VCFPrinter.cpp
new file mode 100644
index 0000000..6b2d98e
--- /dev/null
+++ b/src/print/VCFPrinter.cpp
@@ -0,0 +1,84 @@
+/*
+ * VCFPrinter.cpp
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#include "VCFPrinter.h"
+
+void VCFPrinter::print_header() {
+ fprintf(file, "%s", "##fileformat=VCFv4.1\n");
+ fprintf(file, "%s", "##fileDate=20150221\n"); //TODO change!
+ fprintf(file, "%s", "##ALT=<ID=DEL,Description=\"Deletion\">\n");
+ fprintf(file, "%s", "##ALT=<ID=DUP,Description=\"Duplication\">\n");
+ fprintf(file, "%s", "##ALT=<ID=INV,Description=\"Inversion\">\n");
+ fprintf(file, "%s", "##ALT=<ID=TRA,Description=\"Translocation\">\n");
+ fprintf(file, "%s", "##ALT=<ID=INS,Description=\"Insertion\">\n");
+ //rintf(file, "%s", "##FILTER=<ID=LowQual,Description=\"PE support below 3 or mapping quality below 20.\">\n");
+ //printf(file, "%s", "##INFO=<ID=CIEND,Number=2,Type=Integer,Description=\"PE confidence interval around END\">\n");
+ //fprintf(file, "%s", "##INFO=<ID=CIPOS,Number=2,Type=Integer,Description=\"PE confidence interval around POS\">\n");
+ fprintf(file, "%s", "##INFO=<ID=CHR2,Number=1,Type=String,Description=\"Chromosome for END coordinate in case of a translocation\">\n");
+ fprintf(file, "%s", "##INFO=<ID=END,Number=1,Type=Integer,Description=\"End position of the structural variant\">\n");
+// fprintf(file, "%s", "##INFO=<ID=PE,Number=1,Type=Integer,Description=\"Paired-end support of the structural variant\">\n");
+ fprintf(file, "%s", "##INFO=<ID=MAPQ,Number=1,Type=Integer,Description=\"Median mapping quality of paired-ends\">\n");
+ fprintf(file, "%s", "##INFO=<ID=RE,Number=1,Type=Integer,Description=\"read support\">\n");
+// fprintf(file, "%s", "##INFO=<ID=SRQ,Number=1,Type=Float,Description=\"Split-read consensus alignment quality\">\n");
+// fprintf(file, "%s", "##INFO=<ID=CONSENSUS,Number=1,Type=String,Description=\"Split-read consensus sequence\">\n");
+// fprintf(file, "%s", "##INFO=<ID=CT,Number=1,Type=String,Description=\"Paired-end signature induced connection type\">\n");
+ fprintf(file, "%s", "##INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description=\"Imprecise structural variation\">\n");
+ fprintf(file, "%s", "##INFO=<ID=PRECISE,Number=0,Type=Flag,Description=\"Precise structural variation\">\n");
+ fprintf(file, "%s", "##INFO=<ID=SVLEN,Number=1,Type=Integer,Description=\"Length of the SV\">\n");
+ fprintf(file, "%s", "##INFO=<ID=SUPTYPE,Number=1,Type=String,Description=\"What supports the SV.\">\n");
+ fprintf(file, "%s", "##INFO=<SUPTYPE=SP,Description=\"SV supported by split reads\">\n");
+ fprintf(file, "%s", "##INFO=<SUPTYPE=MD,Description=\"SV supported by MD string\">\n");
+ fprintf(file, "%s", "##INFO=<SUPTYPE=CI,Description=\"SV supported by Cigar string\">\n");
+ //fprintf(file, "%s", "##INFO=<ID=SVTYPE,Number=1,Type=String,Description=\"Type of structural variant\">\n");
+ fprintf(file, "%s", "##INFO=<ID=SVMETHOD,Number=1,Type=String,Description=\"Type of approach used to detect SV\">\n");
+ fprintf(file, "%s", "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=GL,Number=G,Type=Float,Description=\"Log10-scaled genotype likelihoods for RR,RA,AA genotypes\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=FT,Number=1,Type=String,Description=\"Per-sample genotype filter\">\n");
+// fprintf(file, "%s", "##FORMAT=<ID=RC,Number=1,Type=Integer,Description=\"Normalized high-quality read count for the SV\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=DR,Number=1,Type=Integer,Description=\"# high-quality reference pairs\">\n");
+ fprintf(file, "%s", "##FORMAT=<ID=DV,Number=1,Type=Integer,Description=\"# high-quality variant pairs\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=RR,Number=1,Type=Integer,Description=\"# high-quality reference junction reads\">\n");
+ //fprintf(file, "%s", "##FORMAT=<ID=RV,Number=1,Type=Integer,Description=\"# high-quality variant junction reads\">\n");
+ fprintf(file, "%s", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT");
+ for (size_t i = 0; i < Parameter::Instance()->bam_files.size(); i++) {
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%s", Parameter::Instance()->bam_files[i].c_str());
+ }
+ fprintf(file, "%c", '\n');
+}
+void VCFPrinter::print_body(Breakpoint * &SV, RefVector ref) {
+ if (!this->bed_tree.is_in(SV->get_coordinates().start.most_support, this->root) && !this->bed_tree.is_in(SV->get_coordinates().stop.most_support, this->root)) {
+ std::string chr;
+ int pos = IPrinter::calc_pos(SV->get_coordinates().start.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", pos);
+ fprintf(file, "%c", '\t');
+ fprintf(file, "%i", id);
+ id++;
+ fprintf(file, "%s", "\tN\t<");
+ fprintf(file, "%s", IPrinter::get_type(SV->get_SVtype()).c_str());
+ fprintf(file, "%s", ">\t.\tPASS\tIMPRECISE;SVMETHOD=Snifflesv0.0.1;CHR2=");
+ pos = IPrinter::calc_pos(SV->get_coordinates().stop.most_support, ref, chr);
+ fprintf(file, "%s", chr.c_str());
+ fprintf(file, "%s", ";END=");
+ fprintf(file, "%i", pos);
+ fprintf(file, "%s", ";SUPTYPE=");
+ fprintf(file, "%s", SV->get_supporting_types().c_str());
+ fprintf(file, "%s", ";SVLEN=");
+ fprintf(file, "%i", SV->get_length());
+ fprintf(file, "%s", ";STRANDS=");
+ fprintf(file, "%s", SV->get_strand(1).c_str());
+ fprintf(file, "%s", ";RE=");
+ fprintf(file, "%i", SV->get_support());
+ fprintf(file, "%s", "\tGT:DV\t./.:");
+ fprintf(file, "%i", SV->get_support());
+ fprintf(file, "%c", '\n');
+ }
+}
+
diff --git a/src/print/VCFPrinter.h b/src/print/VCFPrinter.h
new file mode 100644
index 0000000..9bd4c4a
--- /dev/null
+++ b/src/print/VCFPrinter.h
@@ -0,0 +1,29 @@
+/*
+ * VCFPrinter.h
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef PRINT_VCFPRINTER_H_
+#define PRINT_VCFPRINTER_H_
+
+#include "IPrinter.h"
+
+class VCFPrinter:public IPrinter{
+private:
+ void print_header();
+ void print_body(Breakpoint * &SV, RefVector ref);
+public:
+ VCFPrinter(){
+
+ }
+
+ ~VCFPrinter(){
+ }
+
+};
+
+
+
+#endif /* PRINT_VCFPRINTER_H_ */
diff --git a/src/realign/IAlignment.h b/src/realign/IAlignment.h
new file mode 100644
index 0000000..2b45bf4
--- /dev/null
+++ b/src/realign/IAlignment.h
@@ -0,0 +1,63 @@
+#ifndef __IALIGNMENT_H__
+#define __IALIGNMENT_H__
+
+struct Align {
+ Align() :
+ pBuffer1(0), pBuffer2(0), ExtendedData(0), PositionOffset(0), QStart(
+ 0), QEnd(0), Score(0.0f), Identity(0.0f), NM(0) {
+ }
+ char * pBuffer1; // = pCigar = pRef
+ char * pBuffer2; // = pMD = pQry
+ void * ExtendedData;
+ int PositionOffset; // Position in Ref, an der das Alignment beginnt
+ int QStart; // Anzahl Basen, die beim Qry am Anfang abgeschnitten wurden
+ int QEnd; // Anzahl Basen, die beim Qry am Ende abgeschnitten wurden
+ float Score;
+ float Identity;
+ int NM;
+};
+
+static int const cCookie = 0x10201130;
+
+/*
+ Anmerkung zum Parameter mode:
+
+ int AlignmentType = mode & 0xFF; // 0..Smith-Waterman, 1..Needleman-Wunsch
+ int ReportType = (mode >> 8) & 0xFF; // 0..Plain alignment (Ref+Qry), 1..SAM (Cigar+MD)
+ bool BSMappingActive = mode & 0x10000;
+
+ Anmerkung BS-Mapping:
+
+ extData zeigt bei BSMappingActive == true auf ein Flag-Array (char*) der L�nge batchSize,
+ wobei bei 0 die TC-Match-Funktion, bei 1 die AG-Match-Funktion verwendet werden soll:
+
+ if (extData[i] == 0) -> TC-Matching f�r ref/qry-Paar i
+ if (extData[i] == 1) -> AG-Matching - "" -
+
+ */
+class IAlignment {
+public:
+ virtual int GetScoreBatchSize() const = 0;
+ virtual int GetAlignBatchSize() const = 0;
+
+ virtual int BatchScore(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, float * const results,
+ void * extData) = 0;
+
+ virtual int SingleAlign(int const mode, int const corridor, char const * const refSeq, char const * const qrySeq, Align & result, void * extData) { return 0; }
+
+ virtual int SingleScore(int const mode, int const corridor, char const * const refSeq, char const * const qrySeq, float & result, void * extData) { return 0; }
+
+ virtual int BatchAlign(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, Align * const results,
+ void * extData) = 0;
+
+ virtual ~IAlignment() {}
+};
+
+typedef IAlignment * (*pfCreateAlignment)(int const gpu_id);
+typedef void (*pfDeleteAlignment)(IAlignment*);
+
+#endif
diff --git a/src/realign/Realign.cpp b/src/realign/Realign.cpp
new file mode 100644
index 0000000..edf0e89
--- /dev/null
+++ b/src/realign/Realign.cpp
@@ -0,0 +1,172 @@
+/*
+ * Realign.cpp
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#include "Realign.h"
+
+void Realigner::init() {
+ //run through ref sequence and store a file * at the begining of each chr;
+ myfile.open(Parameter::Instance()->ref_seq.c_str(), ifstream::in);
+ if (!myfile.good()) {
+ cout << "Fastq Parser: could not open file: " << Parameter::Instance()->ref_seq.c_str() << endl;
+ exit(0);
+ }
+
+ buffer_size = 20000;
+ buffer = new char[buffer_size];
+
+ myfile.getline(buffer, buffer_size);
+
+ long len = 0;
+ while (!myfile.eof()) {
+ if (buffer[0] == '>') {
+ ref_str tmp;
+ tmp.length = len;
+ tmp.file_pos = myfile.tellg();
+ meta_info.push_back(tmp);
+ } else {
+ for (size_t i = 0; i < buffer_size && buffer[i] != '\0'; i++) {
+ len++;
+ }
+ }
+ myfile.getline(buffer, buffer_size);
+ }
+ myfile.close();
+}
+std::string Realigner::read_new_part(long start, long stop) {
+ long pos = start;
+ int i = 0;
+ for (; i < (int) meta_info.size() && start - meta_info[i].length > 0; i++) {
+ }
+ i--; //one step back
+ start -= meta_info[i].length;
+ stop -= meta_info[i].length;
+ myfile.open(Parameter::Instance()->ref_seq.c_str(), ifstream::in);
+ myfile.seekg(meta_info[i].file_pos);
+ myfile.getline(buffer, buffer_size);
+ string seq;
+ pos = 0;
+ while (!myfile.eof() && buffer[0] != '>') {
+ for (size_t i = 0; i < buffer_size && buffer[i] != '\0' && buffer[i] != '\n' && buffer[0] != '>'; i++) {
+ if (pos >= start && pos <= stop) {
+ seq += toupper(buffer[i]);
+ }
+ pos++;
+ }
+ myfile.getline(buffer, buffer_size);
+ }
+ if (buffer[0] != '>') {
+ for (size_t i = 0; i < buffer_size && buffer[i] != '\0' && buffer[i] != '\n' && buffer[0] != '>'; i++) {
+ if (pos >= start && pos <= stop) {
+ seq += toupper(buffer[i]);
+ }
+ pos++;
+ }
+ }
+
+ myfile.close();
+ return seq;
+}
+
+std::string Realigner::read_chr(short id) {
+
+ myfile.open(Parameter::Instance()->ref_seq.c_str(), ifstream::in);
+ myfile.seekg(meta_info[id].file_pos);
+ myfile.getline(buffer, buffer_size);
+ string seq;
+ while (!myfile.eof() && buffer[0] != '>') {
+ for (size_t i = 0; i < buffer_size && buffer[i] != '\0' && buffer[i] != '\n'; i++) {
+ seq += toupper(buffer[i]);
+ }
+ myfile.getline(buffer, buffer_size);
+ }
+ myfile.close();
+ return seq;
+}
+
+void get_coords_DEL(Breakpoint *& sv) {
+ //not much todo unless size is large -> tra???
+ region_ref_str tmp;
+
+ if(sv->get_length()<5000){
+ //take normal ref region;
+ //tmp.start=sv->get_coordinates().start;
+ }else{
+ //chop:
+ }
+}
+void get_coords_DUP(Breakpoint *& sv) {
+ //duplicate the dup region next to each other?? -> define 2 overlapping regions
+}
+void get_coords_TRA(Breakpoint *& sv) {
+ //define 2 regions on the chr
+}
+void get_coords_INS(Breakpoint *& sv) {
+ //nothing todo
+}
+void get_coords_INV(Breakpoint *& sv) {
+ // Estimate 2 breakpoints and set directions!
+}
+void Realigner::align(std::vector<Breakpoint *> sv) {
+ long len = 0;
+ std::string seq;
+ //1 collect regions from ref
+ for (size_t i = 0; i < sv.size(); i++) {//parallel
+
+ //check if split read-> concatenate sequences;
+ //else just a standard alignment.
+
+ //split to ease job:
+ if (sv[i]->get_SVtype() & DEL) {
+ get_coords_DEL(sv[i]);
+ } else if (sv[i]->get_SVtype() & DUP) {
+ get_coords_DUP(sv[i]);
+ } else if (sv[i]->get_SVtype() & TRA) {
+ get_coords_TRA(sv[i]);
+ } else if (sv[i]->get_SVtype() & INV) {
+ get_coords_INV(sv[i]);
+ } else if (sv[i]->get_SVtype() & INS) {
+ get_coords_INS(sv[i]);
+ }
+ }
+ //2: collect regions from ref:
+ //for each chr run through all SV??
+ for (size_t i = 0; i < this->meta_info.size(); i++) {
+ long curr_length=meta_info[i].length;
+
+ std::string ref = this->read_chr(i);
+ long next_length=curr_length+ref.size();
+
+ for(size_t j = 0; j < sv.size(); j++) { //parallel
+ /*for(size_t k=0;k<sv[j]->get_ref_coord().size();k++){
+ if((sv[j]->get_ref_coord()[k].start-curr_length) >0 && (sv[j]->get_ref_coord()[k].stop-next_length)<0){
+ //extract region:
+ int start=sv[j]->get_ref_coord()[k].start-curr_length;
+ int stop=sv[j]->get_ref_coord()[k].stop-curr_length;
+ sv[j]->set_ref_seq(k,ref.substr(start,stop-start));
+ }
+ }*/
+ }
+ }
+ for (size_t i = 0; i < sv.size(); i++) {//parallel
+ //2 send SV+ regions to alignment using OpenMP
+ IAlignment * aligner = new SWCPUCor(0);
+ Align align;
+ align.pBuffer1 = new char[400];
+ align.pBuffer2 = new char[400];
+ char * refSeq = new char[400];
+ char * readSeq = new char[400];
+ int mode = 0;
+ int cigarLength = aligner->SingleAlign(mode, Parameter::Instance()->corridor, refSeq, readSeq, align, 0);
+
+ //3 Backtrack information
+
+ //4 filter out SV if necessary
+
+ }
+
+
+}
diff --git a/src/realign/Realign.h b/src/realign/Realign.h
new file mode 100644
index 0000000..e482317
--- /dev/null
+++ b/src/realign/Realign.h
@@ -0,0 +1,40 @@
+/*
+ * Realign.h
+ *
+ * Created on: Aug 24, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef REALIGN_REALIGN_H_
+#define REALIGN_REALIGN_H_
+#include "../sub/Breakpoint.h"
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include "IAlignment.h"
+#include "SWCPU.h"
+struct ref_str {
+ long length;
+ streampos file_pos;
+};
+
+class Realigner {
+private:
+ std::string read_new_part(long start, long stop);
+ vector<ref_str> meta_info;
+ size_t buffer_size;
+ char*buffer;
+ ifstream myfile;
+ void init();
+ std::string read_chr(short id);
+public:
+ Realigner() {
+ init();
+ }
+ ~Realigner() {
+ delete[] buffer;
+ }
+ void align(std::vector<Breakpoint *> SV);
+};
+
+#endif /* REALIGN_REALIGN_H_ */
diff --git a/src/realign/SWCPU.cpp b/src/realign/SWCPU.cpp
new file mode 100644
index 0000000..35a6ab4
--- /dev/null
+++ b/src/realign/SWCPU.cpp
@@ -0,0 +1,563 @@
+/*
+ * SWCPU.cpp
+ *
+ * Created on: Jun 15, 2011
+ * Author: fritz
+ */
+
+#include "SWCPU.h"
+
+//TODO: hack to pass data for debug output
+Align cur_align;
+
+SWCPUCor::SWCPUCor(int gpu_id) {
+
+// cigar = bool(((gpu_id >> 8) & 0xFF) == 1);
+
+ batch_size = 1;
+
+ mat = 2.0f;
+ mis = -5.0f;
+ gap_open_read = -5.0f;
+ gap_open_ref = -5.0f;
+ gap_ext = -5.0f;
+ gap_decay = 0.05f;
+ gap_ext_min = -1.0f;
+
+ long maxLen = (long) 100000 * (long) 20000;
+ alignMatrix = new MatrixElement[maxLen];
+
+ fprintf(stderr, "Allocationg: %llu\n", maxLen * sizeof(MatrixElement));
+
+ binaryCigar = new int[200000];
+
+// short temp[6][6] = { mat, mis, mis, mis, 0, mis, mis, mat, mis, mis, 0, mis,
+// mis, mis, mat, mis, 0, mis, mis, mis, mis, mat, 0, mis, 0, 0, 0, 0,
+// 0, 0, 0, 0, 0, 0, 0, mat };
+// memcpy(scores, temp, 6 * 6 * sizeof(short));
+
+ fprintf(stderr, "SWCPU initialized\n");
+}
+
+SWCPUCor::~SWCPUCor() {
+ delete[] alignMatrix;
+ alignMatrix = 0;
+ delete[] binaryCigar;
+ binaryCigar = 0;
+}
+
+Score SWCPUCor::SW_Score(char const * const refSeqList,
+ char const * const qrySeqList, int * fwResults, int corr_length,
+ MatrixElement * mat_pointer) {
+
+// memset(local_mat_line, 0, corr_length * sizeof(short));
+ char const * scaff = refSeqList;
+ char const * read = qrySeqList;
+
+ //Init matrix lines
+ MatrixElement * matrix = mat_pointer;
+
+ //Init matrix lines
+ for (int i = 0; i < corr_length; ++i) {
+ //local_mat_line[i] = 0;
+ matrix[i].direction = CIGAR_STOP;
+ matrix[i].indelRun = 0;
+ matrix[i].score = 0;
+ }
+ matrix[corr_length].direction = CIGAR_STOP;
+ matrix[corr_length].indelRun = 0;
+ matrix[corr_length].score = 0;
+
+ Score curr_max = -1.0f;
+ int read_index = 0;
+
+ int x = 0;
+ for (; *read != line_end; ++read) {
+ char read_char_cache = *read;
+ matrix += (corr_length + 1);
+// short left_cell = 0;
+ matrix[0].direction = CIGAR_STOP;
+ matrix[0].indelRun = 0;
+ matrix[0].score = 0;
+
+ for (int ref_index = 0; ref_index < corr_length - 1; ++ref_index) {
+
+ MatrixElement & diag = matrix[-(corr_length + 1) + ref_index + 1];
+ MatrixElement & up = matrix[-(corr_length + 1) + ref_index + 2];
+ MatrixElement & left = matrix[ref_index];
+
+ bool eq = read_char_cache == scaff[ref_index];
+ Score diag_cell = diag.score + ((eq) ? mat : mis);
+
+ Score up_cell = 0;
+ Score left_cell = 0;
+
+ int ins_run = 0;
+ int del_run = 0;
+
+ if (up.direction == CIGAR_I) {
+ ins_run = up.indelRun;
+ if (up.score == 0) {
+ up_cell = 0;
+ } else {
+ up_cell = up.score
+ + std::min(gap_ext_min,
+ gap_ext + ins_run * gap_decay);
+ }
+ } else {
+ up_cell = up.score + gap_open_read;
+ }
+
+ if (left.direction == CIGAR_D) {
+ del_run = left.indelRun;
+ if (left.score == 0) {
+ left_cell = 0;
+ } else {
+ left_cell = left.score
+ + std::min(gap_ext_min,
+ gap_ext + del_run * gap_decay);
+ }
+ } else {
+ left_cell = left.score + gap_open_ref;
+ }
+
+ //find max
+ Score max_cell = 0;
+ max_cell = max(left_cell, max_cell);
+ max_cell = max(diag_cell, max_cell);
+ max_cell = max(up_cell, max_cell);
+
+ MatrixElement & current = matrix[(ref_index + 1)];
+ if (del_run > 0 && max_cell == left_cell) {
+ current.score = max_cell;
+ current.direction = CIGAR_D;
+ current.indelRun = del_run + 1;
+ } else if (ins_run > 0 && max_cell == up_cell) {
+ current.score = max_cell;
+ current.direction = CIGAR_I;
+ current.indelRun = ins_run + 1;
+ } else if (max_cell == diag_cell) {
+ current.score = max_cell;
+ if (eq) {
+ current.direction = CIGAR_EQ;
+ } else {
+ current.direction = CIGAR_X;
+ }
+ current.indelRun = 0;
+ } else if (max_cell == left_cell) {
+ current.score = max_cell;
+ current.direction = CIGAR_D;
+ current.indelRun = 1;
+ } else if (max_cell == up_cell) {
+ current.score = max_cell;
+ current.direction = CIGAR_I;
+ current.indelRun = 1;
+ } else {
+ current.score = 0;
+ current.direction = CIGAR_STOP;
+ current.indelRun = 0;
+ }
+
+ if (max_cell > curr_max) {
+ curr_max = max_cell;
+ fwResults[param_best_ref_index] = ref_index;
+ fwResults[param_best_read_index] = read_index;
+ fwResults[3] = curr_max;
+
+ }
+ }
+ matrix[corr_length].direction = CIGAR_STOP;
+ matrix[corr_length].score = 0;
+ matrix[corr_length].indelRun = 0;
+
+ scaff++;
+ read_index += 1;
+ }
+ fwResults[2] = (read_index - fwResults[0]) - 1;
+ if (read_index == 0) {
+ fwResults[0] = fwResults[1] = 2;
+ }
+ return curr_max;
+}
+
+int SWCPUCor::printCigarElement(char const op, int const length, char * cigar) {
+ int offset = 0;
+ offset = sprintf(cigar, "%d%c", length, op);
+ return offset;
+}
+
+int SWCPUCor::computeCigarMD(Align & result, int const gpuCigarOffset,
+ int const * const gpuCigar, char const * const refSeq, int corr_length,
+ int read_length, int const QStart, int const QEnd) {
+ int alignment_length = corr_length + read_length + 1;
+
+ int finalCigarLength = 0;
+
+ int cigar_offset = 0;
+// int md_offset = 0;
+
+ if (((gpuCigar[gpuCigarOffset] >> 4) + QStart) > 0) {
+ fprintf(stderr, "Adding %d to QSTart\n", QStart);
+ result.QStart = (gpuCigar[gpuCigarOffset] >> 4) + QStart;
+ cigar_offset += printCigarElement('S', result.QStart,
+ result.pRef + cigar_offset);
+ finalCigarLength += result.QStart;
+ }
+
+ int cigar_m_length = 0;
+// int md_eq_length = 0;
+ int ref_index = 0;
+ for (int j = gpuCigarOffset + 1; j < (alignment_length - 1); ++j) {
+ int op = gpuCigar[j] & 15;
+ int length = gpuCigar[j] >> 4;
+
+ //debugCigar(op, length);
+
+ switch (op) {
+ case CIGAR_X:
+ cigar_m_length += length;
+
+ //Produces: [0-9]+(([A-Z]+|\^[A-Z]+)[0-9]+)*
+ //instead of: [0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)*
+// md_offset += sprintf(result.pQry + md_offset, "%d", md_eq_length);
+// for (int k = 0; k < length; ++k) {
+// md_offset += sprintf(result.pQry + md_offset, "%c",
+// refSeq[ref_index++]);
+// }
+// md_eq_length = 0;
+
+ break;
+ case CIGAR_EQ:
+ cigar_m_length += length;
+// md_eq_length += length;
+ ref_index += length;
+ break;
+ case CIGAR_D:
+ if (cigar_m_length > 0) {
+ cigar_offset += printCigarElement('M', cigar_m_length,
+ result.pRef + cigar_offset);
+ finalCigarLength += cigar_m_length;
+ cigar_m_length = 0;
+ }
+ cigar_offset += printCigarElement('D', length,
+ result.pRef + cigar_offset);
+
+// md_offset += sprintf(result.pQry + md_offset, "%d", md_eq_length);
+// md_eq_length = 0;
+// result.pQry[md_offset++] = '^';
+// for (int k = 0; k < length; ++k) {
+// result.pQry[md_offset++] = refSeq[ref_index++];
+// }
+
+ break;
+ case CIGAR_I:
+ if (cigar_m_length > 0) {
+ cigar_offset += printCigarElement('M', cigar_m_length,
+ result.pRef + cigar_offset);
+ finalCigarLength += cigar_m_length;
+ cigar_m_length = 0;
+ }
+ cigar_offset += printCigarElement('I', length,
+ result.pRef + cigar_offset);
+ finalCigarLength += length;
+
+ break;
+ default:
+ fprintf(stderr, "Invalid cigar string: %d\n", op);
+ std::cout << "Offset: " << gpuCigarOffset << std::endl;
+ for (int x = 0; x < alignment_length * 2; ++x) {
+ std::cout << gpuCigar[x] << " ";
+ }
+ std::cout << std::endl;
+ exit(1);
+ }
+ }
+// md_offset += sprintf(result.pQry + md_offset, "%d", md_eq_length);
+ if (cigar_m_length > 0) {
+ cigar_offset += printCigarElement('M', cigar_m_length,
+ result.pRef + cigar_offset);
+ finalCigarLength += cigar_m_length;
+ cigar_m_length = 0;
+ }
+
+ if (((gpuCigar[alignment_length - 1] >> 4) + QEnd) > 0) {
+ fprintf(stderr, "Adding %d to QEnd\n", QEnd);
+ result.QEnd = (gpuCigar[alignment_length - 1] >> 4) + QEnd;
+ cigar_offset += printCigarElement('S', result.QEnd,
+ result.pRef + cigar_offset);
+ finalCigarLength += result.QEnd;
+
+ }
+ //TODO: fix
+ result.Identity = 1.0f;
+ result.pRef[cigar_offset] = '\0';
+// result.pQry[md_offset] = '\0';
+
+ return finalCigarLength;
+}
+
+bool SWCPUCor::Backtracking_CIGAR(char const * const scaff,
+ char const * const read, int *& fwdResults, int *& alignments,
+ int corr_length, int read_length, int alignment_length,
+ MatrixElement * mat_pointer) {
+
+ bool valid = true;
+
+ MatrixElement * matrix = mat_pointer;
+
+ int best_read_index = fwdResults[param_best_read_index];
+ int best_ref_index = fwdResults[param_best_ref_index];
+
+ int cigarLenth = 0;
+ int cigarLengthCheck = 0;
+
+ int totalDelLength = 0;
+ int totalINsLength = 0;
+
+ int minCorridor = corr_length * 0.01f;
+ int maxCorridor = corr_length - minCorridor;
+
+ if (best_read_index > 0) {
+ matrix += (((corr_length + 1) * (best_read_index + 1)));
+
+ int abs_ref_index = best_ref_index + best_read_index;
+ int alignment_index = alignment_length - 1;
+
+ int pointer = CIGAR_STOP;
+ int cigar_element = CIGAR_S;
+ int cigar_length = fwdResults[qend];
+ cigarLenth += fwdResults[qend];
+ while ((pointer = matrix[(best_ref_index + 1)].direction) != CIGAR_STOP) {
+// Log.Message("Best ref index: %d (%d)", best_ref_index + 1, (corr_length + 1));
+// printf("%s\t%d\t%d\t%d\t%d\n", (char *) cur_align.ExtendedData,
+// cur_align.NM, best_read_index, best_ref_index + 1,
+// corr_length + 1);
+ if (best_ref_index <= minCorridor
+ || best_ref_index >= maxCorridor) {
+ fprintf(stderr, "Corridor probably too small\n");
+ valid = false;
+// getchar();
+ }
+
+ if (pointer == CIGAR_X || pointer == CIGAR_EQ) {
+ matrix -= ((corr_length + 1));
+ best_read_index -= 1;
+ abs_ref_index -= 1;
+
+ cigarLenth += 1;
+ } else if (pointer == CIGAR_I) {
+ matrix -= ((corr_length + 1));
+ best_read_index -= 1;
+ best_ref_index += 1;
+
+ cigarLenth += 1;
+ } else if (pointer == CIGAR_D) {
+ best_ref_index -= 1;
+ abs_ref_index -= 1;
+
+ } else {
+ fprintf(stderr, "Error in backtracking. Invalid CIGAR operation found\n");
+ exit(1);
+
+ }
+
+ if (pointer == cigar_element) {
+ cigar_length += 1;
+ } else {
+ alignments[alignment_index--] = (cigar_length << 4
+ | cigar_element);
+ if (cigar_element != CIGAR_D) {
+ cigarLengthCheck += cigar_length;
+ }
+
+ cigar_element = pointer;
+ cigar_length = 1;
+ }
+ }
+ alignments[alignment_index--] = (cigar_length << 4 | cigar_element);
+ if (cigar_element != CIGAR_D) {
+ cigarLengthCheck += cigar_length;
+ }
+
+ alignments[alignment_index] = ((best_read_index + 1) << 4 | CIGAR_S);
+ cigarLengthCheck += (best_read_index + 1);
+ cigarLenth += (best_read_index + 1);
+ fwdResults[ref_position] = abs_ref_index + 1;
+ fwdResults[qstart] = best_read_index + 1;
+ //qend was set by "forward" kernel
+ fwdResults[alignment_offset] = alignment_index;
+
+ if (cigarLenth != cigarLengthCheck) {
+ fprintf(stderr, "Error in CIGAR length: %d vs %d\n", cigarLenth, cigarLengthCheck);
+ } else {
+ if (read_length != cigarLenth) {
+ fprintf(stderr, "Error read length != cigar length: %d vs %d\n", read_length, cigarLenth);
+ exit(1);
+ }
+ }
+ fprintf(stderr, "Read length: %d, CIGAR length: %d\n", read_length, cigarLenth);
+ }
+ return valid;
+}
+
+int SWCPUCor::GetScoreBatchSize() const {
+ return 0;
+}
+int SWCPUCor::GetAlignBatchSize() const {
+ return 0;
+}
+
+int SWCPUCor::BatchAlign(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, Align * const results,
+ void * extData) {
+
+ throw "Not implemented";
+
+ fprintf(stderr, "Unsupported alignment mode %i\n", mode);
+ return 0;
+}
+
+void SWCPUCor::print_matrix(int alignment_length, const char* const refSeq,
+ int read_length, const char* const qrySeq, int corr_length,
+ MatrixElement* mat_pointer) {
+
+ printf(" - ");
+ for (int x = 0; x < alignment_length - 1; ++x) {
+ printf(" %c ", refSeq[x]);
+ }
+ printf("\n");
+ for (size_t row = 0; row < read_length + 1; ++row) {
+ if (row == 0) {
+ printf("-: ");
+ } else {
+ printf("%c: ", qrySeq[row - 1]);
+ }
+ for (int x = 0; x < row; ++x) {
+ printf(" ");
+ }
+ for (size_t col = 0; col < corr_length + 1; ++col) {
+ MatrixElement* cell = mat_pointer + (row * (corr_length + 1) + col);
+ printf("%*d ", 2, cell->indelRun);
+ }
+ printf("\n");
+ }
+
+ printf(" - ");
+ for (int x = 0; x < alignment_length - 1; ++x) {
+ printf(" %c ", refSeq[x]);
+ }
+ printf("\n");
+ for (size_t row = 0; row < read_length + 1; ++row) {
+ if (row == 0) {
+ printf("-: ");
+ } else {
+ printf("%c: ", qrySeq[row - 1]);
+ }
+ for (int x = 0; x < row; ++x) {
+ printf(" ");
+ }
+ for (size_t col = 0; col < corr_length + 1; ++col) {
+ MatrixElement* cell = mat_pointer + (row * (corr_length + 1) + col);
+ printf("%*d ", 2, cell->direction);
+ }
+ printf("\n");
+ }
+
+ printf(" - ");
+ for (int x = 0; x < alignment_length - 1; ++x) {
+ printf(" %c ", refSeq[x]);
+ }
+ printf("\n");
+ for (size_t row = 0; row < read_length + 1; ++row) {
+ if (row == 0) {
+ printf("-: ");
+ } else {
+ printf("%c: ", qrySeq[row - 1]);
+ }
+ for (int x = 0; x < row; ++x) {
+ printf(" ");
+ }
+ for (size_t col = 0; col < corr_length + 1; ++col) {
+ MatrixElement* cell = mat_pointer + (row * (corr_length + 1) + col);
+ printf("%*.*f ", 2, 0, cell->score);
+ }
+ printf("\n");
+ }
+}
+
+int SWCPUCor::SingleAlign(int const mode, int const corridor,
+ char const * const refSeq, char const * const qrySeq, Align & align,
+ void * extData) {
+
+// Log.Message("Aligning: ");
+// Log.Message("%s", refSeq);
+// Log.Message("%s", qrySeq);
+
+ cur_align = align;
+
+ int * clipping = 0;
+ if (extData == 0) {
+ clipping = new int[2];
+ clipping[0] = 0;
+ clipping[1] = 0;
+ } else {
+ clipping = (int *) extData;
+ }
+
+ int read_length = strlen(qrySeq);
+ fprintf(stderr, "Read length (single align) is %d\n", read_length);
+ align.pBuffer1 = new char[read_length * 4];
+ // align.pBuffer2 = new char[read_length * 4];
+ align.pBuffer2 = new char[1];
+ align.pBuffer2[0] = '\0';
+
+ int finalCigarLength = 0;
+
+ int corr_length = corridor;
+ int alignment_length = (corr_length + read_length + 1);
+
+ int * fwdResults = new int[result_number];
+
+ Score score = SW_Score(refSeq, qrySeq, fwdResults, corr_length,
+ alignMatrix);
+
+// print_matrix(alignment_length, refSeq, read_length, qrySeq, corr_length,
+// alignMatrix);
+// Log.Message("%d, %d, %d, %d", fwdResults[0], fwdResults[1], fwdResults[2], fwdResults[3]);
+
+ bool valid = Backtracking_CIGAR(refSeq, qrySeq, fwdResults, binaryCigar,
+ corr_length, read_length, alignment_length, alignMatrix);
+
+ if (valid) {
+ finalCigarLength = computeCigarMD(align, fwdResults[3], binaryCigar,
+ refSeq + fwdResults[0], corr_length, read_length, clipping[0],
+ clipping[1]);
+ align.PositionOffset = fwdResults[0];
+ align.Score = score;
+ }
+
+ delete[] fwdResults;
+
+ if (extData == 0) {
+ delete[] clipping;
+ clipping = 0;
+ }
+
+ if (!valid) {
+ finalCigarLength = -1;
+ }
+
+ return finalCigarLength;
+
+}
+
+int SWCPUCor::BatchScore(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, float * const results,
+ void * extData) {
+ throw "Not implemented";
+ return 0;
+
+}
+
diff --git a/src/realign/SWCPU.h b/src/realign/SWCPU.h
new file mode 100644
index 0000000..20f0050
--- /dev/null
+++ b/src/realign/SWCPU.h
@@ -0,0 +1,119 @@
+/*
+ * SWCPU.h
+ *
+ * Created on: Jun 15, 2011
+ * Author: fritz
+ */
+
+#ifndef SWCPU_H_
+#define SWCPU_H_
+
+#define pRef pBuffer1
+#define pQry pBuffer2
+
+#include "IAlignment.h"
+
+#include <iostream>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+using std::endl;
+using std::cout;
+using std::max;
+
+#define CIGAR_STOP 10
+#define short_min -16000
+#define result_number 4
+#define line_end '\0'
+#define ref_position 0
+#define qstart 1
+#define qend 2
+#define alignment_offset 3
+#define param_best_read_index 0
+#define param_best_ref_index 1
+
+#define CIGAR_M 0
+#define CIGAR_I 1
+#define CIGAR_D 2
+#define CIGAR_N 3
+#define CIGAR_S 4
+#define CIGAR_H 5
+#define CIGAR_P 6
+#define CIGAR_EQ 7
+#define CIGAR_X 8
+
+
+typedef float Score;
+
+struct MatrixElement {
+ Score score;
+ int indelRun;
+ char direction;
+};
+
+
+const char trans[256] = { 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0,
+ 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 5, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 5, 4, 4, 4,
+ 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 };
+
+class SWCPUCor: public IAlignment {
+public:
+ SWCPUCor(int gpu_id);
+ virtual ~SWCPUCor();
+ virtual int GetScoreBatchSize() const;
+ virtual int GetAlignBatchSize() const;
+ virtual int BatchScore(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, float * const results,
+ void * extData);
+ virtual int BatchAlign(int const mode, int const batchSize,
+ char const * const * const refSeqList,
+ char const * const * const qrySeqList, Align * const results,
+ void * extData);
+
+ virtual int SingleAlign(int const mode, int const corridor,
+ char const * const refSeq, char const * const qrySeq,
+ Align & result, void * extData);
+private:
+
+ //bool cigar;
+ //short scores[6][6];
+ Score mat;
+ Score mis;
+ Score gap_open_read;
+ Score gap_open_ref;
+ Score gap_ext;
+ Score gap_ext_min;
+ Score gap_decay;
+
+ MatrixElement * alignMatrix;
+ int * binaryCigar;
+
+ //meta info
+ unsigned int batch_size; //effictive thread number that is started per call
+
+ int printCigarElement(char const op, int const length, char * cigar);
+
+ int computeCigarMD(Align & result, int const gpuCigarOffset,
+ int const * const gpuCigar, char const * const refSeq, int corr_length, int read_length, int const QStart, int const QEnd);
+
+ Score SW_Score(char const * const scaff, char const * const read, int * result, int corr_length, MatrixElement * mat_pointer);
+
+ bool Backtracking_CIGAR(char const * const scaff, char const * const read,
+ int *& result, int *& alignments, int corr_length, int read_length, int alignment_length, MatrixElement * mat_pointer);
+
+ void print_matrix(int alignment_length, const char* const refSeq,
+ int read_length, const char* const qrySeq, int corr_length,
+ MatrixElement* mat_pointer);
+};
+
+#endif /* SWCPU_H_ */
diff --git a/src/sub/Breakpoint.cpp b/src/sub/Breakpoint.cpp
new file mode 100644
index 0000000..3c232b1
--- /dev/null
+++ b/src/sub/Breakpoint.cpp
@@ -0,0 +1,412 @@
+/*
+ * Breakpoint.cpp
+ *
+ * Created on: Sep 1, 2015
+ * Author: fsedlaze
+ */
+
+/*
+ * Breakpoint.h
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+#include "Breakpoint.h"
+
+//TODO define region object and inherit from that. Plus define avoid region objects for mappability problems.
+
+std::string Breakpoint::translate_strand(pair<bool, bool> strand_pair) {
+ if (strand_pair.first && strand_pair.second) {
+ return "++";
+ } else if (strand_pair.first && !strand_pair.second) {
+ return "+-";
+ } else if (!strand_pair.first && strand_pair.second) {
+ return "-+";
+ } else if (!strand_pair.first && !strand_pair.second) {
+ return "--";
+ }
+ return " ";
+}
+
+void Breakpoint::summarize_type(char SV, std::vector<short>& array) {
+ //std::string ss;
+ if (SV & DEL) {
+ // ss += "DEL; ";
+ array[0]++;
+ }
+
+ if (SV & DUP) {
+ // ss += "DUP; ";
+ array[1]++;
+ }
+
+ if (SV & INS) {
+ // ss += "INS; ";
+ array[2]++;
+ }
+
+ if (SV & INV) {
+ // ss += "INV; ";
+ array[3]++;
+ }
+
+ if (SV & TRA) {
+ // ss += "TRA; ";
+ array[4]++;
+ }
+ //return ss;
+}
+
+char Breakpoint::get_SVtype() {
+ if (sv_type == ' ') {
+ std::cout << "was not set" << std::endl;
+ calc_support();
+ predict_SV();
+ }
+ return this->sv_type;
+}
+
+void Breakpoint::calc_support() {
+ std::vector<short> SV;
+ for (size_t i = 0; i < 5; i++) {
+ SV.push_back(0);
+ }
+ //run over all supports and check the majority type:
+ for (std::map<std::string, read_str>::iterator i = positions.support.begin(); i != positions.support.end(); i++) {
+ summarize_type((*i).second.SV, SV);
+ }
+ //given the majority type get the stats:
+ this->sv_type = eval_type(SV);
+}
+char Breakpoint::eval_type(std::vector<short> SV) {
+ std::stringstream ss;
+
+ if (SV[0] != 0) {
+ ss << " DEL(";
+ ss << SV[0];
+ ss << ")";
+ }
+ if (SV[1] != 0) {
+ ss << " DUP(";
+ ss << SV[1];
+ ss << ")";
+ }
+ if (SV[2] != 0) {
+ ss << " INS(";
+ ss << SV[2];
+ ss << ")";
+ }
+ if (SV[3] != 0) {
+ ss << " INV(";
+ ss << SV[3];
+ ss << ")";
+ }
+ if (SV[4] != 0) {
+ ss << " TRA(";
+ ss << SV[4];
+ ss << ")";
+ }
+ this->sv_debug = ss.str(); //only for debug!
+ //std::cout << sv_debug << std::endl;
+
+ int maxim = 0;
+ int id = 0;
+ for (size_t i = 0; i < SV.size(); i++) {
+ if (maxim < SV[i]) {
+ maxim = SV[i];
+ id = i;
+ } else if (maxim == SV[i]) {
+ id = -1;
+ }
+ }
+ this->type_support = maxim;
+ switch (id) {
+ case 0:
+ return DEL;
+ break;
+ case 1:
+ return DUP;
+ break;
+ case 2:
+ return INS;
+ break;
+ case 3:
+ return INV;
+ break;
+ case 4:
+ return TRA;
+ break;
+ }
+ return 'n';
+}
+
+long Breakpoint::overlap(Breakpoint * tmp) {
+ /*if ((*this->positions.support.begin()).second.SV & INS) {
+ std::cout << "\tCurr: " << this->positions.start.min_pos << " " << this->positions.stop.max_pos << endl;
+ std::cout << "\tNEW: " << tmp->get_coordinates().start.min_pos << " " << tmp->get_coordinates().stop.max_pos << endl;
+ }*/
+ if (abs(tmp->get_coordinates().start.min_pos - positions.start.min_pos) < Parameter::Instance()->max_dist && abs(tmp->get_coordinates().stop.max_pos - positions.stop.max_pos) < Parameter::Instance()->max_dist) {
+ /*if ((*this->positions.support.begin()).second.SV & INS) {
+ std::cout << "\tSAME" << std::endl;
+ }*/
+ return 0;
+ }
+ /*else {
+ if ((*this->positions.support.begin()).second.SV & INS) {
+ std::cout << "\tDiff" << std::endl;
+ }
+ }*/
+//as abstraction lets try the start+stop coordinate!
+ return (tmp->get_coordinates().start.min_pos - positions.start.min_pos) + (tmp->get_coordinates().stop.max_pos - positions.stop.max_pos);
+}
+
+void Breakpoint::predict_SV() {
+ bool md = false;
+ bool cigar = false;
+ bool split = false;
+ int num = 0;
+ std::map<long, int> starts;
+ std::map<long, int> stops;
+ std::map<std::string, int> strands;
+
+ for (std::map<std::string, read_str>::iterator i = positions.support.begin(); i != positions.support.end(); i++) {
+ if ((*i).second.SV & this->sv_type) { ///check type
+ if (starts.find((*i).second.coordinates.first) == starts.end()) {
+ starts[(*i).second.coordinates.first] = 1;
+ } else {
+ starts[(*i).second.coordinates.first]++;
+ }
+
+ if (stops.find((*i).second.coordinates.second) == stops.end()) {
+ stops[(*i).second.coordinates.second] = 1;
+ } else {
+ stops[(*i).second.coordinates.second]++;
+ }
+
+ std::string tmp = translate_strand((*i).second.strand); //std::string tmp=
+ //std::cout << tmp << std::endl;
+ if (strands.find(tmp) == strands.end()) {
+ strands[tmp] = 1;
+ } else {
+ strands[tmp]++;
+ }
+
+ if ((*i).second.type == 0) {
+ cigar = true;
+ } else if ((*i).second.type == 1) {
+ md = true;
+ } else if ((*i).second.type == 2) {
+ split = true;
+ } else {
+ std::cerr << "Type " << (*i).second.type << std::endl;
+ }
+ num++;
+ }
+ }
+
+ int maxim = 0;
+ long coord = 0;
+ for (map<long, int>::iterator i = starts.begin(); i != starts.end(); i++) {
+ //cout<<"start:\t"<<(*i).first<<" "<<(*i).second<<endl;
+ if ((*i).second > maxim) {
+ coord = (*i).first;
+ maxim = (*i).second;
+ }
+ }
+ //std::cout<<coord<<"\t"<<maxim<<endl;
+ //if(this->sv_type & INV){
+ // std::cout<<"Break: INV: "<<coord;
+ //}
+ this->positions.start.most_support = coord;
+
+ maxim = 0;
+ coord = 0;
+ for (map<long, int>::iterator i = stops.begin(); i != stops.end(); i++) {
+ if ((*i).second > maxim) {
+ coord = (*i).first;
+ maxim = (*i).second;
+ }
+ }
+ this->positions.stop.most_support = coord;
+
+ if (!(this->get_SVtype() & INS)) {
+ if (this->get_SVtype() & INS) {
+ std::cout << " ERROR " << std::endl;
+ }
+ this->length = this->positions.stop.most_support - this->positions.start.most_support;
+ }
+ starts.clear();
+ stops.clear();
+
+ for (size_t i = 0; i < strands.size(); i++) {
+ maxim = 0;
+ std::string id;
+ for (std::map<std::string, int>::iterator j = strands.begin(); j != strands.end(); j++) {
+ if (maxim < (*j).second) {
+ maxim = (*j).second;
+ id = (*j).first;
+ //std::cout << '\t' << id << std::endl;
+ }
+ }
+ if (maxim > 0) {
+ this->strand.push_back(id);
+ strands[id] = 0;
+ }
+ }
+ strands.clear();
+
+ this->supporting_types = "";
+ if (md) {
+ this->supporting_types += "MD";
+ }
+ if (cigar) {
+ if (!supporting_types.empty()) {
+ this->supporting_types += ",";
+ }
+ this->supporting_types += "CI";
+ }
+ if (split) {
+ if (!supporting_types.empty()) {
+ this->supporting_types += ",";
+ }
+ this->supporting_types += "SR";
+ }
+ //this->strand = eval_strand(strand);
+}
+
+std::string Breakpoint::to_string(RefVector ref) {
+
+ std::stringstream ss;
+ ss << "(";
+ ss << get_chr(get_coordinates().start.min_pos, ref);
+ ss << ":";
+ ss << calc_pos(get_coordinates().start.min_pos, ref);
+ ss << "-";
+ ss << get_chr(get_coordinates().stop.max_pos, ref);
+ ss << ":";
+ ss << calc_pos(get_coordinates().stop.max_pos, ref);
+ ss << " ";
+ ss << positions.support.size();
+ ss << " ";
+ ss << this->sv_debug;
+ ss << " ";
+ ss << this->get_strand(2);
+ ss << "\n";
+ int num = 0;
+ for (std::map<std::string, read_str>::iterator i = positions.support.begin(); i != positions.support.end() && num < Parameter::Instance()->report_n_reads; i++) {
+ ss << "\t";
+ ss << (*i).first;
+ ss << " ";
+ ss << (*i).second.type;
+ if ((*i).second.strand.first) {
+ ss << "+";
+ } else {
+ ss << "-";
+ }
+ if ((*i).second.strand.second) {
+ ss << "+";
+ } else {
+ ss << "-";
+ }
+ num++;
+ ss << "\n";
+ }
+ ss << " ";
+ return ss.str();
+}
+
+void Breakpoint::add_read(Breakpoint * point) {
+ if (point != NULL) {
+
+ this->positions.start.min_pos = min(this->positions.start.min_pos, point->get_coordinates().start.min_pos);
+ this->positions.start.max_pos = max(this->positions.start.max_pos, point->get_coordinates().start.max_pos);
+
+ this->positions.stop.min_pos = min(this->positions.stop.min_pos, point->get_coordinates().stop.min_pos);
+ this->positions.stop.max_pos = max(this->positions.stop.max_pos, point->get_coordinates().stop.max_pos);
+ bool flag = false;
+ std::map<std::string, read_str> support = point->get_coordinates().support;
+ for (std::map<std::string, read_str>::iterator i = support.begin(); i != support.end(); i++) {
+ if (this->positions.support[(*i).first].SV & INS) {
+ flag = true;
+ }
+ this->positions.support[(*i).first] = (*i).second;
+ }
+ if (flag) {
+ this->length = (this->length + point->get_length()) / 2;
+ }
+ }
+}
+
+std::string Breakpoint::get_chr(long pos, RefVector ref) {
+// std::cout << "pos: " << pos << std::endl;
+ size_t id = 0;
+ while (id < ref.size() && pos >= 0) {
+ pos -= (long) ref[id].RefLength;
+ // std::cout << id << std::endl;
+ id++;
+ }
+
+ return ref[id - 1].RefName;
+}
+
+long Breakpoint::calc_pos(long pos, RefVector ref) {
+ size_t i = 0;
+ pos -= ref[i].RefLength;
+ while (i < ref.size() && pos >= 0) {
+ i++;
+ pos -= ref[i].RefLength;
+ }
+ return pos + ref[i].RefLength;
+}
+
+int Breakpoint::get_support() {
+ return type_support;
+}
+char complement(char nuc) {
+ switch (nuc) {
+ case 'A':
+ return 'T';
+ break;
+ case 'C':
+ return 'G';
+ break;
+ case 'G':
+ return 'C';
+ break;
+ case 'T':
+ return 'A';
+ break;
+ default:
+ return nuc;
+ break;
+ }
+}
+std::string Breakpoint::rev_complement(std::string seq) {
+ std::string tmp;
+ for (std::string::reverse_iterator i = seq.rbegin(); i != seq.rend(); i++) {
+ tmp += complement((*i));
+ }
+ return tmp;
+}
+
+std::string Breakpoint::get_strand(int num_best) {
+ //if(this->strand.empty()){
+ // predict_SV();
+ //}
+ if (sv_type == ' ') {
+ // std::cout<<"was not set"<<std::endl;
+ calc_support();
+ predict_SV();
+ }
+ std::string tmp = this->strand[0];
+ for (int i = 1; i < num_best; i++) {
+ tmp += '\t';
+ if (i < (int) this->strand.size()) {
+ tmp += this->strand[i];
+ } else {
+ tmp += ' ';
+ }
+ }
+ return tmp;
+}
+
diff --git a/src/sub/Breakpoint.h b/src/sub/Breakpoint.h
new file mode 100644
index 0000000..fc6e902
--- /dev/null
+++ b/src/sub/Breakpoint.h
@@ -0,0 +1,149 @@
+/*
+ * Breakpoint.h
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef SUB_BREAKPOINT_H_
+#define SUB_BREAKPOINT_H_
+
+#include <string>
+#include <vector>
+#include <iostream>
+#include <sstream>
+#include <math.h>
+#include "../Paramer.h"
+#include "../BamParser.h"
+#include "../tree/BinTree.h"
+
+const unsigned char DEL = 0x01; // hex for 0000 0001
+const unsigned char DUP = 0x02; // hex for 0000 0010
+const unsigned char INS = 0x04; // hex for 0000 0100
+const unsigned char INV = 0x08; // hex for 0000 1000
+const unsigned char TRA = 0x10; // hex for 0001 0000
+
+struct region_ref_str{ //not very nice!
+ std::string read_seq;
+ long read_aln_pos;
+ bool direction;
+ long start;
+ long stop;
+ std::string ref;
+};
+
+struct svs_breakpoint_str{
+ long min_pos;
+ long max_pos;
+ long most_support;
+};
+struct read_str {
+ //to identify
+ std::string name;
+ region_ref_str aln; //maybe we can use this!
+ short type; //split reads, cigar or md string
+ //for later assessment:
+ pair<bool, bool> strand;
+ pair<long,long> coordinates; // I could use the bin tree for that!
+ char SV; // bit vector
+};
+struct position_str {
+ svs_breakpoint_str start;
+ svs_breakpoint_str stop;
+ //int pos; //the chromosomes are encoded over the positions.
+ std::map<std::string,read_str> support;
+ //std::vector<read_str> support; // map?? -> no duplicated reads, easy to catch up which read is included.
+ int coverage;
+ int lowmq_cov;
+ short read_start;
+ short read_stop;
+};
+
+
+//TODO define region object and inherit from that. Plus define avoid region objects for mappability problems.
+class Breakpoint {
+private:
+ position_str positions;
+ std::vector<std::string> strand;
+ std::string supporting_types;
+ char sv_type;
+ std::string sv_debug;
+ std::string ref_seq;
+ std::vector<short> support;
+ double cov;
+ short type_support;
+ //for phasing:
+ BinTree grouped;
+ tree_node * grouped_node;
+ long length;
+
+ void summarize_support(short type);
+ //void summarize_strand(pair<bool, bool> strand, std::vector<short>& array);
+ void summarize_type(char SV, std::vector<short>& array);
+ //std::string translate_strand(short id);
+ char eval_type(std::vector<short> SV);
+ std::string rev_complement(std::string seq);
+ bool is_in(short id);
+ std::string translate_strand(pair<bool, bool> strand);
+public:
+ Breakpoint(position_str sv, int cov,long len) {
+ sv_type=' ';
+ type_support=-1;
+ this->positions = sv;
+ //std::cout << "Break1: " << positions.start << " " << positions.stop << std::endl;
+
+ //std::cout << "Break2: " << positions.start << " " << positions.stop<< std::endl;
+ this->cov = cov;
+ this->grouped_node=NULL;
+ this->set_length(len);
+ }
+ ~Breakpoint() {
+
+ }
+
+ int get_support();
+ long overlap(Breakpoint * tmp);
+ position_str get_coordinates() {
+ return this->positions;
+ }
+ void predict_SV();
+ std::string to_string(RefVector ref);
+
+ void add_read(Breakpoint * point);
+
+ double get_cov() {
+ return cov;
+ }
+ std::string get_chr(long pos, RefVector ref);
+ long calc_pos(long pos, RefVector ref);
+ char get_SVtype();
+ std::string get_strand(int num_best);
+ std::string get_ref_seq() {
+ return this->ref_seq;
+ }
+ void set_ref_seq(std::string seq) {
+ this->ref_seq = seq;
+ }
+ long get_length(){
+ return length;
+ }
+
+ void set_length(long len){
+ this->length=len;
+ }
+ std::string get_supporting_types(){
+ return this->supporting_types;
+ }
+
+ void add_grouped(int id){
+ this->grouped.insert(this->grouped_node, id);
+ }
+ vector<int> get_groupted(){
+ vector<int> tmp;
+ this->grouped.get_nodes(this->grouped_node,tmp);
+ return tmp;
+ }
+ void calc_support();
+};
+
+#endif /* SUB_BREAKPOINT_H_ */
diff --git a/src/sub/Container.h b/src/sub/Container.h
new file mode 100644
index 0000000..f04c3a0
--- /dev/null
+++ b/src/sub/Container.h
@@ -0,0 +1,15 @@
+/*
+ * Container.h
+ *
+ * Created on: Jun 30, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef SUB_CONTAINER_H_
+#define SUB_CONTAINER_H_
+
+
+
+
+
+#endif /* SUB_CONTAINER_H_ */
diff --git a/src/sub/Detect_Breakpoints.cpp b/src/sub/Detect_Breakpoints.cpp
new file mode 100644
index 0000000..87002bf
--- /dev/null
+++ b/src/sub/Detect_Breakpoints.cpp
@@ -0,0 +1,465 @@
+/*
+ * Detect_Breapoints.cpp
+ *
+ * Created on: Jun 19, 2015
+ * Author: fsedlaze
+ */
+
+#include "Detect_Breakpoints.h"
+#include "../tree/IntervallTree.h"
+#include "../tree/TNode.h"
+
+long get_ref_lengths(int id, RefVector ref) {
+ long length = 0;
+
+ for (size_t i = 0; i < (size_t) id && i < ref.size(); i++) {
+ length += ref[i].RefLength + Parameter::Instance()->max_dist;
+ }
+ return length;
+}
+
+bool should_be_stored(Breakpoint *& point) {
+ point->calc_support();
+ if (point->get_SVtype() & TRA) {
+ return (point->get_support() > 2); // this is needed as we take each chr independently and just look at the primary alignment
+ } else {
+ point->predict_SV();
+ return (point->get_support() > Parameter::Instance()->min_support && point->get_length() > Parameter::Instance()->min_length);
+ }
+}
+void flush_tree(IntervallTree & local_tree, TNode *& local_root, IntervallTree & final, TNode *& root_final, long pos) {
+ IntervallTree tmp_tree;
+ TNode *tmp_root = NULL;
+ std::vector<Breakpoint *> points;
+ local_tree.get_breakpoints(local_root, points);
+ clarify(points);
+ for (int i = 0; i < points.size(); i++) {
+ if (abs(pos - points[i]->get_coordinates().start.min_pos) < 100000 || abs(pos - points[i]->get_coordinates().stop.max_pos) < 100000) { //TODO arbitrary threshold!
+ tmp_tree.insert(points[i], tmp_root);
+ } else if (points[i]->get_support() > Parameter::Instance()->min_support && points[i]->get_length() > Parameter::Instance()->min_length) {
+ final.insert(points[i], root_final);
+ }
+ }
+ local_tree.makeempty(local_root);
+ local_tree = tmp_tree;
+ local_root = tmp_root;
+}
+
+void detect_breakpoints(std::string read_filename, IPrinter *& printer) {
+ estimate_parameters(read_filename);
+ BamParser * mapped_file = 0;
+ RefVector ref;
+ if (read_filename.find("bam") != string::npos) {
+ mapped_file = new BamParser(read_filename);
+ ref = mapped_file->get_refInfo();
+ } else {
+ cerr << "File Format not recognized. File must be a sorted .bam file!" << endl;
+ exit(0);
+ }
+//Using PlaneSweep to comp coverage and iterate through reads:
+ //PlaneSweep * sweep = new PlaneSweep();
+ std::cout << "start parsing..." << std::endl;
+//Using Interval tree to store and manage breakpoints:
+
+ IntervallTree final;
+ TNode * root_final = NULL;
+ int current_RefID = 0;
+
+ IntervallTree bst;
+ TNode *root = NULL;
+
+ Alignment * tmp_aln = mapped_file->parseRead(Parameter::Instance()->min_mq);
+ long ref_space = get_ref_lengths(tmp_aln->getRefID(), ref);
+ while (!tmp_aln->getSequence().first.empty()) {
+
+ if ((tmp_aln->getAlignment()->IsPrimaryAlignment()) && (!(tmp_aln->getAlignment()->AlignmentFlag & 0x800) && tmp_aln->get_is_save())) {
+ //flush_tree(local_tree, local_root, final, root_final);
+ if (current_RefID != tmp_aln->getRefID()) { // Regular scan through the SV and move those where the end point lies far behind the current pos or reads. Eg. 1MB?
+ current_RefID = tmp_aln->getRefID();
+ ref_space = get_ref_lengths(tmp_aln->getRefID(), ref);
+ std::cout << "Switch Chr " << ref[tmp_aln->getRefID()].RefName << " " << ref[tmp_aln->getRefID()].RefLength << std::endl;
+ std::vector<Breakpoint *> points;
+ // clarify(points);
+ bst.get_breakpoints(root, points);
+ for (int i = 0; i < points.size(); i++) {
+ if (should_be_stored(points[i])) {
+ if (points[i]->get_SVtype() & TRA) {
+ final.insert(points[i], root_final);
+ } else {
+ printer->printSV(points[i]);
+ }
+ }
+ }
+ bst.makeempty(root);
+ }
+ std::vector<str_event> cigar_event;
+ std::vector<str_event> md_event;
+ std::vector<aln_str> split_events;
+
+ if (tmp_aln->getMappingQual() > Parameter::Instance()->min_mq) {
+#pragma omp parallel // starts a new team
+ {
+#pragma omp sections
+ {
+ {
+ //if (Parameter::Instance()->useMD_CIGAR) {
+ cigar_event = tmp_aln->get_events_CIGAR();
+ //}
+ }
+#pragma omp section
+ {
+ //if (Parameter::Instance()->useMD_CIGAR) {
+ md_event = tmp_aln->get_events_MD(20);
+ //}
+ }
+#pragma omp section
+ {
+ split_events = tmp_aln->getSA(ref);
+ }
+ }
+ }
+ tmp_aln->set_supports_SV((cigar_event.empty() && md_event.empty()) && split_events.empty());
+
+ //sweep->add_read(tmp_aln);
+
+ //maybe flush the tree after each chr.....?
+
+ int cov = 0; //sweep->get_num_reads();
+
+ //maybe just store the extreme intervals for coverage -> If the cov doubled within Xbp or were the coverage is 0.
+ add_events(tmp_aln, cigar_event, 0, ref_space, bst, root, cov, tmp_aln->getQueryBases());
+ add_events(tmp_aln, md_event, 1, ref_space, bst, root, cov, tmp_aln->getQueryBases());
+ add_splits(tmp_aln, split_events, 2, ref, bst, root, cov, tmp_aln->getQueryBases());
+ }
+ }
+ mapped_file->parseReadFast(Parameter::Instance()->min_mq, tmp_aln);
+ }
+ //filter and copy results:
+ std::cout << "Finalizing .." << std::endl;
+ std::vector<Breakpoint *> points;
+ clarify(points);
+ bst.get_breakpoints(root, points);
+ for (int i = 0; i < points.size(); i++) {
+ if (should_be_stored(points[i])) {
+ if (points[i]->get_SVtype() & TRA) {
+ final.insert(points[i], root_final);
+ } else {
+ printer->printSV(points[i]);
+ }
+ }
+ }
+ bst.makeempty(root);
+
+// sweep->finalyze();
+ points.clear();
+ final.get_breakpoints(root_final, points);
+ //std::cout<<"Points: "<<points.size()<<endl;
+ //clarify(points);
+ for(size_t i =0;i<points.size();i++){
+ if (points[i]->get_SVtype() & TRA) {
+ points[i]->calc_support();
+ points[i]->predict_SV();
+ }
+ if (points[i]->get_support() > Parameter::Instance()->min_support && points[i]->get_length() > Parameter::Instance()->min_length) {
+ printer->printSV(points[i]);
+ }
+ }
+}
+
+void add_events(Alignment * tmp, std::vector<str_event> events, short type, long ref_space, IntervallTree & bst, TNode *&root, int cov, std::string read_seq) {
+ bool flag = (strcmp(tmp->getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0);
+ for (size_t i = 0; i < events.size(); i++) {
+ position_str svs;
+ //position_str stop;
+ read_str read;
+ //read.name = tmp->getName();
+ read.type = type;
+ read.SV = 0;
+ //start.support.push_back(read); //not very nice!
+ //stop.support.push_back(read);
+ if (flag) {
+ std::cout << tmp->getName() << " " << tmp->getRefID() << " " << events[i].pos << " " << abs(events[i].length) << std::endl;
+ }
+ svs.start.min_pos = (long) events[i].pos;
+ svs.start.min_pos += ref_space;
+ svs.stop.max_pos = svs.start.min_pos;
+ if (events[i].length > 0) { //length ==- length for insertions!
+ svs.stop.max_pos += events[i].length;
+ }
+ if (tmp->getStrand()) {
+ read.strand.first = (tmp->getStrand());
+ read.strand.second = !(tmp->getStrand());
+ } else {
+ read.strand.first = !(tmp->getStrand());
+ read.strand.second = (tmp->getStrand());
+ }
+ // start.support[0].read_start.min = events[i].read_pos;
+
+ if (type == 0 && events[i].length < 0) {
+ read.SV |= INS; //insertion
+ } else if (type == 0) {
+ read.SV |= DEL; //deletion
+ } else {
+ read.SV |= DEL;
+ read.SV |= INV;
+ }
+
+ if (flag) {
+ std::cout << tmp->getName() << " " << tmp->getRefID() << " " << svs.start.min_pos - ref_space << " " << svs.stop.max_pos - ref_space << std::endl;
+ }
+
+ if (svs.start.min_pos > svs.stop.max_pos) {
+ //can this actually happen?
+ read.coordinates.first = svs.stop.max_pos;
+ read.coordinates.second = svs.start.min_pos;
+ } else {
+ read.coordinates.first = svs.start.min_pos;
+ read.coordinates.second = svs.stop.max_pos;
+ }
+
+ svs.start.max_pos = svs.start.min_pos;
+ svs.stop.min_pos = svs.stop.max_pos;
+
+ if (svs.start.min_pos > svs.stop.max_pos) {
+ //maybe we have to invert the directions???
+ svs_breakpoint_str pos = svs.start;
+ svs.start = svs.stop;
+ svs.stop = pos;
+
+ pair<bool, bool> tmp = read.strand;
+
+ read.strand.first = tmp.second;
+ read.strand.second = tmp.first;
+
+ //read.strand.first = !tmp.first;
+ //read.strand.second = !tmp.second;
+ }
+
+ //TODO: we might not need this:
+ if (svs.start.min_pos > svs.stop.max_pos) {
+ read.coordinates.first = svs.stop.max_pos;
+ read.coordinates.second = svs.start.min_pos;
+ } else {
+ read.coordinates.first = svs.start.min_pos;
+ read.coordinates.second = svs.stop.max_pos;
+ }
+ svs.support[tmp->getName()] = read;
+ Breakpoint * point = new Breakpoint(svs, cov, std::abs(events[i].length));
+ bst.insert(point, root);
+ }
+}
+
+void add_splits(Alignment * tmp, std::vector<aln_str> events, short type, RefVector ref, IntervallTree & bst, TNode *&root, int cov, std::string read_seq) {
+ /* bool flag = false;
+ if (Parameter::Instance()->overlaps(ref[tmp->getRefID()].RefName,
+ tmp->getPosition(), tmp->getPosition() + tmp->getRefLength())) {
+ Parameter::Instance()->read_name = tmp->getName();
+ flag = true;
+ }
+ */
+
+ bool flag = (strcmp(tmp->getName().c_str(), Parameter::Instance()->read_name.c_str()) == 0);
+ for (size_t i = 1; i < events.size() && events.size() < Parameter::Instance()->max_splits; i++) {
+ if (flag) {
+ std::cout << "Genome pos: " << tmp->getName() << " ";
+ if (events[i - 1].strand) {
+ std::cout << "+";
+ } else {
+ std::cout << "-";
+ }
+ std::cout << events[i - 1].pos << " " << events[i - 1].pos + events[i - 1].length << " p2: ";
+
+ if (events[i - 1].strand) {
+ std::cout << "+";
+ } else {
+ std::cout << "-";
+ }
+ std::cout << events[i].pos << " " << events[i].pos + events[i].length << " Ref: " << events[i - 1].RefID << " " << events[i].RefID << std::endl;
+ }
+ position_str svs;
+ //position_str stop;
+ read_str read;
+ //read.name = tmp->getName();
+ read.type = 2;
+ read.SV = 0;
+ //stop.support.push_back(read);
+ //they mimic paired end sequencing:
+ if (events[i].RefID == events[i - 1].RefID) {
+
+ //TODO: changed because of test:
+ if (events[i - 1].strand == events[i].strand) {
+ if (events[i - 1].strand) {
+ read.strand.first = events[i - 1].strand;
+ read.strand.second = !events[i].strand;
+ } else {
+ read.strand.first = !events[i - 1].strand;
+ read.strand.second = events[i].strand;
+ }
+
+ svs.read_start = events[i - 1].read_pos_start + events[i - 1].length;
+ svs.read_stop = events[i].read_pos_start;
+ if (events[i - 1].strand) {
+ svs.start.min_pos = events[i - 1].pos + events[i - 1].length + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + get_ref_lengths(events[i].RefID, ref);
+ } else {
+ svs.start.min_pos = events[i].pos + events[i].length + get_ref_lengths(events[i].RefID, ref);
+ svs.stop.max_pos = events[i - 1].pos + get_ref_lengths(events[i - 1].RefID, ref);
+ }
+
+ if ((svs.start.min_pos - svs.stop.max_pos) > 100) {
+ read.SV |= DUP;
+ //TODO ADDED &&(svs.read_stop - svs.read_start) > (Parameter::Instance()->min_cigar_event * 2)
+ } else if (abs(svs.stop.max_pos - svs.start.min_pos) + (Parameter::Instance()->min_cigar_event * 2) < (svs.read_stop - svs.read_start) && (svs.read_stop - svs.read_start) > (Parameter::Instance()->min_cigar_event * 2)) {
+ read.SV |= INS;
+ } else if (abs(svs.stop.max_pos - svs.start.min_pos) > (svs.read_stop - svs.read_start) + (Parameter::Instance()->min_cigar_event * 2)) {
+ read.SV |= DEL;
+ } else {
+ read.SV = 'n';
+ }
+ } else { // if first part of read is in a different direction as the second part-> INV
+ read.strand.first = events[i - 1].strand;
+ read.strand.second = !events[i].strand;
+ read.SV |= INV;
+ if (events[i - 1].strand) {
+ //std::cout<<events[i].pos<<"\t"<<events[i].RefID<<"\t"<<get_ref_lengths(events[i].RefID, ref)<<endl;
+ svs.start.min_pos = events[i - 1].pos + events[i - 1].length + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + events[i].length + get_ref_lengths(events[i].RefID, ref);
+ } else {
+
+ svs.start.min_pos = events[i - 1].pos + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + get_ref_lengths(events[i].RefID, ref);
+ }
+ }
+ } else { //if not on the same chr-> TRA
+ read.strand.first = events[i - 1].strand;
+ read.strand.second = !events[i].strand;
+ if (events[i - 1].strand == events[i].strand) {
+ if (events[i - 1].strand) {
+ svs.start.min_pos = events[i - 1].pos + events[i - 1].length + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + get_ref_lengths(events[i].RefID, ref);
+ } else {
+ svs.start.min_pos = events[i - 1].pos + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + events[i].length + get_ref_lengths(events[i].RefID, ref);
+ }
+ } else {
+ if (events[i - 1].strand) {
+ svs.start.min_pos = events[i - 1].pos + events[i - 1].length + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + events[i].length + get_ref_lengths(events[i].RefID, ref);
+ } else {
+ svs.start.min_pos = events[i - 1].pos + get_ref_lengths(events[i - 1].RefID, ref);
+ svs.stop.max_pos = events[i].pos + get_ref_lengths(events[i].RefID, ref);
+ }
+ }
+ read.SV |= TRA;
+ }
+
+ if (flag) {
+ std::cout << tmp->getName() << " start: " << svs.start.min_pos << " stop: " << svs.stop.max_pos;
+ if (events[i - 1].strand) {
+ std::cout << " +";
+ } else {
+ std::cout << " -";
+ }
+ if (events[i].strand) {
+ std::cout << " +";
+ } else {
+ std::cout << " -";
+ }
+ std::cout << std::endl;
+ }
+ if (read.SV != 'n') {
+ //std::cout<<"split"<<std::endl;
+
+ svs.start.max_pos = svs.start.min_pos;
+ svs.stop.min_pos = svs.stop.max_pos;
+ if (svs.start.min_pos > svs.stop.max_pos) {
+ //maybe we have to invert the directions???
+ svs_breakpoint_str pos = svs.start;
+ svs.start = svs.stop;
+ svs.stop = pos;
+
+ pair<bool, bool> tmp = read.strand;
+
+ read.strand.first = tmp.second;
+ read.strand.second = tmp.first;
+
+ //read.strand.first = !tmp.first;
+ //read.strand.second = !tmp.second;
+ }
+
+ //TODO: we might not need this:
+ if (svs.start.min_pos > svs.stop.max_pos) {
+ read.coordinates.first = svs.stop.max_pos;
+ read.coordinates.second = svs.start.min_pos;
+ } else {
+ read.coordinates.first = svs.start.min_pos;
+ read.coordinates.second = svs.stop.max_pos;
+ }
+
+ svs.support[tmp->getName()] = read;
+ Breakpoint * point = new Breakpoint(svs, cov, events[i].length);
+ bst.insert(point, root);
+ }
+ }
+}
+
+void clarify(std::vector<Breakpoint *> & points) {
+//if WTF regions next to duplications-> delete!
+ /*for(size_t i=0;i<points.size();i++){
+
+ }*/
+}
+
+void estimate_parameters(std::string read_filename) {
+ BamParser * mapped_file = 0;
+ RefVector ref;
+ if (read_filename.find("bam") != string::npos) {
+ mapped_file = new BamParser(read_filename);
+ ref = mapped_file->get_refInfo();
+ } else {
+ cerr << "File Format not recognized. File must be a sorted .bam file!" << endl;
+ exit(0);
+ }
+
+ Alignment * tmp_aln = mapped_file->parseRead(Parameter::Instance()->min_mq);
+ double num = 0;
+ double avg_score = 0;
+ double avg_mis = 0;
+ double avg_indel = 0;
+
+ while (!tmp_aln->getSequence().first.empty() && num < 1000) {
+
+ if ((tmp_aln->getAlignment()->IsPrimaryAlignment()) && (!(tmp_aln->getAlignment()->AlignmentFlag & 0x800) && tmp_aln->get_is_save())) {
+
+ //get score ratio
+ double score = tmp_aln->get_scrore_ratio();
+ if (score != -1) {
+ avg_score += score;
+ } else {
+ avg_score += avg_score / num;
+ }
+ //cout<<"Para:\t"<<score;
+ //get avg mismatches
+ std::string md = tmp_aln->get_md();
+ if (!md.empty()) {
+ avg_mis += tmp_aln->get_num_mismatches(md);
+ //cout<<"\t"<<tmp_aln->get_num_mismatches(md);
+ }
+ //cigar threshold: (without 1!)
+ avg_indel += tmp_aln->get_avg_indel_length_Cigar();
+ //cout<<"\t"<<tmp_aln->get_avg_indel_length_Cigar()<<endl;
+ num++;
+ }
+ mapped_file->parseReadFast(Parameter::Instance()->min_mq, tmp_aln);
+ }
+ std::cout << avg_indel / num << std::endl;
+
+ Parameter::Instance()->min_num_mismatches = 0.3; //(avg_mis / num) * 0.3; //previously: 0.3
+ Parameter::Instance()->min_cigar_event = 40; //(avg_indel / num) * 20; //previously: 20
+ Parameter::Instance()->score_treshold = 2; //(avg_score / num) / 2; //previously: 2 //2
+
+ std::cout << "score: " << Parameter::Instance()->score_treshold << std::endl;
+ std::cout << "md: " << Parameter::Instance()->min_num_mismatches << std::endl;
+ std::cout << "indel: " << Parameter::Instance()->min_cigar_event << std::endl;
+
+}
diff --git a/src/sub/Detect_Breakpoints.h b/src/sub/Detect_Breakpoints.h
new file mode 100644
index 0000000..d2d50db
--- /dev/null
+++ b/src/sub/Detect_Breakpoints.h
@@ -0,0 +1,29 @@
+/*
+ * Detect_Breakpoints.h
+ *
+ * Created on: Jun 19, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef SUB_DETECT_BREAKPOINTS_H_
+#define SUB_DETECT_BREAKPOINTS_H_
+
+#include "../BamParser.h"
+#include "../Parser.h"
+#include "../Alignment.h"
+#include "../plane-sweep/Plane-sweep.h"
+#include "../tree/IntervallTree.h"
+#include "../tree/TNode.h"
+#include "../Paramer.h"
+#include "../print/IPrinter.h"
+#include <iostream>
+#include <omp.h>
+
+void clarify(std::vector<Breakpoint *> & points);
+void detect_breakpoints(std::string filename, IPrinter *& printer);
+//void screen_for_events(Node * list,IntervallTree & bst ,TNode *&root, int cov, int lowMQ_cov,RefVector ref);
+bool screen_for_events(Alignment * tmp, IntervallTree & bst, TNode *&root, RefVector ref, int cov);
+void add_events(Alignment * tmp, std::vector<str_event> events, short type, long ref_space, IntervallTree & bst, TNode *&root, int cov,std::string read_seq);
+void add_splits(Alignment * tmp, std::vector<aln_str> events, short type, RefVector ref, IntervallTree & bst, TNode *&root, int cov,std::string read_seq);
+void estimate_parameters(std::string read_filename);
+#endif /* SUB_DETECT_BREAKPOINTS_H_ */
diff --git a/src/sub/IRegion.h b/src/sub/IRegion.h
new file mode 100644
index 0000000..a668e8e
--- /dev/null
+++ b/src/sub/IRegion.h
@@ -0,0 +1,56 @@
+/*
+ * IRegion.h
+ *
+ * Created on: Aug 27, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef SUB_IREGION_H_
+#define SUB_IREGION_H_
+#include "../Paramer.h"
+#include "../BamParser.h"
+struct read_str {
+ //to identify
+ std::string name;
+ short type;
+ //for later assessment:
+ pair<bool, bool> strand;
+ char SV; // bits vector
+};
+struct position_str {
+ long start;
+ long stop;
+ //int pos; //the chromosomes are encoded over the positions.
+ std::vector<read_str> support;
+ int coverage;
+ int lowmq_cov;
+ short read_start;
+ short read_stop;
+};
+
+class IRegion {
+protected:
+ position_str start;
+
+public:
+ IRegion(position_str reg) {
+ this->start = reg;
+ //std::cout << "Break1: " << start.start << " " << start.stop << std::endl;
+ if (reg.start > reg.stop) {
+ this->start.start = reg.stop;
+ this->start.stop = reg.start;
+ }
+ }
+ virtual ~IRegion() {
+
+ }
+
+ virtual std::string to_string(RefVector ref)=0;
+ virtual long overlap(IRegion * tmp) =0;
+ position_str get_coordinates() {
+ return this->start;
+ }
+ virtual int support()=0;
+};
+
+#endif /* SUB_IREGION_H_ */
diff --git a/src/tree/BinTree.cpp b/src/tree/BinTree.cpp
new file mode 100644
index 0000000..70ac85d
--- /dev/null
+++ b/src/tree/BinTree.cpp
@@ -0,0 +1,293 @@
+/*
+ * BinTree.cpp
+ *
+ * Created on: Sep 3, 2015
+ * Author: fsedlaze
+ */
+
+#include "BinTree.h"
+
+void BinTree::find(int item, tree_node **par, tree_node **loc) {
+ tree_node *ptr, *ptrsave;
+ if (root == NULL) {
+ *loc = NULL;
+ *par = NULL;
+ return;
+ }
+
+ if (item == root->key) {
+ *loc = root;
+ *par = NULL;
+ return;
+ }
+
+ if (item < root->key) {
+ ptr = root->left;
+ } else {
+ ptr = root->right;
+ }
+ ptrsave = root;
+ while (ptr != NULL) {
+ if (item == ptr->key) {
+ *loc = ptr;
+ *par = ptrsave;
+ return;
+ }
+ ptrsave = ptr;
+ if (item < ptr->key) {
+ ptr = ptr->left;
+ } else {
+ ptr = ptr->right;
+ }
+ }
+ *loc = NULL;
+ *par = ptrsave;
+}
+
+/*
+
+ * Inserting Element into the Tree
+
+ */
+
+void BinTree::insert(tree_node *tree, int value) {
+ if (root == NULL) {
+ root = new tree_node;
+ root->key = value;
+ root->num = 1;
+ root->left = NULL;
+ root->right = NULL;
+ std::cout << "Root tree_node is Added" << std::endl;
+ return;
+ }
+
+ if (tree->key > value) {
+ if (tree->left != NULL) {
+ insert(tree->left, value);
+ } else {
+
+ tree->left = new tree_node;
+ tree->left->key = value;
+ tree->left->num = 1;
+ (tree->left)->left = NULL;
+ (tree->left)->right = NULL;
+ std::cout << "tree_node Added To Left" << std::endl;
+ return;
+ }
+ } else if (tree->key < value) {
+ if (tree->right != NULL) {
+ insert(tree->right, value);
+ } else {
+ tree->right = new tree_node;
+ tree->right->key = value;
+ tree->right->num = 1;
+ (tree->right)->left = NULL;
+ (tree->right)->right = NULL;
+ std::cout << "tree_node Added To Right" << std::endl;
+ return;
+ }
+ } else { // found element -> already exist!
+ tree->num++;
+ }
+}
+
+/*
+
+ * Delete Element from the tree
+
+ */
+
+void BinTree::del(int key) {
+ tree_node *parent, *location;
+
+ if (root == NULL) {
+ std::cout << "Tree empty" << std::endl;
+ return;
+ }
+
+ find(key, &parent, &location);
+ if (location == NULL) {
+ std::cout << "Item not present in tree" << std::endl;
+ return;
+ }
+
+ if (location->left == NULL && location->right == NULL) {
+ case_a(parent, location);
+ }
+ if (location->left != NULL && location->right == NULL) {
+ case_b(parent, location);
+ }
+ if (location->left == NULL && location->right != NULL) {
+ case_b(parent, location);
+ }
+ if (location->left != NULL && location->right != NULL) {
+ case_c(parent, location);
+ }
+ delete location;
+}
+
+/*
+
+ * Case A
+
+ */
+
+void BinTree::case_a(tree_node *par, tree_node *loc) {
+ if (par == NULL) {
+ root = NULL;
+ } else {
+ if (loc == par->left) {
+ par->left = NULL;
+ } else {
+ par->right = NULL;
+ }
+ }
+}
+
+/*
+
+ * Case B
+
+ */
+
+void BinTree::case_b(tree_node *par, tree_node *loc) {
+ tree_node *child;
+ if (loc->left != NULL) {
+ child = loc->left;
+ } else {
+ child = loc->right;
+ }
+ if (par == NULL) {
+ root = child;
+ } else {
+ if (loc == par->left) {
+ par->left = child;
+ } else {
+ par->right = child;
+ }
+ }
+}
+
+/*
+
+ * Case C
+
+ */
+
+void BinTree::case_c(tree_node *par, tree_node *loc) {
+ tree_node *ptr, *ptrsave, *suc, *parsuc;
+ ptrsave = loc;
+ ptr = loc->right;
+ while (ptr->left != NULL) {
+ ptrsave = ptr;
+ ptr = ptr->left;
+ }
+ suc = ptr;
+ parsuc = ptrsave;
+ if (suc->left == NULL && suc->right == NULL) {
+ case_a(parsuc, suc);
+ } else {
+ case_b(parsuc, suc);
+ }
+
+ if (par == NULL) {
+ root = suc;
+ } else {
+ if (loc == par->left) {
+ par->left = suc;
+ } else {
+ par->right = suc;
+ }
+ }
+ suc->left = loc->left;
+ suc->right = loc->right;
+}
+
+void BinTree::get_nodes(tree_node *ptr, std::vector<int> & nodes) {
+ std::cout<<"get_nodes"<<std::endl;
+ if (root == NULL) {
+ std::cout << "Tree is empty" << std::endl;
+ return;
+ }
+ if (ptr != NULL) {
+ nodes.push_back(ptr->key);
+ get_nodes(ptr->left,nodes);
+ get_nodes(ptr->right,nodes);
+ }
+}
+
+/*
+
+ * Pre Order Traversal
+
+ */
+
+void BinTree::preorder(tree_node *ptr) {
+ if (root == NULL) {
+ std::cout << "Tree is empty" << std::endl;
+ return;
+ }
+ if (ptr != NULL) {
+ std::cout << ptr->key << " ";
+ preorder(ptr->left);
+ preorder(ptr->right);
+ }
+}
+
+/*
+
+ * In Order Traversal
+
+ */
+
+void BinTree::inorder(tree_node *ptr) {
+ if (root == NULL) {
+ std::cout << "Tree is empty" << std::endl;
+ return;
+ }
+ if (ptr != NULL) {
+ inorder(ptr->left);
+ std::cout << ptr->key << " ";
+ inorder(ptr->right);
+ }
+}
+
+/*
+
+ * Postorder Traversal
+
+ */
+
+void BinTree::postorder(tree_node *ptr) {
+ if (root == NULL) {
+ std::cout << "Tree is empty" << std::endl;
+ return;
+ }
+ if (ptr != NULL) {
+ postorder(ptr->left);
+ postorder(ptr->right);
+ std::cout << ptr->key << " ";
+ }
+}
+
+/*
+
+ * Display Tree Structure
+
+ */
+void BinTree::display(tree_node *ptr, int level) {
+ int i;
+ if (ptr != NULL) {
+ display(ptr->right, level + 1);
+ std::cout << std::endl;
+ if (ptr == root)
+ std::cout << "Root->: ";
+ else {
+ for (i = 0; i < level; i++) {
+ std::cout << " ";
+ }
+ }
+ std::cout << ptr->key;
+ display(ptr->left, level + 1);
+ }
+}
diff --git a/src/tree/BinTree.h b/src/tree/BinTree.h
new file mode 100644
index 0000000..ca71751
--- /dev/null
+++ b/src/tree/BinTree.h
@@ -0,0 +1,42 @@
+/*
+ * BinTree.h
+ *
+ * Created on: Sep 3, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_BINTREE_H_
+#define TREE_BINTREE_H_
+struct tree_node {
+ int key; // value to store!
+ int num;//times hit 1-> unique
+ struct tree_node *left;
+ struct tree_node *right;
+};
+#include <vector>
+#include <cstddef>
+#include <iostream>
+class BinTree {
+private:
+ tree_node *root;
+public:
+ BinTree() {
+ root = NULL;
+ }
+ ~BinTree(){
+
+ }
+ void find(int item, tree_node **par, tree_node **loc);
+ void insert(tree_node *tree, int value);
+ void del(int key);
+ void case_a(tree_node *par, tree_node *loc);
+ void case_b(tree_node *par, tree_node *loc);
+ void case_c(tree_node *par, tree_node *loc);
+ void preorder(tree_node *ptr);
+ void inorder(tree_node *ptr);
+ void postorder(tree_node *ptr);
+ void display(tree_node *ptr, int);
+ void get_nodes(tree_node *ptr, std::vector<int> & nodes);
+};
+
+#endif /* TREE_BINTREE_H_ */
diff --git a/src/tree/IntervallTree.cpp b/src/tree/IntervallTree.cpp
new file mode 100644
index 0000000..adf8477
--- /dev/null
+++ b/src/tree/IntervallTree.cpp
@@ -0,0 +1,269 @@
+/*
+ * IntervallTree.cpp
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#include "IntervallTree.h"
+
+// Inserting a node
+void IntervallTree::insert(Breakpoint * new_break, TNode *&p) {
+ if (p == NULL) {
+ p = new TNode(new_break); //TODO: this is going to be a problem!
+ if (p == NULL) {
+ std::cout << "Out of Space\n" << std::endl;
+ }
+ } else {
+ long score = p->get_data()->overlap(new_break);
+
+ if (score > 0) {
+ insert(new_break, p->left);
+ if ((bsheight(p->left) - bsheight(p->right)) == 2) {
+ score = p->left->get_data()->overlap(new_break);
+ if (score > 0) {
+ p = srl(p);
+ } else {
+ p = drl(p);
+ }
+ }
+ } else if (score < 0) {
+ insert(new_break, p->right);
+ if ((bsheight(p->right) - bsheight(p->left)) == 2) {
+ score = p->right->get_data()->overlap(new_break);
+ if (score < 0) {
+ p = srr(p);
+ } else {
+ p = drr(p);
+ }
+ }
+ } else { //overlaps!
+ p->get_data()->add_read(new_break);
+ delete new_break;
+ //std::cout << "Breakpoint was detected\n" << std::endl;
+ }
+ }
+ int m, n, d;
+ m = bsheight(p->left);
+ n = bsheight(p->right);
+ d = max(m, n);
+ p->set_height(d + 1);
+}
+// Finding the Smallest
+TNode * IntervallTree::findmin(TNode * p) {
+ if (p == NULL) {
+ std::cout << "The tree is empty\n" << std::endl;
+ return p;
+ } else {
+ while (p->left != NULL) {
+ p = p->left;
+ //return p;
+ }
+ return p;
+ }
+}
+// Finding the Largest node
+TNode * IntervallTree::findmax(TNode * p) {
+ if (p == NULL) {
+ std::cout << "The tree is empty\n" << std::endl;
+ return p;
+ } else {
+ while (p->right != NULL) {
+ p = p->right;
+ //return p;
+ }
+ return p;
+ }
+}
+// Finding an get_value()
+void IntervallTree::find(Breakpoint * point, TNode * &p) {
+ if (p == NULL) {
+ std::cout << "Sorry! get_value() not found\n" << std::endl;
+ } else {
+ long score = p->get_data()->overlap(point);
+ if (score > 0) {
+ find(point, p->left);
+ } else if (score < 0) {
+ find(point, p->right);
+ } else {
+ std::cout << "get_value() found!\n" << std::endl;
+ }
+
+ }
+}
+// Copy a tree
+void IntervallTree::copy(TNode * &p, TNode * &p1) {
+ makeempty(p1);
+ p1 = nodecopy(p);
+}
+// Make a tree empty
+void IntervallTree::makeempty(TNode * &p) {
+ TNode * d;
+ if (p != NULL) {
+ makeempty(p->left);
+ makeempty(p->right);
+ d = p;
+ free(d);
+ p = NULL;
+ }
+}
+// Copy the nodes
+TNode * IntervallTree::nodecopy(TNode * &p) {
+ TNode * temp;
+ if (p == NULL) {
+ return p;
+ } else {
+ temp = new TNode(p->get_data()); //TODO!
+ temp->left = nodecopy(p->left);
+ temp->right = nodecopy(p->right);
+ return temp;
+ }
+}
+
+// Deleting a node
+void IntervallTree::del(Breakpoint * point, TNode * &p) {
+ TNode * d;
+ if (p == NULL) {
+ std::cout << "Sorry! get_value() not found\n" << std::endl;
+ } else {
+ long score = p->get_data()->overlap(point);
+ if (score > 0) {
+ del(point, p->left);
+ } else if (score < 0) {
+ del(point, p->right);
+ } else if ((p->left == NULL) && (p->right == NULL)) {
+ d = p;
+ free(d);
+ p = NULL;
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else if (p->left == NULL) {
+ d = p;
+ free(d);
+ p = p->right;
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else if (p->right == NULL) {
+ d = p;
+ p = p->left;
+ free(d);
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else {
+ //p->set_value(deletemin(p->right));
+ }
+ }
+}
+
+int IntervallTree::deletemin(TNode * &p) {
+ int c;
+ std::cout << "inside deltemin\n" << std::endl;
+ if (p->left == NULL) {
+ //c = p->get_value();
+ p = p->right;
+ return c;
+ } else {
+ c = deletemin(p->left);
+ return c;
+ }
+}
+void IntervallTree::preorder(TNode * p) {
+ if (p != NULL) {
+ //std::cout << p->get_data()->to_string() << "\t";
+ preorder(p->left);
+ preorder(p->right);
+ }
+}
+void IntervallTree::get_breakpoints(TNode *p, std::vector<Breakpoint *> & points) {
+ if (p != NULL) {
+ get_breakpoints(p->left, points);
+ points.push_back(p->get_data());
+ get_breakpoints(p->right, points);
+ }
+}
+
+// Inorder Printing
+void IntervallTree::inorder(TNode * p, TNode * root) {
+ if (p != NULL) {
+ inorder(p->left, root);
+ //std::cout << p->get_data()->to_string();
+ if (p == root) {
+ std::cout << "*\t";
+ } else {
+ std::cout << "\t";
+ }
+ inorder(p->right, root);
+ }
+}
+
+// PostOrder Printing
+void IntervallTree::postorder(TNode * p) {
+ if (p != NULL) {
+ postorder(p->left);
+ postorder(p->right);
+ //std::cout << p->get_data()->to_string() << "\t";
+ }
+}
+
+int IntervallTree::max(int value1, int value2) {
+ return ((value1 > value2) ? value1 : value2);
+}
+int IntervallTree::bsheight(TNode * p) {
+ int t;
+ if (p == NULL) {
+ return -1;
+ } else {
+ t = p->get_height();
+ return t;
+ }
+}
+
+TNode * IntervallTree::srl(TNode * &p1) {
+ TNode * p2;
+ p2 = p1->left;
+ p1->left = p2->right;
+ p2->right = p1;
+ p1->set_height(max(bsheight(p1->left), bsheight(p1->right)) + 1);
+ p2->set_height(max(bsheight(p2->left), p1->get_height()) + 1);
+ return p2;
+}
+TNode * IntervallTree::srr(TNode * &p1) {
+ TNode * p2;
+ p2 = p1->right;
+ p1->right = p2->left;
+ p2->left = p1;
+ p1->set_height(max(bsheight(p1->left), bsheight(p1->right)) + 1);
+ p2->set_height(max(p1->get_height(), bsheight(p2->right)) + 1);
+ return p2;
+}
+TNode * IntervallTree::drl(TNode * &p1) {
+ p1->left = srr(p1->left);
+ return srl(p1);
+}
+TNode * IntervallTree::drr(TNode * &p1) {
+ p1->right = srl(p1->right);
+ return srr(p1);
+}
+
+int IntervallTree::nonodes(TNode * p) {
+ int count = 0;
+ if (p != NULL) {
+ nonodes(p->left);
+ nonodes(p->right);
+ count++;
+ }
+ return count;
+}
+
+void IntervallTree::collapse_intervalls(TNode *&p) {
+ std::cout << "\t Collapse" << std::endl;
+ TNode * new_root = NULL;
+ std::vector<Breakpoint *> points;
+ get_breakpoints(p, points);
+
+ for (size_t i = 0; i < points.size(); i++) {
+ if (points[i]->get_support() > Parameter::Instance()->min_support) {
+ //std::cout << "\tpoints: " << points[i]->to_string(ref) << std::endl;
+ this->insert(points[i], new_root);
+ }
+ }
+ this->makeempty(p);
+ p = new_root;
+}
diff --git a/src/tree/IntervallTree.h b/src/tree/IntervallTree.h
new file mode 100644
index 0000000..9b5b57f
--- /dev/null
+++ b/src/tree/IntervallTree.h
@@ -0,0 +1,40 @@
+/*
+ * IntervallTree.h
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_INTERVALLTREE_H_
+#define TREE_INTERVALLTREE_H_
+
+#include <vector>
+
+#include "TNode.h"
+class IntervallTree {
+private:
+ int max(int, int);
+ TNode * srl(TNode *&);
+ TNode * drl(TNode *&);
+ TNode * srr(TNode *&);
+ TNode * drr(TNode *&);
+public:
+ void insert(Breakpoint * point, TNode *&);
+ void del(Breakpoint * point, TNode *&);
+ int deletemin(TNode *&);
+ void find(Breakpoint * point, TNode *&);
+ TNode * findmin(TNode*);
+ TNode * findmax(TNode*);
+ void makeempty(TNode *&);
+ void copy(TNode * &, TNode *&);
+ TNode * nodecopy(TNode *&);
+ void preorder(TNode*);
+ void inorder(TNode*,TNode * root);
+ void postorder(TNode*);
+ int bsheight(TNode*);
+ void get_breakpoints(TNode *p,std::vector<Breakpoint *> & points);
+ int nonodes(TNode*);
+ void collapse_intervalls(TNode *&p);
+};
+
+#endif /* TREE_INTERVALLTREE_H_ */
diff --git a/src/tree/Intervall_bed.cpp b/src/tree/Intervall_bed.cpp
new file mode 100644
index 0000000..58ae010
--- /dev/null
+++ b/src/tree/Intervall_bed.cpp
@@ -0,0 +1,245 @@
+/*
+ * Intervall_bed.cpp
+ *
+ * Created on: Feb 4, 2016
+ * Author: fsedlaze
+ */
+
+#include "Intervall_bed.h"
+
+
+// Inserting a node
+void IntervallTree_bed::insert(long start, long stop, Leaf *&p) {
+ if (p == NULL) {
+ p = new Leaf(start, stop);
+ if (p == NULL) {
+ std::cout << "Out of Space\n" << std::endl;
+ }
+ } else {
+
+ long score = p->overlap(start, stop);
+
+ if (score > 0) {
+ insert(start, stop, p->left);
+ if ((bsheight(p->left) - bsheight(p->right)) == 2) {
+ score = p->left->overlap(start, stop);
+ if (score > 0) {
+ p = srl(p);
+ } else {
+ p = drl(p);
+ }
+ }
+ } else if (score < 0) {
+ insert(start, stop, p->right);
+ if ((bsheight(p->right) - bsheight(p->left)) == 2) {
+ score = p->right->overlap(start, stop);
+ if (score < 0) {
+ p = srr(p);
+ } else {
+ p = drr(p);
+ }
+ }
+ } else { //overlaps!
+ std::cerr << "Two regions overlap and are thus ignored:" << std::endl;
+ }
+ }
+ int m, n, d;
+ m = bsheight(p->left);
+ n = bsheight(p->right);
+ d = max(m, n);
+ p->set_height(d + 1);
+}
+// Finding the Smallest
+Leaf * IntervallTree_bed::findmin(Leaf * p) {
+ if (p == NULL) {
+ return p;
+ } else {
+ while (p->left != NULL) {
+ p = p->left;
+ //return p;
+ }
+ return p;
+ }
+}
+// Finding the Largest node
+Leaf * IntervallTree_bed::findmax(Leaf * p) {
+ if (p == NULL) {
+ return p;
+ } else {
+ while (p->right != NULL) {
+ p = p->right;
+ //return p;
+ }
+ return p;
+ }
+}
+// Finding an get_value()
+bool IntervallTree_bed::is_in(long position, Leaf * &p) {
+ if (p == NULL) {
+ return false;
+ } else {
+ long score = p->overlap(position);
+ if (score > 0) {
+ is_in(position, p->left);
+ } else if (score < 0) {
+ is_in(position, p->right);
+ } else {
+ return true;
+ }
+
+ }
+}
+// Copy a tree
+void IntervallTree_bed::copy(Leaf * &p, Leaf * &p1) {
+ makeempty(p1);
+ p1 = nodecopy(p);
+}
+// Make a tree empty
+void IntervallTree_bed::makeempty(Leaf * &p) {
+ Leaf * d;
+ if (p != NULL) {
+ makeempty(p->left);
+ makeempty(p->right);
+ d = p;
+ free(d);
+ p = NULL;
+ }
+}
+// Copy the nodes
+Leaf * IntervallTree_bed::nodecopy(Leaf * &p) {
+ Leaf * temp;
+ if (p == NULL) {
+ return p;
+ } else {
+ temp = new Leaf(p->get_start(), p->get_stop()); //TODO!
+ temp->left = nodecopy(p->left);
+ temp->right = nodecopy(p->right);
+ return temp;
+ }
+}
+
+// Deleting a node
+void IntervallTree_bed::del(long start, long stop, Leaf * &p) {
+ Leaf * d;
+ if (p == NULL) {
+ std::cout << "Sorry! get_value() not found\n" << std::endl;
+ } else {
+ long score = p->overlap(start, stop);
+ if (score > 0) {
+ del(start, stop, p->left);
+ } else if (score < 0) {
+ del(start, stop, p->right);
+ } else if ((p->left == NULL) && (p->right == NULL)) {
+ d = p;
+ free(d);
+ p = NULL;
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else if (p->left == NULL) {
+ d = p;
+ free(d);
+ p = p->right;
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else if (p->right == NULL) {
+ d = p;
+ p = p->left;
+ free(d);
+ std::cout << "get_value() deleted successfully\n" << std::endl;
+ } else {
+ //p->set_value(deletemin(p->right));
+ }
+ }
+}
+
+int IntervallTree_bed::deletemin(Leaf * &p) {
+ int c;
+ std::cout << "inside deltemin\n" << std::endl;
+ if (p->left == NULL) {
+ //c = p->get_value();
+ p = p->right;
+ return c;
+ } else {
+ c = deletemin(p->left);
+ return c;
+ }
+}
+void IntervallTree_bed::preorder(Leaf * p) {
+ if (p != NULL) {
+ //std::cout << p->get_data()->to_string() << "\t";
+ preorder(p->left);
+ preorder(p->right);
+ }
+}
+
+// Inorder Printing
+void IntervallTree_bed::inorder(Leaf * p, Leaf * root) {
+ if (p != NULL) {
+ inorder(p->left, root);
+ //std::cout << p->get_data()->to_string();
+ if (p == root) {
+ std::cout << "*\t";
+ } else {
+ std::cout << "\t";
+ }
+ inorder(p->right, root);
+ }
+}
+
+// PostOrder Printing
+void IntervallTree_bed::postorder(Leaf * p) {
+ if (p != NULL) {
+ postorder(p->left);
+ postorder(p->right);
+ std::cout << p->get_start()<<" "<<p->get_stop()<< "\t";
+ }
+}
+
+int IntervallTree_bed::max(int value1, int value2) {
+ return ((value1 > value2) ? value1 : value2);
+}
+int IntervallTree_bed::bsheight(Leaf * p) {
+ int t;
+ if (p == NULL) {
+ return -1;
+ } else {
+ t = p->get_height();
+ return t;
+ }
+}
+
+Leaf * IntervallTree_bed::srl(Leaf * &p1) {
+ Leaf * p2;
+ p2 = p1->left;
+ p1->left = p2->right;
+ p2->right = p1;
+ p1->set_height(max(bsheight(p1->left), bsheight(p1->right)) + 1);
+ p2->set_height(max(bsheight(p2->left), p1->get_height()) + 1);
+ return p2;
+}
+Leaf * IntervallTree_bed::srr(Leaf * &p1) {
+ Leaf * p2;
+ p2 = p1->right;
+ p1->right = p2->left;
+ p2->left = p1;
+ p1->set_height(max(bsheight(p1->left), bsheight(p1->right)) + 1);
+ p2->set_height(max(p1->get_height(), bsheight(p2->right)) + 1);
+ return p2;
+}
+Leaf * IntervallTree_bed::drl(Leaf * &p1) {
+ p1->left = srr(p1->left);
+ return srl(p1);
+}
+Leaf * IntervallTree_bed::drr(Leaf * &p1) {
+ p1->right = srl(p1->right);
+ return srr(p1);
+}
+
+int IntervallTree_bed::nonodes(Leaf * p) {
+ int count = 0;
+ if (p != NULL) {
+ nonodes(p->left);
+ nonodes(p->right);
+ count++;
+ }
+ return count;
+}
+
diff --git a/src/tree/Intervall_bed.h b/src/tree/Intervall_bed.h
new file mode 100644
index 0000000..83ce14d
--- /dev/null
+++ b/src/tree/Intervall_bed.h
@@ -0,0 +1,39 @@
+/*
+ * Intervall_bed.h
+ *
+ * Created on: Feb 4, 2016
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_INTERVALL_BED_H_
+#define TREE_INTERVALL_BED_H_
+
+#include "Leaf.h"
+#include <vector>
+#include "../Paramer.h"
+
+class IntervallTree_bed {
+private:
+ int max(int, int);
+ Leaf * srl(Leaf *&);
+ Leaf * drl(Leaf *&);
+ Leaf * srr(Leaf *&);
+ Leaf * drr(Leaf *&);
+public:
+ void insert(long start, long stop, Leaf *&);
+ int deletemin(Leaf *&);
+ bool is_in(long pos, Leaf *&); //true if found
+ Leaf * findmin(Leaf*);
+ Leaf * findmax(Leaf*);
+ void makeempty(Leaf *&);
+ void copy(Leaf * &, Leaf *&);
+ Leaf * nodecopy(Leaf *&);
+ void preorder(Leaf*);
+ void inorder(Leaf*, Leaf * root);
+ void postorder(Leaf*);
+ int bsheight(Leaf*);
+ int nonodes(Leaf*);
+ void del(long start, long stop, Leaf * &p);
+};
+
+#endif /* TREE_INTERVALL_BED_H_ */
diff --git a/src/tree/Leaf.h b/src/tree/Leaf.h
new file mode 100644
index 0000000..4a7bd19
--- /dev/null
+++ b/src/tree/Leaf.h
@@ -0,0 +1,64 @@
+/*
+ * Leaf.h
+ *
+ * Created on: Feb 4, 2016
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_LEAF_H_
+#define TREE_LEAF_H_
+#include "../Paramer.h"
+#include <iostream>
+#include <cstdlib>
+#include <cmath>
+class Leaf {
+private:
+ long start;
+ long stop;
+ int height;
+ void init() {
+ height = 0;
+ this->parent = NULL;
+ this->left = NULL;
+ this->right = NULL;
+ }
+
+public:
+ Leaf * parent;
+ Leaf * left;
+ Leaf * right;
+
+ Leaf(long start, long stop) {
+ this->start = start;
+ this->stop = stop;
+ init();
+ }
+ int overlap(long position) {
+ if (abs(position - get_start()) < Parameter::Instance()->max_dist && abs(position - get_stop()) < Parameter::Instance()->max_dist) {
+ return 0;
+ }
+ return (position - start); //((start < position) && (stop > position));
+ }
+ int overlap(long start, long stop) {
+ if (abs(start - get_start()) < Parameter::Instance()->max_dist && abs(stop - get_stop()) < Parameter::Instance()->max_dist) {
+ return 0;
+ }
+ //as abstraction lets try the start+stop coordinate!
+ return (start - get_start()); // + (stop-get_stop());
+ }
+
+ long get_start() {
+ return start;
+ }
+ long get_stop() {
+ return stop;
+ }
+ int get_height() {
+ return height;
+ }
+ void set_height(int val) {
+ this->height = val;
+ }
+};
+
+#endif /* TREE_LEAF_H_ */
diff --git a/src/tree/Scapegoat.cpp b/src/tree/Scapegoat.cpp
new file mode 100644
index 0000000..2b6f9dd
--- /dev/null
+++ b/src/tree/Scapegoat.cpp
@@ -0,0 +1,10 @@
+/*
+ * Scapegoat.cpp
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#include "Scapegoat.h"
+
+
diff --git a/src/tree/Scapegoat.h b/src/tree/Scapegoat.h
new file mode 100644
index 0000000..94a625c
--- /dev/null
+++ b/src/tree/Scapegoat.h
@@ -0,0 +1,234 @@
+/*
+ * Scapegoat.h
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_SCAPEGOAT_H_
+#define TREE_SCAPEGOAT_H_
+#include "TNode.h"
+/*
+struct TNode {
+ TNode * parent;
+ TNode * left;
+ TNode * right;
+ Breakpoint * data;
+ int value;
+};*/
+
+/*
+class Scapegoat {
+private:
+ TNode * root;
+ int n, q;
+public:
+ Scapegoat() {
+ root = NULL;
+ n = 0;
+ }
+
+ bool isEmpty() {
+ return root == NULL;
+ }
+
+ void makeEmpty() {
+ root = NULL;
+ n = 0;
+ }
+
+ int size(TNode *r) {
+ if (r == NULL)
+ return 0;
+ else {
+ int l = 1;
+ l += size(r->left);
+ l += size(r->right);
+ return l;
+ }
+ }
+
+ bool search(int value, Breakpoint * point) {
+ return search(root,value, point);
+ }
+ bool search(TNode *r, int value, Breakpoint * point) {
+ bool found = false;
+ while ((r != NULL) && !found) {
+
+ //int score=r->get_data()->overlap(point);
+ //if (score>0) {
+
+ if(r->get_value() > value){
+ r = r->left;
+ // } else if (score<0) {
+ }else if (r->get_value() < value){
+ r = r->right;
+ } else {
+ found = true;
+ break;
+ }
+ found = search(r, value,point);
+ }
+ return found;
+ }
+
+ int size() {
+ return n;
+ }
+ void inorder() {
+ inorder(root);
+ }
+
+ void inorder(TNode *r) {
+ if (r != NULL) {
+ inorder(r->left);
+ std::cout << r->get_value() ;//->get_data()->to_string() << " ";
+ if(r==this->root){
+ std::cout<<"* ";
+ }else{
+ std::cout<<" ";
+ }
+ inorder(r->right);
+ }
+ }
+
+
+ void preorder() {
+ preorder(root);
+ }
+
+ void preorder(TNode *r) {
+
+ if (r != NULL) {
+ std::cout << r->get_value() <<" ";//get_data()->to_string() << " ";
+ preorder(r->left);
+ preorder(r->right);
+ }
+
+ }
+
+ void postorder() {
+ postorder(root);
+ }
+
+ void postorder(TNode *r) {
+ if (r != NULL) {
+ postorder(r->left);
+ postorder(r->right);
+ std::cout << r->get_value() <<" ";//->get_data()->to_string() << " ";
+ }
+ }
+
+ static int const log32(int q) {
+ double const log23 = 2.4663034623764317;
+ return (int) ceil(log23 * log(q));
+ }
+
+ bool add(int value,Breakpoint * point) {
+ TNode *u = new TNode(value,point);
+ int d = addWithDepth(u);
+ if (d > log32(q)) {
+ TNode *w = u->parent;
+ while (3 * size(w) <= 2 * size(w->parent)) {
+ w = w->parent;
+ }
+ rebuild(w->parent);
+ }
+ return d >= 0;
+ }
+
+
+ void rebuild(TNode *u) {
+ int ns = size(u);
+ TNode *p = u->parent;
+ TNode **a = new TNode*[ns];
+ packIntoArray(u, a, 0);
+ if (p == NULL) {
+ root = buildBalanced(a, 0, ns);
+ root->parent = NULL;
+ } else if (p->right == u) {
+ p->right = buildBalanced(a, 0, ns);
+ p->right->parent = p;
+ } else {
+ p->left = buildBalanced(a, 0, ns);
+ p->left->parent = p;
+ }
+ }
+
+
+ int packIntoArray(TNode *u, TNode *a[], int i) {
+
+ if (u == NULL) {
+ return i;
+ }
+ i = packIntoArray(u->left, a, i);
+ a[i++] = u;
+ return packIntoArray(u->right, a, i);
+ }
+
+
+ TNode *buildBalanced(TNode **a, int i, int ns)
+
+ {
+ if (ns == 0) {
+ return NULL;
+ }
+ int m = ns / 2;
+ a[i + m]->left = buildBalanced(a, i, m);
+ if (a[i + m]->left != NULL) {
+ a[i + m]->left->parent = a[i + m];
+ }
+ a[i + m]->right = buildBalanced(a, i + m + 1, ns - m - 1);\
+ if (a[i + m]->right != NULL) {
+ a[i + m]->right->parent = a[i + m];
+ }
+ return a[i + m];
+ }
+
+
+ int addWithDepth(TNode *u) {
+ TNode *w = root;
+ if (w == NULL) {
+ root = u;
+ n++;
+ q++;
+ return 0;
+ }
+ bool done = false;
+ int d = 0;
+ do {
+ //int score=u->get_data()->overlap(w->get_data());
+ //if (score>0) {
+
+ if(w->get_value() > u->get_value()){
+ if (w->left == NULL) {
+ w->left = u;
+ u->parent = w;
+ done = true;
+ } else {
+ w = w->left;
+ }
+ //} else if (score<0) {
+ }else if(w->get_value() < u->get_value()){
+ if (w->right == NULL) {
+ w->right = u;
+ u->parent = w;
+ done = true;
+ } else {
+
+ w = w->right;
+ }
+ } else {
+ //are equal!
+ return -1;
+ }
+ d++;
+ } while (!done);
+ n++;
+ q++;
+ return d;
+ }
+
+};
+*/
+#endif /* TREE_SCAPEGOAT_H_ */
diff --git a/src/tree/TNode.h b/src/tree/TNode.h
new file mode 100644
index 0000000..7183d25
--- /dev/null
+++ b/src/tree/TNode.h
@@ -0,0 +1,68 @@
+/*
+ * TNode.h
+ *
+ * Created on: Jun 23, 2015
+ * Author: fsedlaze
+ */
+
+#ifndef TREE_TNODE_H_
+#define TREE_TNODE_H_
+
+#include <iostream>
+#include <cstdlib>
+#include <cmath>
+#include "../sub/Breakpoint.h"
+//#include "TNode.h"
+class TNode {
+private:
+
+ Breakpoint * data;
+ //int value;
+ int height;
+ int MAX_DIST;
+ void init() {
+ this->parent = NULL;
+ this->left = NULL;
+ this->right = NULL;
+ MAX_DIST=500;
+ }
+public:
+ TNode * parent;
+ TNode * left;
+ TNode * right;
+ TNode() {
+ height=0;
+ init();
+ this->data=NULL;
+ }
+ TNode(Breakpoint * point) {
+ init();
+ this->data=point;
+ height=0;
+ }
+
+ ~TNode() {
+
+ }
+
+ Breakpoint * get_data() {
+ return data;
+ }
+ int get_height() {
+ return height;
+ }
+ void set_height(int val) {
+ this->height = val;
+ }
+ /*int overlap(TNode * tmp) {
+ int score = 0; //(tmp->get_type() - type); // 0 if both points are on the same chr!
+ if (abs(tmp->get_data()->get_coordinates().start - data->get_coordinates().start) < MAX_DIST
+ && abs(tmp->get_data()->get_coordinates().stop - data->get_coordinates().stop) < MAX_DIST) {
+ return score;
+ }
+ //as abstraction lets try the start coordinate only!
+ return abs(tmp->get_data()->get_coordinates().start - data->get_coordinates().start);
+ }*/
+};
+
+#endif /* TREE_TNODE_H_ */
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/sniffles.git
More information about the debian-med-commit
mailing list