Merge branch 'master' into develop

.github/ISSUE_TEMPLATE.rb

@@ -0,0 +1,22 @@
+# Please report
+
+- [ ] version of ABySS with `abyss-pe version`
+- [ ] distribution of Linux with `lsb_release -d`
+
+# Assembly error
+
+- [ ] complete `abyss-pe` command line
+- [ ] last 20 lines of the output of `abyss-pe`
+- [ ] number of sequenced bases
+- [ ] estimated genome size and ploidy
+- [ ] estimated sequencing depth of coverage
+
+# Build error
+
+Consider installing ABySS using [Linuxbrew](https://linuxbrew.sh) on Linux or [Homebrew](https://brew.sh) on macOS with `brew install abyss`, or using [Bioconda](https://bioconda.github.io) with `conda install abyss`.
+
+- [ ] Have you tried installing ABySS using Brew or Bioconda?
+- [ ] version of GCC or compiler with `gcc --version`
+- [ ] complete `./configure` command line
+- [ ] last 20 lines of the output of `./configure`
+- [ ] last 20 lines of the output of `make`

Sealer/sealer.cc

@@ -111,6 +111,7 @@ static const char USAGE_MESSAGE[] =
 "  -s, --search-mem=N           mem limit for graph searches; multiply by the\n"
 "                               number of threads (-j) to get the total mem used\n"
 "                               for graph traversal [500M]\n"
+"  -g, --gap-file=FILE          write sealed gaps to FILE\n"
 "  -t, --trace-file=FILE        write graph search stats to FILE\n"
 "  -v, --verbose                display verbose output\n"
 "      --help                   display this help and exit\n"
@@ -205,6 +206,9 @@ namespace opt {
 	/** Output file for graph search stats */
 	static string tracefilePath;
 
+	/** Output file for sealed gaps */
+	static string gapfilePath;
+
 	/** Mask bases not in flanks */
 	static int mask = 0;
 
@@ -237,7 +241,7 @@ struct Counters {
 	size_t skipped;
 };
 
-static const char shortopts[] = "S:L:b:B:d:ef:F:G:i:Ij:k:lm:M:no:P:q:r:s:t:v";
+static const char shortopts[] = "S:L:b:B:d:ef:F:G:g:i:Ij:k:lm:M:no:P:q:r:s:t:v";
 
 enum { OPT_HELP = 1, OPT_VERSION };
 
@@ -274,6 +278,7 @@ static const struct option longopts[] = {
 	{ "read-name",        required_argument, NULL, 'r' },
 	{ "search-mem",       required_argument, NULL, 's' },
 	{ "trace-file",       required_argument, NULL, 't' },
+	{ "gap-file",         required_argument, NULL, 'g' },
 	{ "verbose",          no_argument, NULL, 'v' },
 	{ "help",             no_argument, NULL, OPT_HELP },
 	{ "version",          no_argument, NULL, OPT_VERSION },
@@ -416,8 +421,8 @@ template <typename Graph>
 string merge(const Graph& g,
 	unsigned k,
 	const Gap& gap,
-	FastaRecord &read1,
-	FastaRecord &read2,
+	const FastaRecord &read1,
+	const FastaRecord &read2,
 	const ConnectPairsParams& params,
 	Counters& g_count,
 	ofstream& traceStream)
@@ -573,12 +578,12 @@ void kRun(const ConnectPairsParams& params,
 	map<FastaRecord, map<FastaRecord, Gap> > &flanks,
 	unsigned &gapsclosed,
 	ofstream &logStream,
-	ofstream &traceStream)
+	ofstream &traceStream,
+	ofstream &gapStream)
 {
 	map<FastaRecord, map<FastaRecord, Gap> >::iterator read1_it;
 	map<FastaRecord, Gap>::iterator read2_it;
 	unsigned uniqueGapsClosed = 0;
-	bool success;
 
 	Counters g_count;
 	g_count.noStartOrGoalKmer = 0;
@@ -598,31 +603,50 @@ void kRun(const ConnectPairsParams& params,
 
 	printLog(logStream, "Flanks inserted into k run = " + IntToString(flanks.size()) + "\n");
 
-	for (read1_it = flanks.begin(); read1_it != flanks.end();) {
-		success = false;
+	int counter = 0;
+	vector<map<FastaRecord, map<FastaRecord, Gap> >::iterator> flanks_closed;
+#pragma omp parallel private(read1_it, read2_it) firstprivate(counter)
+	for (read1_it = flanks.begin(); read1_it != flanks.end(); ++read1_it) {
 		FastaRecord read1 = read1_it->first;
-		for (read2_it = flanks[read1].begin(); read2_it != flanks[read1].end(); read2_it++) {
+		bool success = false;
+		for (read2_it = flanks[read1].begin(); read2_it != flanks[read1].end(); ++read2_it, ++counter) {
+#if _OPENMP
+			if (counter % omp_get_num_threads() != omp_get_thread_num())
+				continue;
+#endif
 			FastaRecord read2 = read2_it->first;
 
 			int startposition = read2_it->second.gapStart();
 			string tempSeq = merge(g, k, read2_it->second, read1, read2, params, g_count, traceStream);
 			if (!tempSeq.empty()) {
 				success = true;
+#pragma omp critical (allmerged)
 				allmerged[read1.id.substr(0,read1.id.length()-2)][startposition]
 					= ClosedGap(read2_it->second, tempSeq);
-//#pragma omp atomic
-				gapsclosed++;
-//#pragma omp atomic
-				uniqueGapsClosed++;
-				if (gapsclosed % 100 == 0)
+#pragma omp atomic
+				++uniqueGapsClosed;
+#pragma omp critical (gapsclosed)
+				if (++gapsclosed % 100 == 0)
 					printLog(logStream, IntToString(gapsclosed) + " gaps closed so far\n");
+
+				if (!opt::gapfilePath.empty())
+#pragma omp critical (gapStream)
+					gapStream << ">" << read1.id.substr(0,read1.id.length()-2)
+						  << "_" << read2_it->second.gapStart() << "-" << read2_it->second.gapEnd()
+						  << " LN:i:" << tempSeq.length() << '\n'
+						  << tempSeq << '\n';
 			}
 		}
 		if (success) {
-			flanks.erase(read1_it++);
+#pragma omp critical (flanks_closed)
+			flanks_closed.push_back(read1_it);
 		}
-		else
-			read1_it++;
+	}
+
+	for (vector<map<FastaRecord, map<FastaRecord, Gap> >::iterator>::iterator it = flanks_closed.begin();
+			it != flanks_closed.end(); ++it) {
+		if (flanks.count((*it)->first) > 0)
+		    flanks.erase(*it);
 	}
 
 	printLog(logStream, IntToString(uniqueGapsClosed) + " unique gaps closed for k" + IntToString(k) + "\n");
@@ -772,6 +796,8 @@ int main(int argc, char** argv)
 			opt::searchMem = SIToBytes(arg); break;
 		  case 't':
 			arg >> opt::tracefilePath; break;
+		  case 'g':
+		    arg >> opt::gapfilePath; break;
 		  case 'v':
 			opt::verbose++; break;
 		  case OPT_HELP:
@@ -875,6 +901,13 @@ int main(int argc, char** argv)
 		assert_good(traceStream, opt::tracefilePath);
 	}
 
+	ofstream gapStream;
+	if (!opt::gapfilePath.empty()) {
+		gapStream.open(opt::gapfilePath.c_str());
+		assert(gapStream.is_open());
+		assert_good(gapStream, opt::gapfilePath);
+	}
+
 	string logOutputPath(opt::outputPrefix);
 	logOutputPath.append("_log.txt");
 	ofstream logStream(logOutputPath.c_str());
@@ -1021,7 +1054,7 @@ int main(int argc, char** argv)
 		temp = "Starting K run with k = " + IntToString(opt::k) + "\n";
 		printLog(logStream, temp);
 
-		kRun(params, opt::k, g, allmerged, flanks, gapsclosed, logStream, traceStream);
+		kRun(params, opt::k, g, allmerged, flanks, gapsclosed, logStream, traceStream, gapStream);
 
 		temp = "k" + IntToString(opt::k) + " run complete\n"
 				+ "Total gaps closed so far = " + IntToString(gapsclosed) + "\n\n";
@@ -1074,5 +1107,10 @@ int main(int argc, char** argv)
 		traceStream.close();
 	}
 
+	if (!opt::gapfilePath.empty()) {
+		assert_good(gapStream, opt::gapfilePath);
+		gapStream.close();
+	}
+
 	return 0;
 }