Skip to main content
  • Home
  • Development
  • Documentation
  • Donate
  • Operational login
  • Browse the archive

swh logo
SoftwareHeritage
Software
Heritage
Archive
Features
  • Search

  • Downloads

  • Save code now

  • Add forge now

  • Help

swh:1:snp:76adc8cf123d92909f09e46db9ddfcbecf4cd1aa
  • Code
  • Branches (12)
  • Releases (0)
    • Branches
    • Releases
    • HEAD
    • refs/heads/azure-macOS
    • refs/heads/low_kmer_sample_check
    • refs/heads/master
    • refs/heads/refactor
    • refs/heads/spacedseeds
    • refs/tags/1.0.0
    • refs/tags/1.0.1
    • refs/tags/1.2.0
    • refs/tags/1.2.1
    • refs/tags/1.2.2
    • refs/tags/v1.1.0
    • refs/tags/v1.1.1
    No releases to show
  • 6673230
  • /
  • ntcard.cpp
Raw File Download
Permalinks

To reference or cite the objects present in the Software Heritage archive, permalinks based on SoftWare Hash IDentifiers (SWHIDs) must be used.
Select below a type of object currently browsed in order to display its associated SWHID and permalink.

  • content
  • directory
  • revision
  • snapshot
content badge Iframe embedding
swh:1:cnt:328abeacae8edd44cba5b39c2e91077b4ed20fc0
directory badge Iframe embedding
swh:1:dir:6673230adfb9c6ca6edd729c0ab8ba16b9e9de00
revision badge
swh:1:rev:9433d03872380c13bb131abf1e388e94862a4f15
snapshot badge
swh:1:snp:76adc8cf123d92909f09e46db9ddfcbecf4cd1aa
Citations

This interface enables to generate software citations, provided that the root directory of browsed objects contains a citation.cff or codemeta.json file.
Select below a type of object currently browsed in order to generate citations for them.

  • content
  • directory
  • revision
  • snapshot
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Tip revision: 9433d03872380c13bb131abf1e388e94862a4f15 authored by vlad0x00 on 05 January 2022, 21:31:33 UTC
Update formatting
Tip revision: 9433d03
ntcard.cpp
#include "Uncompress.h"
#include "vendor/ntHash/ntHashIterator.hpp"
#include "vendor/ntHash/stHashIterator.hpp"

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <getopt.h>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

#ifdef _OPENMP
#include <omp.h>
#endif

#define PROGRAM "ntCard"

static const char VERSION_MESSAGE[] =
    PROGRAM " 1.2.2 \n"
            "Written by Hamid Mohamadi.\n"
            "Copyright 2018 Hamid Mohamadi, Licensed under MIT License\n";

static const char USAGE_MESSAGE[] =
    "Usage: " PROGRAM " [OPTION]... FILE(S)...\n"
    "Estimates k-mer coverage histogram in FILE(S).\n\n"
    "Acceptable file formats: fastq, fasta, sam, bam and in compressed formats gz, bz, zip, xz.\n"
    "A list of files containing file names in each row can be passed with @ prefix.\n"

    "\n"
    " Options:\n"
    "\n"
    "  -t, --threads=N	use N parallel threads [1] (N>=2 should be used when input files are >=2)\n"
    "  -k, --kmer=N	the length of kmer \n"
    "  -g, --gap=N	the length of gap in the gap seed [0]. g mod 2 must equal k mod 2 unless g == "
    "0\n"
    "           	-g does not support multiple k currently.\n"
    "  -c, --cov=N	the maximum coverage of kmer in output [1000]\n"
    "  -p, --pref=STRING    the prefix for output file name(s)\n"
    "  -o, --output=STRING	the name for output file name (used when output should be a single "
    "file)\n"
    "      --help	display this help and exit\n"
    "      --version	output version information and exit\n"
    "\n"
    "Report bugs to https://github.com/bcgsc/ntCard/issues\n";

using namespace std;

namespace opt {
unsigned nThrd = 1;
unsigned kmLen = 64;
unsigned gap = 0;
size_t rBuck;
unsigned rBits = 27;
unsigned sBits = 11;
unsigned sMask;
unsigned covMax = 1000;
size_t nSamp = 2;
size_t nK = 0;
string prefix;
string output;
bool samH = true;
std::vector<std::vector<unsigned> > seedSet;
}

static const char shortopts[] = "t:s:r:k:c:l:p:f:o:g:";

enum
{
	OPT_HELP = 1,
	OPT_VERSION
};

static const struct option longopts[] = { { "threads", required_argument, NULL, 't' },
	                                      { "kmer", required_argument, NULL, 'k' },
	                                      { "gap", required_argument, NULL, 'g' },
	                                      { "cov", required_argument, NULL, 'c' },
	                                      { "rbit", required_argument, NULL, 'r' },
	                                      { "sbit", required_argument, NULL, 's' },
	                                      { "output", required_argument, NULL, 'o' },
	                                      { "pref", required_argument, NULL, 'p' },
	                                      { "help", no_argument, NULL, OPT_HELP },
	                                      { "version", no_argument, NULL, OPT_VERSION },
	                                      { NULL, 0, NULL, 0 } };

size_t
getInf(const char* inFile)
{
	std::ifstream in(inFile, std::ifstream::ate | std::ifstream::binary);
	return in.tellg();
}

bool
isNumber(const string& seq)
{
	string::const_iterator itr = seq.begin();
	while (itr != seq.end() && isdigit(*itr))
		++itr;
	return (itr == seq.end() && !seq.empty());
}

unsigned
getftype(std::ifstream& in, std::string& samSeq)
{
	std::string hseq;
	getline(in, hseq);
	if (hseq[0] == '>') {
		return 1;
	}
	if (hseq[0] == '@') {
		if ((hseq[1] == 'H' && hseq[2] == 'D') || (hseq[1] == 'S' && hseq[2] == 'Q') ||
		    (hseq[1] == 'R' && hseq[2] == 'G') || (hseq[1] == 'P' && hseq[2] == 'G') ||
		    (hseq[1] == 'C' && hseq[2] == 'O')) {
			return 2;
		} else
			return 0;
	}
	std::istringstream alnSec(hseq);
	std::string s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11;
	alnSec >> s1 >> s2 >> s3 >> s4 >> s5 >> s6 >> s7 >> s8 >> s9 >> s10 >> s11;
	if (isNumber(s2) && isNumber(s5)) {
		opt::samH = false;
		samSeq = hseq;
		return 2;
	}
	return 3;
}

inline void
ntComp(const uint64_t hVal, uint16_t* t_Counter)
{
	uint64_t indBit = opt::nSamp;
	if (hVal >> (63 - opt::sBits) == 1)
		indBit = 0;
	if (hVal >> (64 - opt::sBits) == opt::sMask)
		indBit = 1;
	if (indBit < opt::nSamp) {
		size_t shVal = hVal & (opt::rBuck - 1);
#pragma omp atomic
		++t_Counter[indBit * opt::rBuck + shVal];
	}
}

inline void
ntRead(const string& seq, const std::vector<unsigned>& kList, uint16_t* t_Counter, size_t totKmer[])
{
	for (unsigned k = 0; k < kList.size(); k++) {
		ntHashIterator itr(seq, 1, kList[k]);
		while (itr != itr.end()) {
			ntComp((*itr)[0], t_Counter + k * opt::nSamp * opt::rBuck);
			++itr;
			++totKmer[k];
		}
	}
}

inline void
stRead(const string& seq, const std::vector<unsigned>& kList, uint16_t* t_Counter, size_t totKmer[])
{
	for (unsigned k = 0; k < kList.size(); k++) {
		stHashIterator itr(seq, opt::seedSet, 1, 1, kList[k]);
		while (itr != itr.end()) {
			ntComp((*itr)[0], t_Counter + k * opt::nSamp * opt::rBuck);
			++itr;
			++totKmer[k];
		}
	}
}

void
getEfq(std::ifstream& in, const std::vector<unsigned>& kList, uint16_t* t_Counter, size_t totKmer[])
{
	bool good = true;
	for (string seq, hseq; good;) {
		good = static_cast<bool>(getline(in, seq));
		good = static_cast<bool>(getline(in, hseq));
		good = static_cast<bool>(getline(in, hseq));
		if (good && opt::gap == 0) {
			ntRead(seq, kList, t_Counter, totKmer);
		}
		if (good && opt::gap != 0) {
			stRead(seq, kList, t_Counter, totKmer);
		}
		good = static_cast<bool>(getline(in, hseq));
	}
}

void
getEfa(std::ifstream& in, const std::vector<unsigned>& kList, uint16_t* t_Counter, size_t totKmer[])
{
	bool good = true;
	for (string seq, hseq; good;) {
		string line;
		good = static_cast<bool>(getline(in, seq));
		while (good && seq[0] != '>') {
			line += seq;
			good = static_cast<bool>(getline(in, seq));
		}
		if (opt::gap == 0) {
			ntRead(line, kList, t_Counter, totKmer);
		} else {
			stRead(line, kList, t_Counter, totKmer);
		}
	}
}

void
getEsm(
    std::ifstream& in,
    const std::vector<unsigned>& kList,
    const std::string& samSeq,
    uint16_t* t_Counter,
    size_t totKmer[])
{
	std::string samLine, seq;
	std::string s1, s2, s3, s4, s5, s6, s7, s8, s9, s11;
	if (opt::samH) {
		while (getline(in, samLine))
			if (samLine[0] != '@')
				break;
	} else
		samLine = samSeq;
	do {
		std::istringstream iss(samLine);
		iss >> s1 >> s2 >> s3 >> s4 >> s5 >> s6 >> s7 >> s8 >> s9 >> seq >> s11;
		if (opt::gap == 0) {
			ntRead(seq, kList, t_Counter, totKmer);
		} else {
			stRead(seq, kList, t_Counter, totKmer);
		}
	} while (getline(in, samLine));
}

void
compEst(const uint16_t* t_Counter, double& F0Mean, double fMean[])
{
	unsigned p[opt::nSamp][65536];
	for (size_t i = 0; i < opt::nSamp; i++)
		for (size_t j = 0; j < 65536; j++)
			p[i][j] = 0;

	for (size_t i = 0; i < opt::nSamp; i++)
		for (size_t j = 0; j < opt::rBuck; j++)
			++p[i][t_Counter[i * opt::rBuck + j]];

	double pMean[65536];
	for (size_t i = 0; i < 65536; i++)
		pMean[i] = 0.0;
	for (size_t i = 0; i < 65536; i++) {
		for (size_t j = 0; j < opt::nSamp; j++)
			pMean[i] += p[j][i];
		pMean[i] /= 1.0 * opt::nSamp;
	}

	F0Mean = (ssize_t)(
	    (opt::rBits * log(2) - log(pMean[0])) * 1.0 * ((uint64_t)1 << (opt::sBits + opt::rBits)));
	for (size_t i = 0; i < 65536; i++)
		fMean[i] = 0;
	if (pMean[0] * (log(pMean[0]) - opt::rBits * log(2)) == 0) {
		return;
	}
	fMean[1] = -1.0 * pMean[1] / (pMean[0] * (log(pMean[0]) - opt::rBits * log(2)));
	for (size_t i = 2; i < 65536; i++) {
		double sum = 0.0;
		for (size_t j = 1; j < i; j++)
			sum += j * pMean[i - j] * fMean[j];
		fMean[i] = -1.0 * pMean[i] / (pMean[0] * (log(pMean[0]) - opt::rBits * log(2))) -
		           sum / (i * pMean[0]);
	}
	for (size_t i = 1; i < 65536; i++)
		fMean[i] = abs((ssize_t)(fMean[i] * F0Mean));
}

void
outDefault(const std::vector<unsigned>& kList, const size_t totalKmers[], const uint16_t* t_Counter)
{
	std::ofstream histFiles[opt::nK];
	for (unsigned k = 0; k < opt::nK; k++) {
		std::stringstream hstm;
		hstm << opt::prefix << "_k" << kList[k] << ".hist";
		histFiles[k].open((hstm.str()).c_str());
	}
#pragma omp parallel for schedule(dynamic)
	for (unsigned k = 0; k < kList.size(); k++) {
		double F0Mean = 0.0;
		double fMean[65536];
		compEst(t_Counter + k * opt::nSamp * opt::rBuck, F0Mean, fMean);
		histFiles[k] << "F1\t" << totalKmers[k] << "\n";
		histFiles[k] << "F0\t" << (uint64_t)F0Mean << "\n";
		for (size_t i = 1; i <= opt::covMax; i++)
			histFiles[k] << i << "\t" << (uint64_t)fMean[i] << "\n";
	}
	for (unsigned k = 0; k < opt::nK; k++)
		histFiles[k].close();
}

void
outCompact(const std::vector<unsigned>& kList, const size_t totalKmers[], const uint16_t* t_Counter)
{
	ofstream histFile(opt::output.c_str());
	histFile << "k\tf\tn\n";
	for (unsigned k = 0; k < kList.size(); k++) {
		double F0Mean = 0.0;
		double fMean[65536];
		compEst(t_Counter + k * opt::nSamp * opt::rBuck, F0Mean, fMean);
		std::cerr << "k=" << kList[k] << "\tF1\t" << totalKmers[k] << "\n";
		std::cerr << "k=" << kList[k] << "\tF0\t" << (uint64_t)F0Mean << "\n";
		for (size_t i = 1; i <= opt::covMax; i++)
			histFile << kList[k] << "\t" << i << "\t" << (uint64_t)fMean[i] << "\n";
	}
	histFile.close();
}

int
main(int argc, char** argv)
{

	double sTime = omp_get_wtime();

	vector<unsigned> kList;
	bool die = false;
	for (int c; (c = getopt_long(argc, argv, shortopts, longopts, NULL)) != -1;) {
		std::istringstream arg(optarg != NULL ? optarg : "");
		switch (c) {
		case '?':
			die = true;
			break;
		case 't':
			arg >> opt::nThrd;
			break;
		case 's':
			arg >> opt::sBits;
			break;
		case 'r':
			arg >> opt::rBits;
			break;
		case 'c':
			arg >> opt::covMax;
			if (opt::covMax > 65535)
				opt::covMax = 65535;
			break;
		case 'p':
			arg >> opt::prefix;
			break;
		case 'o':
			arg >> opt::output;
			break;
		case 'g':
			arg >> opt::gap;
			break;
		case 'k': {
			std::string token;
			while (getline(arg, token, ',')) {
				unsigned myK;
				std::stringstream ss(token);
				ss >> myK;
				kList.push_back(myK);
				++opt::nK;
			}
			break;
		}
		case OPT_HELP:
			std::cerr << USAGE_MESSAGE;
			exit(EXIT_SUCCESS);
		case OPT_VERSION:
			std::cerr << VERSION_MESSAGE;
			exit(EXIT_SUCCESS);
		}
		if (optarg != NULL && !arg.eof()) {
			std::cerr << PROGRAM ": invalid option: `-" << (char)c << optarg << "'\n";
			exit(EXIT_FAILURE);
		}
	}
	if (argc - optind < 1) {
		std::cerr << PROGRAM ": missing arguments\n";
		die = true;
	}

	if (opt::gap != 0 && (opt::gap % 2 != kList[0] % 2)) {
		std::cerr << PROGRAM "Gap size and kmer must have the same modulus\n";
		die = true;
	}

	if (opt::nK == 0) {
		std::cerr << PROGRAM ": missing argument -k ... \n";
		die = true;
	}

	if (opt::prefix.empty() && opt::output.empty()) {
		std::cerr << PROGRAM ": missing argument -p/-o ... \n";
		die = true;
	}

	if (opt::gap != 0 && kList.size() != 1) {
		std::cerr << PROGRAM ": -g does not support multiple k currently.\n";
		die = true;
	}

	if (die) {
		std::cerr << "Try `" << PROGRAM << " --help' for more information.\n";
		exit(EXIT_FAILURE);
	}

	if (opt::gap != 0) {
		std::string gap(opt::gap, '0');
		std::string nonGap((kList[0] - opt::gap) / 2, '1');
		std::vector<std::string> seedString;
		seedString.push_back(nonGap + gap + nonGap);
		opt::seedSet = stHashIterator::parseSeed(seedString);
	}

	vector<string> inFiles;
	for (int i = optind; i < argc; ++i) {
		string file(argv[i]);
		if (file[0] == '@') {
			string inName;
			ifstream inList(file.substr(1, file.length()).c_str());
			while (getline(inList, inName))
				inFiles.push_back(inName);
		} else
			inFiles.push_back(file);
	}

	uint64_t totalSize = 0;
	for (unsigned file_i = 0; file_i < inFiles.size(); ++file_i)
		totalSize += getInf(inFiles[file_i].c_str());
	if (totalSize < 50000000000)
		opt::sBits = 7;

	size_t totalKmers[kList.size()];
	for (unsigned k = 0; k < kList.size(); k++)
		totalKmers[k] = 0;

	opt::rBuck = ((size_t)1) << opt::rBits;
	opt::sMask = (((size_t)1) << (opt::sBits - 1)) - 1;
	uint16_t* t_Counter = new uint16_t[opt::nK * opt::nSamp * opt::rBuck]();

#ifdef _OPENMP
	omp_set_num_threads(opt::nThrd);
#endif

#pragma omp parallel for schedule(dynamic)
	for (unsigned file_i = 0; file_i < inFiles.size(); ++file_i) {
		size_t totKmer[kList.size()];
		for (unsigned k = 0; k < kList.size(); k++)
			totKmer[k] = 0;
		std::ifstream in(inFiles[file_i].c_str());
		std::string samSeq;
		unsigned ftype = getftype(in, samSeq);
		if (ftype == 0)
			getEfq(in, kList, t_Counter, totKmer);
		else if (ftype == 1)
			getEfa(in, kList, t_Counter, totKmer);
		else if (ftype == 2)
			getEsm(in, kList, samSeq, t_Counter, totKmer);
		else {
			std::cerr << "Error in reading file: " << inFiles[file_i] << std::endl;
			exit(EXIT_FAILURE);
		}
		in.close();
		for (unsigned k = 0; k < kList.size(); k++)
#pragma omp atomic
			totalKmers[k] += totKmer[k];
	}

	if (opt::output.empty())
		outDefault(kList, totalKmers, t_Counter);
	else
		outCompact(kList, totalKmers, t_Counter);

	delete[] t_Counter;

	std::cerr << "Runtime(sec): " << setprecision(4) << fixed << omp_get_wtime() - sTime << "\n";
	return 0;
}

Software Heritage — Copyright (C) 2015–2025, The Software Heritage developers. License: GNU AGPLv3+.
The source code of Software Heritage itself is available on our development forge.
The source code files archived by Software Heritage are available under their own copyright and licenses.
Terms of use: Archive access, API— Contact— JavaScript license information— Web API

back to top