Commit cd703265 authored by Afif Elghraoui's avatar Afif Elghraoui
Browse files

New upstream version 1.0.8+ds

parent 6b4aed99

CMakeLists.txt

+4 −2
Original line number Diff line number Diff line 
cmake_minimum_required(VERSION 2.8)

project(Sniffles)



option(STATIC "Build static binary" OFF)



 set( SNIF_VERSION_MAJOR 1 )

 set( SNIF_VERSION_MINOR 0 )

IF(CMAKE_BUILD_TYPE STREQUAL "Debug")

    message(STATUS "Building in debug mode!")

    set( SNIF_VERSION_BUILD 7-debug )

    set( SNIF_VERSION_BUILD 8-debug )

else()

	set( SNIF_VERSION_BUILD 7 )

	set( SNIF_VERSION_BUILD 8 )

ENDIF()

README.md

+5 −8
Original line number Diff line number Diff line 
# Sniffles

Sniffles is a structural variation caller using third generation sequencing (PacBio or Oxford Nanopore). It detects all types of SVs (10bp+) using evidence from split-read alignments, high-mismatch regions, and coverage analysis. Please note the current version of Sniffles requires sorted output from BWA-MEM (use -M and -x parameter) or NGM-LR with the optional SAM attributes enabled! If you experience problems or have suggestions please contact: fritz.sedlazeck@gmail.com

Sniffles is a structural variation caller using third generation sequencing (PacBio or Oxford Nanopore). It detects all types of SVs (10bp+) using evidence from split-read alignments, high-mismatch regions, and coverage analysis. Please note the current version of Sniffles requires sorted output from BWA-MEM (use -M and -x parameter) or NGMLR with the optional SAM attributes enabled! If you experience problems or have suggestions please contact: fritz.sedlazeck@gmail.com





Please see our github wiki for more information (https://github.com/fritzsedlazeck/Sniffles/wiki)



**************************************



## NextGenMap-LR: (NGM-LR)





Sniffles performs best with the mappings of NGM-LR our novel long read mapping method. 

## NGMLR

Sniffles performs best with the mappings of NGMLR our novel long read mapping method. 

Please see:

https://github.com/philres/nextgenmap-lr

https://github.com/philres/ngmlr



****************************************

## Citation:
@@ -24,7 +21,7 @@ http://www.biorxiv.org/content/early/2017/07/28/169557 
[Accurate and fast detection of complex and nested structural variations using long read technologies](http://schatzlab.cshl.edu/presentations/2016/2016.10.28.BIODATA.PacBioSV.pdf)

Biological Data Science, Cold Spring Harbor Laboratory, Cold Spring Harbor, NY, 26 - 29.10.2016



[NextGenMap-LR: Highly accurate read mapping of third generation sequencing reads for improved structural variation analysis](http://www.cibiv.at/~philipp_/files/gi2016_poster_phr.pdf) 

[NGMLR: Highly accurate read mapping of third generation sequencing reads for improved structural variation analysis](http://www.cibiv.at/~philipp_/files/gi2016_poster_phr.pdf) 

Genome Informatics 2016, Wellcome Genome Campus Conference Centre, Hinxton, Cambridge, UK, 19.09.-2.09.2016



**************************************

src/Alignment.cpp

+132 −22
Original line number Diff line number Diff line
@@ -51,6 +51,81 @@ void add_event(int pos, size_t & i, vector<differences_str> & events) { 
	events.insert(events.begin() + i, ev);

}



/*

 * :6-ata:10+gtc:4*at:3,

 * -ata: del

 * +gtc: ins

 * *at a->t

 * :10 match bp

 * cs:Z::5+g:7+t:4+c:7-a:6*tg:15*tg:2+g:6+c:1+c:1+c:4+a:15+g:4+c:12*gt*ta:2+c:7+t:8+t:17+t:8+g:5+g:10-g:4+g*tc:5+t*ag:10+a:4-ttt

 */



vector<differences_str> Alignment::summarize_csstring(std::vector<indel_str> &dels) {

	string cs = this->get_cs();

	int pos = this->getPosition();

	int corr = 0;

	int ref_pos = 0;

	size_t pos_events = 0;

	int max_size = (this->getRefLength() * 0.9) + getPosition();

	//	comp_aln = clock();



	indel_str del;

	del.sequence = "";

	del.pos = -1;



	vector<differences_str> events;

	differences_str ev;

	std::cout<<"CS: "<<cs<<std::endl;

	for (size_t i = 0; i < cs.size() && pos < max_size; i++) {

		ev.position = pos;

		if (cs[i] == ':') { //match (can be numbers or bp!)

			i++;

			int num = atoi(&cs[i]); //is 0 if letter!

			bool is_long = ((num == 0 && cs[i] != '0'));

			while ((i < cs.size() && cs[i] != '+') && (cs[i] != '-' && cs[i] != '+')) {

				i++;

				if (is_long) {

					num++;

				}

			}

			cout<<"\tMatch: "<<num<<std::endl;

			pos += num;

		} else if (cs[i] == '*') { //mismatch (ref,alt) pairs

			//only every second char counts!

			add_event(pos, pos_events, events);

			pos++;

			i += 2;

			ev.type = 0; //mismatch

			cout<<"\tMiss: "<<1<<std::endl;

		} else if (cs[i] == '-') { //del

			//collet del seq in dels!

			indel_str del;

			del.sequence = "";

			del.pos = pos;

			while ((i < cs.size() && cs[i] != '+') && (cs[i] != '-' && cs[i] != '+')) {

				del.sequence += cs[i];

			}

			dels.push_back(del);

			pos += del.sequence.size();

			ev.type = del.sequence.size();

			cout<<"\tDEL: "<<del.sequence.size()<<std::endl;



		} else if (cs[i] == '+') { //ins

			int num=0;

			while ((i < cs.size() && cs[i] != '+') && (cs[i] != '-' && cs[i] != '+')) {

				i++;

				num++;

			}

			ev.type =num*-1;

			cout<<"\tINS: "<<num<<std::endl;

		}

		events.push_back(ev);



	}

	std::cout<<"end CS: "<<cs<<std::endl;

	return events;

}



vector<differences_str> Alignment::summarizeAlignment(std::vector<indel_str> &dels) {

//	clock_t comp_aln = clock();

	vector<differences_str> events;
@@ -69,16 +144,19 @@ vector<differences_str> Alignment::summarizeAlignment(std::vector<indel_str> &de 
			ev.position = pos;

			ev.type = al->CigarData[i].Length; //deletion

			ev.readposition = read_pos;

			ev.resolved = true;

			events.push_back(ev);

			pos += al->CigarData[i].Length;

		} else if (al->CigarData[i].Type == 'I') {

			ev.position = pos;

			ev.resolved = true;

			ev.readposition = read_pos;

			ev.type = al->CigarData[i].Length * -1; //insertion

			events.push_back(ev);

			read_pos += al->CigarData[i].Length;

		} else if (al->CigarData[i].Type == 'N') {

			pos += al->CigarData[i].Length;

			ev.resolved = true;

			read_pos += al->CigarData[i].Length;

		} else if (al->CigarData[i].Type == 'S' && al->CigarData[i].Length > Parameter::Instance()->huge_ins) { /// Used for reads ranging into an inser

			string sa;
@@ -91,6 +169,7 @@ vector<differences_str> Alignment::summarizeAlignment(std::vector<indel_str> &de 
				} else {

					ev.readposition = read_pos;

				}

				ev.resolved = false;

				ev.type = Parameter::Instance()->huge_ins * -1; //insertion: WE have to fix the length since we cannot estimate it!]

				events.push_back(ev);

			}
@@ -842,9 +921,24 @@ std::string Alignment::get_md() { 
	std::string md;

	if (al->GetTag("MD", md)) {

		return md;

	}else{

		std::cerr<<"No MD string detected! Check bam file! Otherwise generate using e.g. samtools."<<std::endl;

		exit(0);

	}

	return md;

}



std::string Alignment::get_cs() {

	std::string cs;

	if (al->GetTag("cs", cs)) {

		return cs;

	}else{

		std::cerr<<"No CS string detected! Check bam file!"<<std::endl;

		exit(0);

	}

	return cs;

}



vector<str_event> Alignment::get_events_MD(int min_mis) {

	vector<str_event> events;

	/*std::string md;
@@ -1013,7 +1107,13 @@ vector<str_event> Alignment::get_events_Aln() { 


//clock_t comp_aln = clock();

	std::vector<indel_str> dels;

	vector<differences_str> event_aln = summarizeAlignment(dels);

	vector<differences_str> event_aln;

	if (Parameter::Instance()->cs_string) {

		cout<<"run cs check "<<std::endl;

		event_aln = summarize_csstring(dels);

	} else {

		event_aln = summarizeAlignment(dels);

	}

//double time2 = Parameter::Instance()->meassure_time(comp_aln, "\tcompAln Events: ");



	vector<str_event> events;
@@ -1032,7 +1132,7 @@ vector<str_event> Alignment::get_events_Aln() { 
	for (size_t i = 0; i < event_aln.size(); i++) {

		pair_str tmp;

		tmp.position = -1;

		if (event_aln[i].type == 0) {

		if (event_aln[i].type == 0) { //substitutions.

			tmp = plane->add_mut(event_aln[i].position, 1, Parameter::Instance()->window_thresh);

		} else {

			tmp = plane->add_mut(event_aln[i].position, 1, Parameter::Instance()->window_thresh);	// abs(event_aln[i].type)
@@ -1049,7 +1149,6 @@ vector<str_event> Alignment::get_events_Aln() { 
	int stop = 0;

	size_t start = 0;

	for (size_t i = 0; i < profile.size() && stop < event_aln.size(); i++) {



		if (profile[i].position >= event_aln[stop].position) {

			//find the postion:

			size_t pos = 0;
@@ -1072,8 +1171,10 @@ vector<str_event> Alignment::get_events_Aln() { 
				}

				prev += prev_type;

			}

			start++; //we are running one too far!



			if (start + 1 < event_aln.size()) { //TODO do some testing!

				start++; //we are running one too far!

			}

			//run forward to identify the stop:

			prev = event_aln[pos].position;

			stop = pos;
@@ -1091,8 +1192,11 @@ vector<str_event> Alignment::get_events_Aln() { 
			if (stop > 0) {

				stop--;

			}



			//	cout<<start<<" events: "<<event_aln[start].type <<" pos "<<event_aln[start].readposition<<endl;

			int insert_max_pos = 0;

			int insert_max = 0;



			if (event_aln[start].type < 0) {

				insert_max_pos = event_aln[start].position;

				insert_max = abs(event_aln[start].type);
@@ -1101,6 +1205,13 @@ vector<str_event> Alignment::get_events_Aln() { 
			int del_max = 0;

			int del_max_pos = 0;



			if (event_aln[start].type > 0) {

				//	cout<<"HIT"<<endl;

				del_max_pos = event_aln[start].position;

				del_max = event_aln[start].type;



			}



			double insert = 0;

			double del = 0;

			double mismatch = 0;
@@ -1122,6 +1233,8 @@ vector<str_event> Alignment::get_events_Aln() { 
					}

				}

			}



			//	cout << "DELMAX: " << del_max << " " << Parameter::Instance()->avg_del << endl;

			str_event tmp;

			tmp.pos = event_aln[start].position;
@@ -1148,9 +1261,6 @@ vector<str_event> Alignment::get_events_Aln() { 
					}

					tmp.read_pos = event_aln[start].readposition;

					if (Parameter::Instance()->print_seq) {

						//if (tmp.read_pos + tmp.length > this->getAlignment()->QueryBases.size() || tmp.read_pos<0) {

						//	cerr << "BUG! ALN event INS: " << this->getName() << " " << tmp.read_pos << " " << tmp.length << " " << this->getAlignment()->QueryBases.size() << endl;

						//	}

						if (flag) {

							std::cout << "Seq+:" << this->getAlignment()->QueryBases.substr(tmp.read_pos, tmp.length) << std::endl;

src/Alignment.h

+3 −0
Original line number Diff line number Diff line
@@ -40,6 +40,7 @@ struct differences_str{ 
	int position;

	int readposition;

	short type;

	bool resolved;

};

struct indel_str{

	int pos;
@@ -81,6 +82,7 @@ private: 
	 size_t get_length(std::vector<CigarOp> CigarData);

	 int get_id(RefVector ref, std::string chr);

	 vector<differences_str> summarizeAlignment(std::vector<indel_str> &dels);

	 vector<differences_str> summarize_csstring (std::vector<indel_str> &dels) ;

	 void sort_insert(aln_str tmp, vector<aln_str> &entries);



	 void sort_insert_ref(aln_str tmp, vector<aln_str> &entries);
@@ -141,6 +143,7 @@ public: 
	 double get_num_mismatches(std::string md);

	 double get_scrore_ratio();

	 std::string get_md();

	 std::string get_cs();

	 double get_avg_indel_length_Cigar();

	 vector<int> get_avg_diff(double & dist,double & avg_del, double & avg_len);

src/ArgParseOutput.h

0 → 100644
+53 −0
Original line number Diff line number Diff line 
/**

 * Contact: philipp.rescheneder@gmail.com

 */



#include <iostream>

#include <string>



#include <tclap/CmdLine.h>



using std::cerr;

using std::cout;

using std::endl;



class ArgParseOutput : public TCLAP::StdOutput

{

private:



	std::string usageStr;



	std::string versionStr;



public:



	ArgParseOutput(std::string usage, std::string version) {

		usageStr = usage;

		versionStr = version;

	}



	virtual ~ArgParseOutput() {



	}



	virtual void failure(TCLAP::CmdLineInterface& c, TCLAP::ArgException& e) {

		cerr << "Error:" << endl;

		cerr << "         " << e.error() << endl;

		cerr << endl;

		cerr << "Short usage:" << endl;

		cerr << "       sniffles [options] -m <sorted.bam> -v <output.vcf> " << endl;

		cerr << endl;

		cerr << "For complete USAGE and HELP type:" << endl;

		cerr << "    sniffles --help" << endl;

		cerr << endl;

		exit(1);

	}



	virtual void usage(TCLAP::CmdLineInterface& c) {

		cerr << usageStr << std::endl;

	}



	virtual void version(TCLAP::CmdLineInterface& c) {



	}

};