roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Add copy buttons to all \u003cpre\u003e\u003ccode\u003e blocks\n(function() {\n function addCopyButtons() {\n document.querySelectorAll('pre code').forEach(function(codeBlock) {\n if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;\n codeBlock.parentElement.setAttribute('data-copy-added', 'true');\n \n var btn = document.createElement('button');\n btn.textContent = 'Copy';\n btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';\n btn.onmouseover = function() { this.style.opacity = '1'; };\n btn.onmouseout = function() { this.style.opacity = '0.7'; };\n btn.onclick = function() {\n navigator.clipboard.writeText(codeBlock.textContent).then(function() {\n btn.textContent = 'Copied!';\n setTimeout(function() { btn.textContent = 'Copy'; }, 1500);\n });\n };\n codeBlock.parentElement.style.position = 'relative';\n codeBlock.parentElement.appendChild(btn);\n });\n }\n \n addCopyButtons();\n \n // Re-run on dynamic content\n var observer = new MutationObserver(addCopyButtons);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Add Copy Buttons to Code Blocks"); } } catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); } })(); (function(){ try { var __m = "github.com"; var __re = new RegExp('^' + "github\\.com" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Force GitHub README to respect dark mode\n(function() {\n var style = document.createElement('style');\n style.textContent = '\n .markdown-body {\n color-scheme: dark light;\n }\n .markdown-body pre { background: #161b22 !important; }\n .markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }\n .markdown-body table th, .markdown-body table td { border-color: #30363d !important; }\n .markdown-body img { background: #0d1117; }\n .markdown-body blockquote { border-left-color: #8b949e; }\n .markdown-body hr { border-color: #30363d; }\n ';\n document.head.appendChild(style);\n})();", "GitHub Dark Mode README Fix"); } } catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Highlight search terms from Google/DuckDuckGo/Bing referrer\n(function() {\n var ref = document.referrer;\n var terms = [];\n \n if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {\n var url = new URL(ref);\n var q = url.searchParams.get('q') || url.searchParams.get('p');\n if (q) {\n terms = q.split(/\\s+/).filter(function(t) { return t.length \u003e 2; });\n }\n }\n \n if (terms.length === 0) return;\n \n var style = document.createElement('style');\n style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';\n document.head.appendChild(style);\n \n function highlight(node) {\n if (node.nodeType === 3) { // text node\n var text = node.textContent;\n var found = false;\n terms.forEach(function(term) {\n var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\\]\\\\]/g, '\\\\') + ')', 'gi');\n if (regex.test(text)) {\n found = true;\n var frag = document.createDocumentFragment();\n var parts = text.split(regex);\n parts.forEach(function(part, i) {\n if (i % 2 === 0) {\n frag.appendChild(document.createTextNode(part));\n } else {\n var span = document.createElement('span');\n span.className = 'userscript-highlight';\n span.textContent = part;\n frag.appendChild(span);\n }\n });\n node.parentNode.replaceChild(frag, node);\n }\n });\n } else if (node.nodeType === 1 && node.childNodes) { // element\n var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];\n if (!skipTags.includes(node.tagName)) {\n Array.from(node.childNodes).forEach(highlight);\n }\n }\n }\n \n highlight(document.body);\n \n // Re-highlight on dynamic content\n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1 || node.nodeType === 3) highlight(node);\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Highlight Search Terms"); } } catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Strip utm_, fbclid, gclid, etc. from all links on page\n(function() {\n var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',\n 'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',\n 'ref', 'ref_src', 'source', 'medium', 'campaign'];\n \n function cleanUrl(url) {\n try {\n var u = new URL(url, window.location.origin);\n var changed = false;\n trackingParams.forEach(function(p) {\n if (u.searchParams.has(p)) {\n u.searchParams.delete(p);\n changed = true;\n }\n });\n return changed ? u.toString() : url;\n } catch (e) {\n return url;\n }\n }\n \n function cleanLinks() {\n document.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n \n cleanLinks();\n \n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1) {\n if (node.tagName === 'A') cleanLinks();\n node.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Remove Tracking Parameters from Links"); } } catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); } })(); (function(){ try { var __m = "youtube.com"; var __re = new RegExp('^' + "youtube\\.com" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Auto-enable theater mode on YouTube\n(function() {\n function tryTheater() {\n var btn = document.querySelector('button[aria-label=\"Theater mode\"], ytd-player #player button[title=\"Theater mode\"]');\n if (btn && !btn.classList.contains('activated')) {\n btn.click();\n }\n }\n \n // Try immediately\n tryTheater();\n \n // Try after navigation (SPA)\n var lastUrl = location.href;\n setInterval(function() {\n if (location.href !== lastUrl) {\n lastUrl = location.href;\n setTimeout(tryTheater, 500);\n }\n }, 1000);\n \n // Also try on player load\n var observer = new MutationObserver(tryTheater);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "YouTube Theater Mode Default"); } } catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Remove or un-stick sticky/fixed headers that block content\n(function() {\n function unstick() {\n document.querySelectorAll('header, nav, [role=\"banner\"], .header, .navbar, .sticky, .fixed-top, [style*=\"position: fixed\"], [style*=\"position:sticky\"]').forEach(function(el) {\n if (el.style.position === 'fixed' || el.style.position === 'sticky' || \n getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {\n el.style.position = 'static';\n el.style.top = 'auto';\n el.style.zIndex = 'auto';\n }\n });\n }\n \n unstick();\n \n var observer = new MutationObserver(unstick);\n observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });\n})();", "Kill Sticky Headers"); } } catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Universal Dark Mode - works on any site\n(function() {\n var enabled = true;\n \n function applyDarkMode() {\n if (!enabled) return;\n \n // Create style element if it doesn't exist\n var style = document.getElementById('universal-dark-mode-style');\n if (!style) {\n style = document.createElement('style');\n style.id = 'universal-dark-mode-style';\n document.head.appendChild(style);\n }\n \n // Dark mode CSS - inverts colors but preserves images/video\n style.textContent = '\n /* Invert everything except media */\n html {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #1a1a2e !important;\n }\n \n /* Restore images, videos, iframes, canvas */\n img, video, iframe, canvas, svg, picture, [style*=\"background-image\"] {\n filter: invert(1) hue-rotate(180deg) !important;\n }\n \n /* Preserve specific elements that should not be inverted */\n .no-dark-mode, .no-dark-mode *,\n [data-theme=\"light\"], [data-theme=\"light\"],\n .ace_editor, .ace_editor *,\n .CodeMirror, .CodeMirror *,\n .monaco-editor, .monaco-editor *,\n .markdown-body pre, .markdown-body pre *,\n .highlight, .highlight *,\n pre code, pre code * {\n filter: none !important;\n }\n \n /* Fix common UI elements */\n .modal, .popup, .dropdown-menu, .tooltip, .popover {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #2d2d44 !important;\n border-color: #444 !important;\n }\n \n /* Scrollbars */\n ::-webkit-scrollbar { background: #1a1a2e !important; }\n ::-webkit-scrollbar-thumb { background: #444 !important; }\n ::-webkit-scrollbar-thumb:hover { background: #555 !important; }\n \n /* Selection */\n ::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ';\n }\n \n function removeDarkMode() {\n var style = document.getElementById('universal-dark-mode-style');\n if (style) style.remove();\n }\n \n // Toggle with Alt+Shift+D\n document.addEventListener('keydown', function(e) {\n if (e.altKey && e.shiftKey && e.key === 'D') {\n e.preventDefault();\n enabled = !enabled;\n if (enabled) {\n applyDarkMode();\n console.log('[Universal Dark Mode] Enabled');\n } else {\n removeDarkMode();\n console.log('[Universal Dark Mode] Disabled');\n }\n }\n });\n \n // Apply on load\n applyDarkMode();\n \n // Re-apply on dynamic content\n var observer = new MutationObserver(function(mutations) {\n if (enabled && !document.getElementById('universal-dark-mode-style')) {\n applyDarkMode();\n }\n });\n observer.observe(document.head, { childList: true });\n \n console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');\n})();", "Universal Dark Mode"); } } catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); } })(); })();
Skip to content

roi accepts xsym input#235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Jun 23, 2024
Merged

roi accepts xsym input #235

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 164 additions & 42 deletions src/analysis/roimethstat.cpp
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
/* roimethstat: average methylation in each of a set of regions
*
* Copyright (C) 2014-2023 Andrew D. Smith
* Copyright (C) 2014-2024 Andrew D. Smith
*
* Authors: Andrew D. Smith and Masaru Nakajima
*
Expand All@@ -24,16 +24,19 @@
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <filesystem>
#include <charconv>

#include "GenomicRegion.hpp"
#include "LevelsCounter.hpp"
#include "MSite.hpp"
#include "OptionParser.hpp"
#include "bsutils.hpp"
#include "smithlab_utils.hpp"
#include "xcounts_utils.hpp"

using std::cerr;
using std::cout;
Expand All@@ -46,12 +49,144 @@ using std::runtime_error;
using std::string;
using std::to_string;
using std::unordered_map;
using std::unordered_set;
using std::vector;
using std::from_chars;
using std::size;
using std::cend;
using std::ostream;
using std::size;

using bamxx::bgzf_file;

namespace fs = std::filesystem;


static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}


struct genomic_interval {
string chrom{};
uint64_t start_pos{};
uint64_t end_pos{};
};


static void
update(LevelsCounter &lc, const xcounts_entry &xse) {
const uint64_t n_reads = xse.n_meth + xse.n_unmeth;
if (n_reads > 0) {
++lc.sites_covered;
lc.max_depth = std::max(lc.max_depth, n_reads);
lc.total_c += xse.n_meth;
lc.total_t += xse.n_unmeth;
const auto meth = static_cast<double>(xse.n_unmeth) / n_reads;
lc.total_meth += meth;
double lower = 0.0, upper = 0.0;
wilson_ci_for_binomial(lc.alpha, n_reads, meth, lower, upper);
lc.called_meth += (lower > 0.5);
lc.called_unmeth += (upper < 0.5);
}
++lc.total_sites;
}

static void
process_chrom(const bool report_more_info, const char level_code,
const vector<GenomicRegion> &intervals,
const vector<xcounts_entry> &sites,
ostream &out) {

uint64_t j = 0;
for (auto i = 0ul; i < intervals.size(); ++i) {
while (j < size(sites) && sites[j].pos < intervals[i].get_start()) ++j;

LevelsCounter lc;
while (j < size(sites) && sites[j].pos < intervals[i].get_end())
update(lc, sites[j++]);

GenomicRegion r(intervals[i]);
r.set_score(level_code == 'w' ? lc.mean_meth_weighted()
: (level_code == 'u' ? lc.mean_meth()
: lc.fractional_meth()));
r.set_name("X_" +
std::to_string((level_code == 'w'
? lc.coverage()
: (level_code == 'u' ? lc.sites_covered
: lc.total_called()))));
out << r;
if (report_more_info) out << '\t' << format_levels_counter(lc);
out << '\n';
}
}

static void
process_chrom(const bool report_more_info,
const vector<GenomicRegion> &intervals,
ostream &out) {
LevelsCounter lc;
const string lc_formatted = format_levels_counter(lc);
for (const auto &r: intervals) {
out << r;
if (report_more_info) out << '\t' << lc_formatted;
out << '\n';
}
}


static void
process_from_xcounts(const uint32_t n_threads,
const bool report_more_info,
const char level_code, const string &xsym_file,
const vector<GenomicRegion> &intervals,
ostream &out) {

const auto sites_by_chrom = read_xcounts_by_chrom(n_threads, xsym_file);
// const auto intervals = get_GenomicRegions(intervals_file);

vector<vector<GenomicRegion>> intervals_by_chrom;
string prev_chrom;
for (auto i = 0u; i < size(intervals); ++i) {
if (intervals[i].get_chrom() != prev_chrom) {
intervals_by_chrom.push_back(vector<GenomicRegion>());
prev_chrom = intervals[i].get_chrom();
}
intervals_by_chrom.back().push_back(intervals[i]);
}

for (const auto &intervals : intervals_by_chrom) {
const auto chrom_name = intervals.front().get_chrom();
const auto sites = sites_by_chrom.find(chrom_name);
if (sites != cend(sites_by_chrom))
process_chrom(report_more_info, level_code,
intervals, sites->second, out);
else
process_chrom(report_more_info, intervals, out);
}
}



bool
cmp_within_chrom(const GenomicRegion &r1, const GenomicRegion &r2) {
return (r1.get_start() < r2.get_start() ||
Expand DownExpand Up@@ -142,28 +277,6 @@ region_bounds(const unordered_map<string, uint32_t> &chrom_order,
return {lower_bound(first, last, a, cmp), lower_bound(first, last, b, cmp)};
}

static string
format_levels_counter(const LevelsCounter &lc) {
// ...
// (7) weighted mean methylation
// (8) unweighted mean methylation
// (9) fractional methylation
// (10) number of sites in the region
// (11) number of sites covered at least once
// (12) number of observations in reads indicating methylation
// (13) total number of observations from reads in the region
std::ostringstream oss;
// clang-format off
oss << lc.mean_meth_weighted() << '\t'
<< lc.mean_meth() << '\t'
<< lc.fractional_meth() << '\t'
<< lc.total_sites << '\t'
<< lc.sites_covered << '\t'
<< lc.total_c << '\t'
<< (lc.total_c + lc.total_t);
// clang-format on
return oss.str();
}

static bool
is_sorted_within_chrom(const vector<MSite> &sites) {
Expand DownExpand Up@@ -194,10 +307,10 @@ read_sites(const string &filename) {
}

static void
process_preloaded(const bool VERBOSE, const bool report_more_information,
process_preloaded(const bool VERBOSE, const bool report_more_info,
const char level_code, const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {

const auto sites = read_sites(sites_file);
if (sites.empty()) throw runtime_error("failed to read sites: " + sites_file);
Expand All@@ -219,7 +332,7 @@ process_preloaded(const bool VERBOSE, const bool report_more_information,
GenomicRegion r_scored{r};
r_scored.set_score(score);
out << r_scored;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -266,10 +379,10 @@ calc_site_stats(ifstream &sites_in, const GenomicRegion &region,
}

static void
process_on_disk(const bool report_more_information, const char level_code,
process_on_disk(const bool report_more_info, const char level_code,
const string &sites_file,
const unordered_map<string, uint32_t> &chrom_order,
const vector<GenomicRegion> &regions, std::ostream &out) {
const vector<GenomicRegion> &regions, ostream &out) {
ifstream in(sites_file);
if (!in) throw runtime_error("failed to open file: " + sites_file);

Expand All@@ -282,7 +395,7 @@ process_on_disk(const bool report_more_information, const char level_code,
GenomicRegion r{region};
r.set_score(score);
out << r;
if (report_more_information)
if (report_more_info)
out << '\t' << format_levels_counter(lc);
out << '\n';
}
Expand DownExpand Up@@ -328,8 +441,9 @@ Columns (beyond the first 6) in the BED format output:
bool VERBOSE = false;
bool print_numeric_only = false;
bool preload = false;
bool report_more_information = false;
bool report_more_info = false;
bool sort_data_if_needed = false;
uint32_t n_threads = 1;

string level_code = "w";

Expand All@@ -351,7 +465,9 @@ Columns (beyond the first 6) in the BED format output:
"in bed format output (w, u or f)",
false, level_code);
opt_parse.add_opt("more-levels", 'M', "report more methylation information",
false, report_more_information);
false, report_more_info);
opt_parse.add_opt("threads", 't', "threads to use (if input compressed)",
false, n_threads);
opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
vector<string> leftover_args;
opt_parse.parse(argc, argv, leftover_args);
Expand DownExpand Up@@ -380,14 +496,17 @@ Columns (beyond the first 6) in the BED format output:
const string sites_file = leftover_args.back();
/****************** END COMMAND LINE OPTIONS *****************/

if (!is_msite_file(sites_file))
throw runtime_error("dnmtools counts format required: " + sites_file);
const bool is_xcounts = get_is_xcounts_file(sites_file);
if (!is_msite_file(sites_file) && !is_xcounts)
throw runtime_error("dnmtools counts or xcounts format required: " +
sites_file);

// make a map that specifies their order; otherwise we can't
// ensure regions are sorted in the same way
unordered_map<string, uint32_t> chrom_order;
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());
if (!is_xcounts)
for (auto &i : get_chroms(sites_file))
chrom_order.emplace(i, chrom_order.size());

if (VERBOSE) cerr << "loading regions" << endl;

Expand DownExpand Up@@ -427,14 +546,17 @@ Columns (beyond the first 6) in the BED format output:
of.open(outfile);
if (!of) throw runtime_error("failed to open outfile: " + outfile);
}
std::ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (preload)
process_preloaded(VERBOSE, report_more_information, level_code[0],
sites_file, chrom_order, regions, out);
ostream out(outfile.empty() ? cout.rdbuf() : of.rdbuf());

if (is_xcounts)
process_from_xcounts(n_threads, report_more_info, level_code[0],
sites_file, regions, out);
else if (preload)
process_preloaded(VERBOSE, report_more_info, level_code[0], sites_file,
chrom_order, regions, out);
else
process_on_disk(report_more_information, level_code[0], sites_file,
chrom_order, regions, out);
process_on_disk(report_more_info, level_code[0], sites_file, chrom_order,
regions, out);
}
catch (const std::exception &e) {
cerr << e.what() << endl;
Expand Down