OSCR

Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches
  1. [1] § Methods › Alignment, QC, and CNV calling ↔ src/coral.h, lines 500–584 · score 0.60 · CNV calling, mapping quality, ploidy, Somatic, depth, genotyped
  2. [2] § Methods › Alignment, QC, and CNV calling ↔ workflow/scripts/haplotagging_scripts/haplotagTable.R, lines 16–63 · score 0.57 · findOverlaps, GRanges, exported, BED, Overlapping, filtered
  3. [3] § Methods › Demultiplexing and alignment ↔ transloc_pipeline_v3/bin/TranslocPipeline.pl, lines 1078–1147 · score 0.52 · TranslocWrapper, Bowtie2, pipeline, adapters, prey, bait

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

C/C++ header · 933 lines · 35 KB · BSD-3-Clause · 1 match

  1. #ifndef CORAL_H
  2. #define CORAL_H
  4. #include <limits>
  5. #include <iomanip>
  7. #include <boost/icl/split_interval_map.hpp>
  8. #include <boost/dynamic_bitset.hpp>
  9. #include <boost/unordered_map.hpp>
  10. #include <boost/date_time/posix_time/posix_time.hpp>
  11. #include <boost/date_time/gregorian/gregorian.hpp>
  12. #include <boost/math/special_functions/round.hpp>
  14. #include <htslib/sam.h>
  15. #include <htslib/faidx.h>
  17. #include "bed.h"
  18. #include "scan.h"
  19. #include "gcbias.h"
  20. #include "cnv.h"
  21. #include "version.h"
  23. namespace torali
  24. {
  26. struct CountDNAConfig {
  27. bool basecov;
  28. bool somatic;
  29. bool adaptive;
  30. bool hasStatsFile;
  31. bool hasScanFile;
  32. bool noScanWindowSelection;
  33. bool regionalGc;
  34. bool hasSegFile;
  35. bool hasGenoFile;
  36. bool hasExcludeFile;
  37. uint32_t nchr;
  38. uint32_t minClip;
  39. uint32_t minRefSep;
  40. uint32_t minBpSupport;
  41. float cnMergeTol;
  42. float penalty;
  43. uint32_t meanisize;
  44. uint32_t window_size;
  45. uint32_t window_offset;
  46. uint32_t scanWindow;
  47. uint32_t minChrLen;
  48. uint32_t minCnvSize;
  49. uint32_t targetReads;
  50. uint32_t minSegWin;
  51. double targetExpCov;
  52. uint16_t minQual;
  53. uint16_t mapqUniq;
  54. uint16_t mad;
  55. float ploidy;
  56. float ctrlPloidy;
  57. float expectedCN;
  58. float purity;
  59. float exclgc;
  60. float uniqueToTotalCovRatio;
  61. float fracWindow;
  62. float fragmentUnique;
  63. int32_t cnvMinQual;
  64. float cnvDelConfirm;
  65. float cnvDelRatio;
  66. std::string sampleName;
  67. boost::filesystem::path segfile;
  68. boost::filesystem::path genofile;
  69. boost::filesystem::path outfile;
  70. boost::filesystem::path covfile;
  71. boost::filesystem::path genome;
  72. boost::filesystem::path statsFile;
  73. boost::filesystem::path bamFile;
  74. boost::filesystem::path scanFile;
  75. boost::filesystem::path exclude;
  76. std::set<int32_t> refIdx;
  77. };
  79. struct CovWin {
  80. uint32_t start;
  81. uint32_t end;
  82. uint32_t winlen;
  83. double covsum;
  84. double expcov;
  85. double ucov;
  86. double tcov;
  87. double aall;
  88. double eall;
  89. bool valid;
  91. CovWin(uint32_t const s, uint32_t const e, uint32_t const w, double const cs, double const ec, double const uc, double const tc, double const aa, double const ea, bool vld) : start(s), end(e), winlen(w), covsum(cs), expcov(ec), ucov(uc), tcov(tc), aall(aa), eall(ea), valid(vld) {}
  92. };
  94. struct CountDNAConfigLib {
  95. uint16_t madCutoff;
  96. uint16_t madNormalCutoff;
  97. boost::filesystem::path genome;
  98. std::vector<boost::filesystem::path> files;
  99. };
  101. template<typename TConfig>
  102. inline int32_t
  103. bamCount(TConfig const& c, LibraryInfo const& li, std::vector<GcBias> const& gcbias, std::pair<uint32_t, uint32_t> const& gcbound, std::vector<double> const& regcorr, uint32_t const regWin) {
  104. // Load bam file
  105. samFile* samfile = sam_open(c.bamFile.string().c_str(), "r");
  106. hts_set_fai_filename(samfile, c.genome.string().c_str());
  107. hts_idx_t* idx = sam_index_load(samfile, c.bamFile.string().c_str());
  108. bam_hdr_t* hdr = sam_hdr_read(samfile);
  110. // Parse BAM file
  111. boost::posix_time::ptime now = boost::posix_time::second_clock::local_time();
  112. std::cerr << '[' << boost::posix_time::to_simple_string(now) << "] " << "Count fragments" << std::endl;
  114. // Open output files
  115. boost::iostreams::filtering_ostream dataOut;
  116. if (!c.covfile.empty()) {
  117. dataOut.push(boost::iostreams::gzip_compressor());
  118. dataOut.push(boost::iostreams::file_sink(c.covfile.string(), std::ios_base::out | std::ios_base::binary));
  119. dataOut << "chr\tstart\tend\t" << c.sampleName << "_uniqfrac\t" << c.sampleName << "_logR\t" << c.sampleName << "_CN" << std::endl;
  120. }
  122. // CNVs
  123. std::vector<CNV> cnvs;
  124. if (c.hasGenoFile) parseVcfCNV(c, hdr, cnvs);
  126. // Iterate chromosomes
  127. faidx_t* faiRef = fai_load(c.genome.string().c_str());
  128. for(int32_t refIndex=0; refIndex < (int32_t) hdr->n_targets; ++refIndex) {
  129. if ((!c.hasGenoFile) && (chrNoData(c, refIndex, idx))) continue;
  131. // Haploid chromosome?
  132. float chrCtrlPloidy = c.ctrlPloidy;
  133. float chrPloidy = c.ploidy;
  134. if (c.refIdx.find(refIndex) != c.refIdx.end()) {
  135. chrCtrlPloidy = c.ctrlPloidy - 1;
  136. chrPloidy = c.ploidy - 1;
  137. }
  139. // Reference sequence
  140. std::string tname(hdr->target_name[refIndex]);
  141. int32_t seqlen = faidx_seq_len(faiRef, tname.c_str());
  142. if (seqlen == - 1) continue;
  143. else seqlen = -1;
  144. char* ref = faidx_fetch_seq(faiRef, tname.c_str(), 0, faidx_seq_len(faiRef, tname.c_str()), &seqlen);
  145. if (ref == NULL) continue;
  147. // Get GC content
  148. std::vector<uint16_t> uniqContent(hdr->target_len[refIndex], 0);
  149. std::vector<uint16_t> gcContent(hdr->target_len[refIndex], 0);
  150. {
  151. // GC map
  152. typedef boost::dynamic_bitset<> TBitSet;
  153. TBitSet gcref(hdr->target_len[refIndex], false);
  154. for(uint32_t i = 0; i < hdr->target_len[refIndex]; ++i) {
  155. if ((ref[i] == 'c') || (ref[i] == 'C') || (ref[i] == 'g') || (ref[i] == 'G')) gcref[i] = 1;
  156. }
  158. // Sum across fragment
  159. int32_t halfwin = (int32_t) (c.meanisize / 2);
  160. int32_t gcsum = 0;
  161. for(int32_t pos = halfwin; pos < (int32_t) hdr->target_len[refIndex] - halfwin; ++pos) {
  162. if (pos == halfwin) {
  163. for(int32_t i = pos - halfwin; i<=pos+halfwin; ++i) gcsum += gcref[i];
  164. } else {
  165. gcsum -= gcref[pos - halfwin - 1];
  166. gcsum += gcref[pos + halfwin];
  167. }
  168. gcContent[pos] = gcsum;
  169. }
  170. }
  172. // Broad window for replication wave
  173. std::vector<float> tileFac;
  174. if ((!regcorr.empty()) && (regWin > 0)) {
  175. uint32_t ntile = hdr->target_len[refIndex] / regWin + 1;
  176. tileFac.resize(ntile, 1.0);
  177. for(uint32_t t = 0; t < ntile; ++t) {
  178. uint32_t s = t * regWin;
  179. uint32_t e = std::min((uint32_t) hdr->target_len[refIndex], s + regWin);
  180. double gcnum = 0;
  181. uint32_t winlen = 0;
  182. for(uint32_t pos = s; pos < e; ++pos) {
  183. if ((gcContent[pos] > gcbound.first) && (gcContent[pos] < gcbound.second)) { gcnum += gcContent[pos]; ++winlen; }
  184. }
  185. if (winlen > 0) tileFac[t] = (float) regCorrFactor(regcorr, (gcnum / (double) winlen) / (double) c.meanisize);
  186. }
  187. }
  189. // Coverage track
  190. typedef uint16_t TCount;
  191. uint32_t maxCoverage = std::numeric_limits<TCount>::max();
  192. typedef std::vector<TCount> TCoverage;
  193. TCoverage cov(hdr->target_len[refIndex], 0);
  194. TCoverage covUniq(hdr->target_len[refIndex], 0);
  195. TCoverage covAll(hdr->target_len[refIndex], 0);
  196. TCoverage covTot;
  197. if (!c.basecov) covTot.resize(hdr->target_len[refIndex], 0);
  198. TCoverage& covMap = (!c.basecov) ? covTot : cov;
  200. // Split-read breakpoints
  201. std::vector<int32_t> clips;
  202. {
  203. // Mate map
  204. typedef boost::unordered_map<std::size_t, bool> TMateMap;
  205. TMateMap mateMap;
  207. // Count reads
  208. hts_itr_t* iter = sam_itr_queryi(idx, refIndex, 0, hdr->target_len[refIndex]);
  209. bam1_t* rec = bam_init1();
  210. int32_t lastAlignedPos = 0;
  211. std::set<std::size_t> lastAlignedPosReads;
  212. while (sam_itr_next(samfile, iter, rec) >= 0) {
  213. if (rec->core.flag & (BAM_FQCFAIL | BAM_FDUP | BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) continue;
  214. if ((rec->core.flag & BAM_FPAIRED) && ((rec->core.flag & BAM_FMUNMAP) || (rec->core.tid != rec->core.mtid))) continue;
  216. // Base coverage
  217. addBaseCoverage3(rec, covAll, covMap, covUniq, c.minQual, c.mapqUniq, hdr->target_len[refIndex], maxCoverage);
  219. // Filter by quality
  220. if (rec->core.qual < c.minQual) continue;
  222. // Collect split-read breakpoints
  223. if (rec->core.qual >= c.mapqUniq) addSplitReadBreakpoints(rec, c.minClip, c.minRefSep, hdr->target_len[refIndex], clips);
  225. // Base-level counting done
  226. if (c.basecov) continue;
  228. // Fragment coverage
  229. int32_t midPoint = rec->core.pos + halfAlignmentLength(rec);
  230. if (rec->core.flag & BAM_FPAIRED) {
  231. std::size_t seed = hash_sr(rec);
  233. // Clean-up the read store for identical alignment positions
  234. if (rec->core.pos > lastAlignedPos) {
  235. lastAlignedPosReads.clear();
  236. lastAlignedPos = rec->core.pos;
  237. }
  239. if (_firstPairObs(rec, lastAlignedPosReads)) {
  240. // First read
  241. lastAlignedPosReads.insert(seed);
  242. std::size_t hv = hash_pair(rec);
  243. mateMap[hv] = true;
  244. continue;
  245. } else {
  246. // Second read
  247. std::size_t hv = hash_pair_mate(rec);
  248. auto itMM = mateMap.find(hv);
  249. if ((itMM == mateMap.end()) || (!itMM->second)) continue; // Mate discarded
  250. mateMap.erase(itMM);
  251. }
  253. // update midpoint
  254. int32_t isize = (rec->core.pos + alignmentLength(rec)) - rec->core.mpos;
  255. if ((li.minNormalISize < isize) && (isize < li.maxNormalISize)) midPoint = rec->core.mpos + (int32_t) (isize/2);
  256. }
  258. // Count fragment
  259. if ((midPoint >= 0) && (midPoint < (int32_t) hdr->target_len[refIndex]) && (cov[midPoint] < maxCoverage - 1)) ++cov[midPoint];
  260. }
  261. // Clean-up
  262. bam_destroy1(rec);
  263. hts_itr_destroy(iter);
  264. }
  266. // Callable positions: unique if mostly high-MAPQ reads
  267. for(uint32_t pos = 0; pos < hdr->target_len[refIndex]; ++pos) {
  268. bool u;
  269. if (covMap[pos] == 0) u = ((ref[pos] != 'N') && (ref[pos] != 'n'));
  270. else u = (2 * (uint32_t) covUniq[pos] >= (uint32_t) covMap[pos]);
  271. uniqContent[pos] = (u ? (uint16_t) c.meanisize : 0);
  272. }
  273. // hom-del or unmappable?
  274. uint32_t maxHomDel = 1000000;
  275. uint32_t rstart = 0;
  276. while (rstart < hdr->target_len[refIndex]) {
  277. if (covMap[rstart] == 0) {
  278. uint32_t rend = rstart;
  279. while ((rend < hdr->target_len[refIndex]) && (covMap[rend] == 0)) ++rend;
  280. bool leftOK = (rstart > 0) && (uniqContent[rstart - 1] > 0);
  281. bool rightOK = (rend < hdr->target_len[refIndex]) && (uniqContent[rend] > 0);
  282. if ((!leftOK) || (!rightOK) || (rend - rstart > maxHomDel)) {
  283. for(uint32_t k = rstart; k < rend; ++k) uniqContent[k] = 0;
  284. }
  285. rstart = rend;
  286. } else ++rstart;
  287. }
  288. // Genome-wide read-depth windows
  289. DepthTrack const dt = uniqueTrack();
  290. std::vector<CovWin> wins;
  291. if (c.adaptive) {
  292. double covsum = 0;
  293. double expraw = 0;
  294. double expcor = 0;
  295. double ucov = 0;
  296. double tcov = 0;
  297. double aall = 0;
  298. double eall = 0;
  299. uint32_t winlen = 0;
  300. uint32_t start = 0;
  301. for(uint32_t pos = 0; pos < hdr->target_len[refIndex]; ++pos) {
  302. ucov += covUniq[pos];
  303. tcov += covMap[pos];
  304. if ((gcContent[pos] > gcbound.first) && (gcContent[pos] < gcbound.second)) {
  305. aall += covAll[pos];
  306. eall += gcbias[gcContent[pos]].coverageTotal;
  307. }
  308. if (_posUsed(dt, c, gcContent, uniqContent, gcbound, pos)) {
  309. double e1 = _expCov(dt, gcbias, gcContent[pos]);
  310. covsum += cov[pos];
  311. expraw += e1;
  312. expcor += e1 * (tileFac.empty() ? 1.0 : (double) tileFac[pos / regWin]);
  313. ++winlen;
  314. if (expraw >= c.targetExpCov) {
  315. wins.push_back(CovWin(start, pos + 1, winlen, covsum, expcor, ucov, tcov, aall, eall, true));
  316. covsum = 0;
  317. expraw = 0;
  318. expcor = 0;
  319. ucov = 0;
  320. tcov = 0;
  321. aall = 0;
  322. eall = 0;
  323. winlen = 0;
  324. start = pos + 1;
  325. }
  326. }
  327. }
  328. } else {
  329. // Fixed windows (non-overlapping)
  330. for(uint32_t start = 0; start < hdr->target_len[refIndex]; start = start + c.window_offset) {
  331. if (start + c.window_size < hdr->target_len[refIndex]) {
  332. double covsum = 0;
  333. double expcov = 0;
  334. double ucov = 0;
  335. double tcov = 0;
  336. double aall = 0;
  337. double eall = 0;
  338. uint32_t winlen = 0;
  339. for(uint32_t pos = start; pos < start + c.window_size; ++pos) {
  340. ucov += covUniq[pos];
  341. tcov += covMap[pos];
  342. if ((gcContent[pos] > gcbound.first) && (gcContent[pos] < gcbound.second)) {
  343. aall += covAll[pos];
  344. eall += gcbias[gcContent[pos]].coverageTotal;
  345. }
  346. if (_posUsed(dt, c, gcContent, uniqContent, gcbound, pos)) {
  347. covsum += cov[pos];
  348. expcov += _expCov(dt, gcbias, gcContent[pos]) * (tileFac.empty() ? 1.0 : (double) tileFac[pos / regWin]);
  349. ++winlen;
  350. }
  351. }
  352. bool valid = (winlen >= c.fracWindow * c.window_size);
  353. wins.push_back(CovWin(start, start + c.window_size, winlen, covsum, expcov, ucov, tcov, aall, eall, valid));
  354. }
  355. }
  356. }
  358. // Separate true hom. dels from unmappable
  359. uint32_t nw = wins.size();
  360. std::vector<bool> naFlag(nw, false);
  361. std::vector<bool> suspect(nw, false);
  362. std::vector<bool> strong(nw, false);
  363. double lowFrac = 0.1; // close to CN0
  364. double flankFrac = 0.5; // clear no CN0
  365. uint32_t maxHomDelWin = 1000000; // max size for hom. DEL
  366. for(uint32_t i = 0; i < nw; ++i) {
  367. if ((!wins[i].valid) || (wins[i].expcov <= 0)) {
  368. naFlag[i] = true;
  369. continue;
  370. }
  371. double r = wins[i].covsum / wins[i].expcov;
  372. suspect[i] = (r < lowFrac);
  373. strong[i] = (r >= flankFrac);
  374. }
  375. for(uint32_t i = 0; i < nw; ) {
  376. if (naFlag[i] || (!suspect[i])) {
  377. ++i;
  378. continue;
  379. }
  380. uint32_t a = i;
  381. uint32_t b = i;
  382. while ((b + 1 < nw) && (!naFlag[b+1]) && suspect[b+1]) ++b;
  383. uint32_t runBp = wins[b].end - wins[a].start;
  384. bool leftStrong = (a > 0) && (!naFlag[a-1]) && strong[a-1];
  385. bool rightStrong = (b + 1 < nw) && (!naFlag[b+1]) && strong[b+1];
  386. bool keepDel = leftStrong && rightStrong && (runBp <= maxHomDelWin);
  387. if (!keepDel) {
  388. for(uint32_t k = a; k <= b; ++k) naFlag[k] = true;
  389. }
  390. i = b + 1;
  391. }
  393. // Total-depth deletion signal
  394. double delRatio = c.cnvDelRatio;
  396. // Flag non-unique windows
  397. bool uniqGate = c.basecov;
  398. if (uniqGate) {
  399. for(uint32_t i = 0; i < nw; ++i) {
  400. if (naFlag[i]) continue;
  401. double rTot = (wins[i].eall > 0) ? (wins[i].aall / wins[i].eall) : 1.0;
  402. if (rTot < delRatio) continue;
  403. if ((wins[i].tcov > 0) && (wins[i].ucov <= c.uniqueToTotalCovRatio * wins[i].tcov)) naFlag[i] = true;
  404. }
  405. }
  407. // Flag low callable windows
  408. if ((c.adaptive) && (nw > 4)) {
  409. std::vector<int32_t> spans(nw);
  410. for(uint32_t i = 0; i < nw; ++i) spans[i] = (int32_t) (wins[i].end - wins[i].start);
  411. std::vector<int32_t> tmp(spans);
  412. std::sort(tmp.begin(), tmp.end());
  413. int32_t medSpan = tmp[tmp.size() / 2];
  414. for(uint32_t i = 0; i < nw; ++i) tmp[i] = std::abs(spans[i] - medSpan);
  415. std::sort(tmp.begin(), tmp.end());
  416. int32_t madSpan = tmp[tmp.size() / 2];
  417. int32_t maxSpan = std::max(medSpan + (int32_t) c.mad * madSpan, 2 * medSpan);
  418. for(uint32_t i = 0; i < nw; ++i) {
  419. if (spans[i] > maxSpan) {
  420. double rTot = (wins[i].eall > 0) ? (wins[i].aall / wins[i].eall) : 1.0;
  421. if (rTot >= delRatio) naFlag[i] = true;
  422. }
  423. }
  424. }
  426. // Exclude NA windows
  427. std::vector<std::pair<int32_t, int32_t> > naiv;
  428. for(uint32_t i = 0; i < nw; ++i) {
  429. if (naFlag[i]) {
  430. for(uint32_t k = wins[i].start; k < wins[i].end; ++k) uniqContent[k] = 0;
  431. int32_t nas = (int32_t) wins[i].start;
  432. int32_t nae = (int32_t) wins[i].end;
  433. if ((!naiv.empty()) && (naiv.back().second >= nas)) naiv.back().second = std::max(naiv.back().second, nae);
  434. else naiv.push_back(std::make_pair(nas, nae));
  435. }
  436. }
  438. // Split-read breakpoints
  439. std::vector<SVBreakpoint> chrbp;
  440. collectBreakpoints(c, gcbound, gcContent, uniqContent, gcbias, cov, hdr, refIndex, clips, chrbp);
  442. // CNV discovery
  443. if (!c.hasGenoFile) segmentRD(c, gcbound, gcContent, uniqContent, gcbias, tileFac, regWin, cov, hdr, refIndex, chrbp, uniqueTrack(), naiv, cnvs);
  445. // CNV genotyping
  446. genotypeCNVs(c, gcbound, gcContent, uniqContent, gcbias, tileFac, regWin, cov, covUniq, covMap, covAll, ref, hdr, refIndex, cnvs);
  447. if (ref != NULL) free(ref);
  449. // Write windows
  450. if (!c.covfile.empty()) {
  451. std::string chrn(hdr->target_name[refIndex]);
  452. for(uint32_t i = 0; i < nw; ++i) {
  453. double uniqFrac;
  454. if (uniqGate) uniqFrac = (wins[i].tcov > 0) ? (wins[i].ucov / wins[i].tcov) : -1.0;
  455. else uniqFrac = (wins[i].end > wins[i].start) ? ((double) wins[i].winlen / (double) (wins[i].end - wins[i].start)) : -1.0;
  456. if (naFlag[i]) {
  457. dataOut << chrn << "\t" << wins[i].start << "\t" << wins[i].end << "\t" << uniqFrac << "\tNA\tNA" << std::endl;
  458. } else {
  459. double cn = chrPloidy;
  460. double logR = 0;
  461. if (wins[i].expcov > 0) {
  462. cn = (c.expectedCN * wins[i].covsum / wins[i].expcov - chrCtrlPloidy * (1 - c.purity)) / c.purity;
  463. logR = std::log2((wins[i].covsum + 1.0) / (wins[i].expcov + 1.0));
  464. }
  465. dataOut << chrn << "\t" << wins[i].start << "\t" << wins[i].end << "\t" << uniqFrac << "\t" << logR << "\t" << cn << std::endl;
  466. }
  467. }
  468. }
  469. }
  471. // Sort CNVs
  472. sort(cnvs.begin(), cnvs.end());
  474. // Merge CNVs
  475. if (!c.hasGenoFile) mergeAdjacentSameCN(cnvs, c.cnMergeTol);
  477. // Exclude regions
  478. if (c.hasExcludeFile) {
  479. typedef boost::icl::interval_set<uint32_t> TChrIntervals;
  480. typedef TChrIntervals::interval_type TIVal;
  481. std::vector<TChrIntervals> validRegions;
  482. if (_parseExcludeIntervals(c, hdr, validRegions)) {
  483. std::vector<CNV> keep;
  484. keep.reserve(cnvs.size());
  485. for(uint32_t i = 0; i < cnvs.size(); ++i) {
  486. int32_t s = cnvs[i].start;
  487. int32_t e = cnvs[i].end;
  488. if ((e <= s) || (cnvs[i].chr < 0) || (cnvs[i].chr >= (int32_t) validRegions.size())) { keep.push_back(cnvs[i]); continue; }
  489. TChrIntervals ci;
  490. ci.insert(TIVal::right_open((uint32_t) s, (uint32_t) e));
  491. TChrIntervals inter = ci & validRegions[cnvs[i].chr];
  492. uint32_t validbp = 0;
  493. for(TChrIntervals::const_iterator it = inter.begin(); it != inter.end(); ++it) validbp += (it->upper() - it->lower());
  494. if ((double) validbp >= 0.5 * (double) (e - s)) keep.push_back(cnvs[i]);
  495. }
  496. cnvs.swap(keep);
  497. }
  498. }
  500. // Genotype CNVs
  501. cnvVCF(c, cnvs);
  503. // clean-up
  504. fai_destroy(faiRef);
  505. bam_hdr_destroy(hdr);
  506. hts_idx_destroy(idx);
  507. sam_close(samfile);
  508. if (!c.covfile.empty()) {
  509. dataOut.pop();
  510. dataOut.pop();
  511. }
  513. return 0;
  514. }
  517. int coral(int argc, char **argv) {
  518. CountDNAConfig c;
  519. std::string haploidChr;
  520. std::string mode;
  522. // Parameter
  523. boost::program_options::options_description generic("Generic options");
  524. generic.add_options()
  525. ("help,?", "show help message")
  526. ("genome,g", boost::program_options::value<boost::filesystem::path>(&c.genome), "genome file")
  527. ("exclude,x", boost::program_options::value<boost::filesystem::path>(&c.exclude), "file with regions to exclude")
  528. ("quality,q", boost::program_options::value<uint16_t>(&c.minQual)->default_value(10), "min. mapping quality")
  529. ("outfile,o", boost::program_options::value<boost::filesystem::path>(&c.outfile), "BCF output file")
  530. ("covfile,c", boost::program_options::value<boost::filesystem::path>(&c.covfile), "gzipped coverage file")
  531. ("segmentation,u", boost::program_options::value<boost::filesystem::path>(&c.segfile), "segmentation BED output file")
  532. ;
  534. boost::program_options::options_description cnv("CNV calling");
  535. cnv.add_options()
  536. ("cnv-size,z", boost::program_options::value<uint32_t>(&c.minCnvSize)->default_value(1000), "min. CNV size")
  537. ("vcffile,v", boost::program_options::value<boost::filesystem::path>(&c.genofile), "input CNV BCF file for re-genotyping")
  538. ("minclip", boost::program_options::value<uint32_t>(&c.minClip)->default_value(25), "min. clipping length")
  539. ("minrefsep", boost::program_options::value<uint32_t>(&c.minRefSep)->default_value(30), "min. reference separation")
  540. ("min-bp-support", boost::program_options::value<uint32_t>(&c.minBpSupport)->default_value(3), "min. split-read support")
  541. ("penalty", boost::program_options::value<float>(&c.penalty)->default_value(3), "segmentation penalty")
  542. ("cnv-merge", boost::program_options::value<float>(&c.cnMergeTol)->default_value(0.25), "min. log2 ratio to separate CNVs")
  543. ("cnv-qual", boost::program_options::value<int32_t>(&c.cnvMinQual)->default_value(5), "min. quality for PASS")
  544. ;
  546. boost::program_options::options_description cancer("Ploidy/purity correction");
  547. cancer.add_options()
  548. ("mode,m", boost::program_options::value<std::string>(&mode)->default_value("germline"), "CNV calling mode [germline|somatic]")
  549. ("ploidy,y", boost::program_options::value<float>(&c.ploidy)->default_value(2), "sample ploidy")
  550. ("purity,p", boost::program_options::value<float>(&c.purity)->default_value(1), "sample purity [0.1, 1]")
  551. ("ctrl-ploidy", boost::program_options::value<float>(&c.ctrlPloidy)->default_value(2), "control ploidy")
  552. ("haploid-chr", boost::program_options::value<std::string>(&haploidChr), "haploid chromosomes, e.g. chrX,chrY")
  553. ;
  555. boost::program_options::options_description window("Read-depth windows");
  556. window.add_options()
  557. ("window,w", boost::program_options::value<uint32_t>(&c.window_size)->default_value(0), "window size in bp (0: automatic)")
  558. ("fraction-unique", boost::program_options::value<float>(&c.uniqueToTotalCovRatio)->default_value(0.8), "uniqueness filter [0,1]")
  559. ("cnv-del-confirm", boost::program_options::value<float>(&c.cnvDelConfirm)->default_value(1.5), "max. total-depth CN for DEL")
  560. ("cnv-del-ratio", boost::program_options::value<float>(&c.cnvDelRatio)->default_value(0.75), "total-depth ratio")
  561. ("basecov", "force base-level counting")
  562. ("fragmentcov", "force fragment-level counting")
  563. ("no-regional-gc", "disable broad GC correction")
  564. ;
  566. boost::program_options::options_description hidden("Hidden options");
  567. hidden.add_options()
  568. ("input-file", boost::program_options::value<boost::filesystem::path>(&c.bamFile), "input BAM/CRAM file")
  569. ("fragment", boost::program_options::value<float>(&c.fragmentUnique)->default_value(0.97), "min. fragment uniqueness [0,1]")
  570. ("statsfile", boost::program_options::value<boost::filesystem::path>(&c.statsFile), "gzipped stats output file (optional)")
  571. ("window-offset", boost::program_options::value<uint32_t>(&c.window_offset)->default_value(0), "window offset (0: window size)")
  572. ("fraction-window", boost::program_options::value<float>(&c.fracWindow)->default_value(0.25), "min. callable window fraction [0,1]")
  573. ("mapq-uniq", boost::program_options::value<uint16_t>(&c.mapqUniq)->default_value(20), "min. MAPQ for a uniquely-placed read")
  574. ("target-reads", boost::program_options::value<uint32_t>(&c.targetReads)->default_value(150), "target reads/window")
  575. ("min-windows", boost::program_options::value<uint32_t>(&c.minSegWin)->default_value(2), "min. windows per segment")
  576. ("scan-window", boost::program_options::value<uint32_t>(&c.scanWindow)->default_value(10000), "GC scanning window size")
  577. ("scan-regions", boost::program_options::value<boost::filesystem::path>(&c.scanFile), "GC scanning regions in BED format")
  578. ("mad-cutoff", boost::program_options::value<uint16_t>(&c.mad)->default_value(3), "median + 3 * mad count cutoff")
  579. ("percentile", boost::program_options::value<float>(&c.exclgc)->default_value(0.0005), "excl. extreme GC fraction")
  580. ("no-window-selection", "no scan window selection")
  581. ;
  583. boost::program_options::positional_options_description pos_args;
  584. pos_args.add("input-file", -1);
  586. // Set the visibility
  587. boost::program_options::options_description cmdline_options;
  588. cmdline_options.add(generic).add(cnv).add(cancer).add(window).add(hidden);
  589. boost::program_options::options_description visible_options;
  590. visible_options.add(generic).add(cnv).add(cancer).add(window);
  592. // Parse command-line
  593. boost::program_options::variables_map vm;
  594. boost::program_options::store(boost::program_options::command_line_parser(argc, argv).options(cmdline_options).positional(pos_args).run(), vm);
  595. boost::program_options::notify(vm);
  597. // Check command line arguments
  598. if ((vm.count("help")) || (!vm.count("input-file")) || (!vm.count("genome"))) {
  599. std::cerr << std::endl;
  600. std::cerr << "Usage: delly " << argv[0] << " [OPTIONS] -g <genome.fa> <aligned.bam>" << std::endl;
  601. std::cerr << visible_options << "\n";
  602. return 1;
  603. }
  605. // Show cmd
  606. boost::posix_time::ptime now = boost::posix_time::second_clock::local_time();
  607. std::cerr << '[' << boost::posix_time::to_simple_string(now) << "] ";
  608. std::cerr << "delly ";
  609. for(int i=0; i<argc; ++i) { std::cerr << argv[i] << ' '; }
  610. std::cerr << std::endl;
  612. // Stats file
  613. if (vm.count("statsfile")) c.hasStatsFile = true;
  614. else c.hasStatsFile = false;
  616. // Scan regions
  617. if (vm.count("scan-regions")) c.hasScanFile = true;
  618. else c.hasScanFile = false;
  620. // Exclude regions
  621. if (vm.count("exclude")) c.hasExcludeFile = true;
  622. else c.hasExcludeFile = false;
  624. // Scan window selection
  625. if (vm.count("no-window-selection")) c.noScanWindowSelection = true;
  626. else c.noScanWindowSelection = false;
  627. c.regionalGc = (!vm.count("no-regional-gc"));
  629. // Adaptive windows
  630. c.adaptive = false;
  631. c.targetExpCov = 0;
  632. if (c.window_size == 0) c.adaptive = true;
  633. if (c.targetReads == 0) c.targetReads = 150;
  635. // Calling mode
  636. c.somatic = (mode == "somatic");
  637. if (c.somatic) {
  638. // More sensitive segmentation
  639. if (vm["penalty"].defaulted()) c.penalty = 1.25;
  640. if (vm["cnv-merge"].defaulted()) c.cnMergeTol = 0;
  641. }
  643. // Purity
  644. if (c.purity > 1) c.purity = 1;
  645. if (c.purity < 0.1) c.purity = 0.1;
  647. // Expected ploidy
  648. c.expectedCN = c.purity * c.ploidy + (1.0 - c.purity) * c.ctrlPloidy;
  650. // Segmentation
  651. if (vm.count("segmentation")) c.hasSegFile = true;
  652. else c.hasSegFile = false;
  654. // Window offset
  655. if ((c.window_offset == 0) || (c.window_offset > c.window_size)) c.window_offset = c.window_size;
  657. // Check input VCF file (CNV genotyping)
  658. if (vm.count("vcffile")) {
  659. if (!(boost::filesystem::exists(c.genofile) && boost::filesystem::is_regular_file(c.genofile) && boost::filesystem::file_size(c.genofile))) {
  660. std::cerr << "Input VCF/BCF file is missing: " << c.genofile.string() << std::endl;
  661. return 1;
  662. }
  663. htsFile* ifile = bcf_open(c.genofile.string().c_str(), "r");
  664. if (ifile == NULL) {
  665. std::cerr << "Fail to open file " << c.genofile.string() << std::endl;
  666. return 1;
  667. }
  668. bcf_hdr_t* hdr = bcf_hdr_read(ifile);
  669. if (hdr == NULL) {
  670. std::cerr << "Fail to open index file " << c.genofile.string() << std::endl;
  671. return 1;
  672. }
  673. bcf_hdr_destroy(hdr);
  674. bcf_close(ifile);
  675. c.hasGenoFile = true;
  676. } else c.hasGenoFile = false;
  678. // Check outfile
  679. if (!vm.count("outfile")) c.outfile = "-";
  680. else {
  681. if (c.outfile.string() != "-") {
  682. if (!_outfileValid(c.outfile)) return 1;
  683. }
  684. }
  686. // Check bam file
  687. LibraryInfo li;
  688. bool pairedLib = false;
  689. if (!(boost::filesystem::exists(c.bamFile) && boost::filesystem::is_regular_file(c.bamFile) && boost::filesystem::file_size(c.bamFile))) {
  690. std::cerr << "Alignment file is missing: " << c.bamFile.string() << std::endl;
  691. return 1;
  692. } else {
  693. // Get scan regions
  694. typedef boost::icl::interval_set<uint32_t> TChrIntervals;
  695. typedef typename TChrIntervals::interval_type TIVal;
  696. typedef std::vector<TChrIntervals> TRegionsGenome;
  697. TRegionsGenome scanRegions;
  699. // Open BAM file
  700. samFile* samfile = sam_open(c.bamFile.string().c_str(), "r");
  701. if (samfile == NULL) {
  702. std::cerr << "Fail to open file " << c.bamFile.string() << std::endl;
  703. return 1;
  704. }
  705. hts_idx_t* idx = sam_index_load(samfile, c.bamFile.string().c_str());
  706. if (idx == NULL) {
  707. if (bam_index_build(c.bamFile.string().c_str(), 0) != 0) {
  708. std::cerr << "Fail to open index for " << c.bamFile.string() << std::endl;
  709. return 1;
  710. }
  711. }
  712. bam_hdr_t* hdr = sam_hdr_read(samfile);
  713. if (hdr == NULL) {
  714. std::cerr << "Fail to open header for " << c.bamFile.string() << std::endl;
  715. return 1;
  716. }
  717. c.nchr = hdr->n_targets;
  718. c.minChrLen = setMinChrLen(hdr, 0.95);
  719. std::string sampleName = "unknown";
  720. getSMTag(std::string(hdr->text), c.bamFile.stem().string(), sampleName);
  721. c.sampleName = sampleName;
  723. // Check special chromosomes
  724. if (haploidChr.size()) _selectedRefIndices(c, hdr, haploidChr);
  726. // Check matching chromosome names
  727. faidx_t* faiRef = fai_load(c.genome.string().c_str());
  728. uint32_t refFound = 0;
  729. for(int32_t refIndex=0; refIndex < hdr->n_targets; ++refIndex) {
  730. std::string tname(hdr->target_name[refIndex]);
  731. if (faidx_has_seq(faiRef, tname.c_str())) ++refFound;
  732. else {
  733. std::cerr << "Warning: BAM chromosome " << tname << " not present in reference genome!" << std::endl;
  734. }
  735. }
  736. fai_destroy(faiRef);
  737. if (!refFound) {
  738. std::cerr << "Reference genome chromosome naming disagrees with BAM file!" << std::endl;
  739. return 1;
  740. }
  742. // Estimate library params
  743. if (c.hasScanFile) {
  744. if (!_parseBedIntervals(c.scanFile.string(), c.hasScanFile, hdr, scanRegions)) {
  745. std::cerr << "Warning: Couldn't parse BED intervals. Do the chromosome names match?" << std::endl;
  746. return 1;
  747. }
  748. } else {
  749. scanRegions.resize(hdr->n_targets);
  750. for (int32_t refIndex = 0; refIndex < hdr->n_targets; ++refIndex) {
  751. scanRegions[refIndex].insert(TIVal::right_open(0, hdr->target_len[refIndex]));
  752. }
  753. }
  754. typedef std::vector<LibraryInfo> TSampleLibrary;
  755. TSampleLibrary sampleLib(1, LibraryInfo());
  756. CountDNAConfigLib dellyConf;
  757. dellyConf.genome = c.genome;
  758. dellyConf.files.push_back(c.bamFile);
  759. dellyConf.madCutoff = 9;
  760. dellyConf.madNormalCutoff = c.mad;
  761. getLibraryParams(dellyConf, scanRegions, sampleLib);
  762. li = sampleLib[0];
  763. pairedLib = (li.median > 0);
  764. if (!li.median) {
  765. li.median = 250;
  766. li.mad = 15;
  767. li.minNormalISize = 0;
  768. li.maxNormalISize = 400;
  769. }
  770. c.meanisize = ((int32_t) (li.median / 2)) * 2 + 1;
  772. // Coverage-aware GC scan window
  773. bool autoScanWin = (!vm.count("scan-window")) || (vm["scan-window"].defaulted());
  774. if ((autoScanWin) && (idx != NULL)) {
  775. uint64_t totalMapped = 0;
  776. uint64_t genomeLen = 0;
  777. for(int32_t refIndex = 0; refIndex < hdr->n_targets; ++refIndex) {
  778. uint64_t mapped = 0;
  779. uint64_t unmapped = 0;
  780. hts_idx_get_stat(idx, refIndex, &mapped, &unmapped);
  781. if (mapped > 0) {
  782. totalMapped += mapped;
  783. genomeLen += hdr->target_len[refIndex];
  784. }
  785. }
  786. if (pairedLib) totalMapped /= 2; // fragments
  787. if ((totalMapped > 0) && (genomeLen > 0)) {
  788. double fragPerBp = (double) totalMapped / (double) genomeLen;
  789. uint32_t targetScanReads = 30;
  790. uint32_t autoScan = (uint32_t) (targetScanReads / fragPerBp);
  791. if (autoScan < c.scanWindow) autoScan = c.scanWindow;
  792. uint32_t maxScan = 1000000;
  793. if (autoScan > maxScan) autoScan = maxScan;
  794. if (autoScan > c.scanWindow) c.scanWindow = autoScan;
  795. }
  796. }
  798. // Clean-up
  799. bam_hdr_destroy(hdr);
  800. hts_idx_destroy(idx);
  801. sam_close(samfile);
  802. }
  804. // Counting model
  805. if (vm.count("basecov")) c.basecov = true;
  806. else if (vm.count("fragmentcov")) c.basecov = false;
  807. else c.basecov = ((!pairedLib) && (li.rs >= 500));
  808. if ((c.basecov) && (!vm.count("cnv-del-confirm"))) c.cnvDelConfirm = std::numeric_limits<float>::max();
  810. // GC bias estimation
  811. typedef std::pair<uint32_t, uint32_t> TGCBound;
  812. TGCBound gcbound;
  813. std::vector<GcBias> gcbias(c.meanisize + 1, GcBias());
  814. typedef std::vector<ScanWindow> TWindowCounts;
  815. typedef std::vector<TWindowCounts> TGenomicWindowCounts;
  816. TGenomicWindowCounts scanCounts(c.nchr, TWindowCounts());
  817. {
  818. // Scan genomic windows
  819. scan(c, li, scanCounts);
  821. // Check coverage
  822. {
  823. std::vector<uint32_t> sampleScanVec;
  824. for(uint32_t i = 0; i < scanCounts.size(); ++i) {
  825. for(uint32_t j = 0; j < scanCounts[i].size(); ++j) {
  826. sampleScanVec.push_back(scanCounts[i][j].cov);
  827. }
  828. if (sampleScanVec.size() > 1000000) break;
  829. }
  830. std::sort(sampleScanVec.begin(), sampleScanVec.end());
  831. if (sampleScanVec.empty()) {
  832. std::cerr << "Not enough windows!" << std::endl;
  833. return 1;
  834. }
  835. if (sampleScanVec[sampleScanVec.size()/2] < 5) {
  836. std::cerr << "Coverage in the GC scan window is too low." << std::endl;
  837. return 1;
  838. }
  839. }
  841. // Select stable windows
  842. selectWindows(c, scanCounts);
  844. // Estimate GC bias
  845. gcBias(c, scanCounts, li, gcbias, gcbound);
  847. // Statistics output
  848. if (c.hasStatsFile) {
  849. // Open stats file
  850. boost::iostreams::filtering_ostream statsOut;
  851. statsOut.push(boost::iostreams::gzip_compressor());
  852. statsOut.push(boost::iostreams::file_sink(c.statsFile.string(), std::ios_base::out | std::ios_base::binary));
  854. // Library Info
  855. statsOut << "LP\t" << li.rs << ',' << li.median << ',' << li.mad << ',' << li.minNormalISize << ',' << li.maxNormalISize << std::endl;
  857. // Scan window summry
  858. samFile* samfile = sam_open(c.bamFile.string().c_str(), "r");
  859. bam_hdr_t* hdr = sam_hdr_read(samfile);
  860. statsOut << "SW\tchrom\tstart\tend\tselected\tcoverage\tuniqcov" << std::endl;
  861. for(uint32_t refIndex = 0; refIndex < (uint32_t) hdr->n_targets; ++refIndex) {
  862. for(uint32_t i = 0; i < scanCounts[refIndex].size(); ++i) {
  863. statsOut << "SW\t" << hdr->target_name[refIndex] << '\t' << scanCounts[refIndex][i].start << '\t' << scanCounts[refIndex][i].end << '\t' << scanCounts[refIndex][i].select << '\t' << scanCounts[refIndex][i].cov << '\t' << scanCounts[refIndex][i].uniqcov << std::endl;
  864. }
  865. }
  866. bam_hdr_destroy(hdr);
  867. sam_close(samfile);
  869. // GC bias summary
  870. statsOut << "GC\tgcsum\tsample\treference\tpercentileSample\tpercentileReference\tfractionSample\tfractionReference\tobsexp\tmeancoverage\tmeancoverageTotal\treferenceTotal" << std::endl;
  871. for(uint32_t i = 0; i < gcbias.size(); ++i) statsOut << "GC\t" << i << "\t" << gcbias[i].sample << "\t" << gcbias[i].reference << "\t" << gcbias[i].percentileSample << "\t" << gcbias[i].percentileReference << "\t" << gcbias[i].fractionSample << "\t" << gcbias[i].fractionReference << "\t" << gcbias[i].obsexp << "\t" << gcbias[i].coverage << "\t" << gcbias[i].coverageTotal << "\t" << gcbias[i].referenceTotal << std::endl;
  872. statsOut << "BoundsGC\t" << gcbound.first << "," << gcbound.second << std::endl;
  873. statsOut.pop();
  874. statsOut.pop();
  875. }
  876. }
  878. // Coverage-aware window size
  879. uint32_t effWin = (c.window_size > 0) ? c.window_size : 50000;
  880. if (c.adaptive) {
  881. // Mean expected coverage at CN2 over the callable GC range
  882. double covMean = 0;
  883. uint64_t refCnt = 0;
  884. for(uint32_t i = gcbound.first + 1; i < gcbound.second; ++i) {
  885. covMean += gcbias[i].coverage * (double) gcbias[i].reference;
  886. refCnt += gcbias[i].reference;
  887. }
  888. if (refCnt) covMean /= (double) refCnt;
  889. if (covMean <= 0) {
  890. // Fixed windows
  891. c.adaptive = false;
  892. c.window_size = 10000;
  893. c.window_offset = c.window_size;
  894. } else {
  895. double readLen = (li.rs > 0) ? (double) li.rs : (double) c.meanisize;
  896. double molPerBp = c.basecov ? (covMean / readLen) : covMean;
  897. if (molPerBp <= 0) molPerBp = 1e-9;
  898. double winBp = (double) c.targetReads / molPerBp;
  899. double minWin = std::max(100.0, 4.0 * readLen);
  900. double maxWin = 2000000.0;
  901. if (winBp < minWin) winBp = minWin;
  902. if (winBp > maxWin) winBp = maxWin;
  903. c.targetExpCov = covMean * winBp;
  904. effWin = (uint32_t) winBp;
  905. double effReads = molPerBp * winBp;
  906. double covDepth = c.basecov ? covMean : (covMean * readLen);
  907. now = boost::posix_time::second_clock::local_time();
  908. std::cerr << '[' << boost::posix_time::to_simple_string(now) << "] " << "Auto window size: " << (uint32_t) winBp << " bp, " << (uint32_t) effReads << " reads/window (" << (c.basecov ? "base-level" : "fragment") << ", coverage " << std::fixed << std::setprecision(2) << covDepth << std::defaultfloat << "x)" << std::endl;
  909. }
  910. }
  912. // Regional GC correction
  913. std::vector<double> regcorr;
  914. uint32_t regWin = std::max((uint32_t) 50000, effWin);
  915. if (c.regionalGc) estimateRegionalGc(c, gcbound, gcbias, scanCounts, regWin, regcorr);
  917. // Count reads
  918. if (bamCount(c, li, gcbias, gcbound, regcorr, regWin)) {
  919. std::cerr << "Read counting error!" << std::endl;
  920. return 1;
  921. }
  923. // Done
  924. now = boost::posix_time::second_clock::local_time();
  925. std::cerr << '[' << boost::posix_time::to_simple_string(now) << "] " << "Done." << std::endl;
  927. return 0;
  928. }
  931. }
  933. #endif

coral.h at commit c6474a2, under BSD-3-Clause · at the source

Overview

Authors: Lorenzo Corazzi1,2, Alex Ing1, Eva Benito3, Marco Raffaele Cosenza3, Patrick Hasenfeld3, Thomas Weber4, Anna J M Marx1,2, Vivien S Ionasz1, Nathan Trausch1,2, Sarah Benedetto5, Giulia Di Muzio1,2, Boyu Ding1,6, Jana Berlanda1,2, Marco Giaisi1, Nina Claudino5, Thomas Höfer5, Jan O Korbel3,7, Pei-Chi Wei1,2
ORCID iDs: Lorenzo Corazzi, Marco Raffaele Cosenza, Vivien S Ionasz, Boyu Ding, Jana Berlanda, Thomas Höfer, Pei-Chi Wei
  1. Brain Mosaicism and Tumorigenesis Laboratory, German Cancer Research Center, Heidelberg, Germany
  2. Faculty of Bioscience, Ruprecht-Karl-University of Heidelberg, Heidelberg, Germany
  3. European Molecular Biology Laboratory (EMBL), Genome Biology Unit, Heidelberg, Germany
  4. Data Science Centre, European Molecular Biology Laboratory (EMBL), Heidelberg, Germany
  5. Division of Theoretical Systems Biology, German Cancer Research Center, Heidelberg, Germany
  6. Faculty of Medicine, Ruprecht-Karl-University of Heidelberg, Heidelberg, Germany
  7. Bridging Research Division on Mechanisms of Genomic Variation and Data Science, German Cancer Research Center, Heidelberg, Germany
Institutions: German Cancer Research Center (Germany); Heidelberg University (Germany); European Molecular Biology Laboratory (Germany); European Molecular Biology Laboratory (Germany)
Journal: Nature communications, volume 17, issue 1, article 3627
Dates: received 20 August 2025; accepted 25 March 2026; published online 20 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-71790-5 · PMID 42009664 · PMCID PMC13096501 · OpenAlex W7154939394
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), mouse (organism), cellular / molecular (subfield)
Methods: Statistics, Preprocessing
Keywords: Genomic instability, Genomics, DNA damage and repair, DNA replication
MeSH: DNA Breaks*, DNA Copy Number Variations*, DNA Replication*, Animals, CRISPR-Cas Systems, DNA Polymerase theta, DNA Repair, DNA-Directed DNA Polymerase, Genome, Mice, Neural Stem Cells, Whole Genome Sequencing (* major topic)
Topic: Genomic variations and chromosomal abnormalities (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: European Research Council (101098056, 949990); Helmholtz Association (YIP-DKFZ)
Citations: cited by 3 papers (Europe PMC); 63 references in the paper

Abstract

Copy number variants (CNVs) are strongly implicated in neurological and psychiatric disorders and brain cancer, yet the process by which replication stress generates CNVs—and why some recur while others remain rare—remains poorly understood. Here, we show that recurrent DNA-break clusters (RDCs) act as common initiating lesions that drive both recurrent and non-recurrent CNVs. In murine neural progenitor cells subjected to chemically induced replication stress, bulk whole-genome sequencing identifies recurrent CNVs enriched at late-replicating RDCs within actively transcribed genes. Single-cell genome sequencing further uncovers frequent, non-recurrent CNVs associated with RDCs that arise during the transition from early to late DNA replication. These CNVs represent stable, heritable structural variants with breakpoints consistently enriched at RDCs. CRISPR/Cas9-mediated transcriptional suppression abolishes both RDC formation and CNV generation, establishing RDC-associated breaks as a shared upstream source. Mechanistically, CNV formation depends on DNA repair context: CNVs are Pol θ-dependent in NHEJ-deficient cells but arise independently of Pol θ in NHEJ-proficient cells. Together, these findings define RDCs as central drivers of replication-stress-induced genome diversification.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

DKFZ-ODCF/AlignmentAndQCWorkflows

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 2701b22ad5391ce53747345fff80b87f521cc60f, 9 May 2023
Languages: Shell (29), Java (25), Python (22), Perl (12), R (9)
Size: 148 files, 97 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (resources/analysisTools/qcPipeline/environments/conda.yml), tests, documentation
Not found: CITATION.cff, continuous integration
Tools: pysam (7 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
99 files
At the source: github.com/DKFZ-ODCF/AlignmentAndQCWorkflows/

CL-CHEN-Lab/RepliSeq

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 359c213d846061eba0a688412f203dfb05369bf1, 6 September 2021
Languages: R (13)
Size: 64 files, 13 scripts
Software Heritage: not archived
Found in: the text, “Two-fraction replication sequencing (Repli-Seq)”
Holds: README, license file, environment (DESCRIPTION), documentation, 1 notebook
Not found: CITATION.cff, tests, continuous integration
Tools: tidyverse (3 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
15 files
At the source: github.com/CL-CHEN-Lab/RepliSeq

brainbreaks/GROseq

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 724b09bd564e0cb2ea41a6855ec055c86dd35e11, 19 March 2024
Languages: Python (11)
Size: 18 files, 11 scripts
Software Heritage: not archived
Found in: the text, “Global run-on sequencing (GRO-seq)”
Holds: README, environment (Dockerfile)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (2 files), pandas (2 files), BEDTools (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
12 files
At the source: github.com/brainbreaks/GROseq

dellytools/delly

License: BSD-3-Clause
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: c6474a25a125982a85829749ab7c14c63eaf17a2, 24 September 2026
Languages: C/C++ (35), R (3), C++ (2), Python (1)
Size: 76 files, 41 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, environment (Dockerfile, Dockerfile.staticbuild, pyproject.toml, singularity/delly.def), continuous integration
Not found: CITATION.cff, tests, documentation
Tools: ggplot2 (3 files), reshape2 (2 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
43 files
At the source: github.com/dellytools/delly

friendsofstrandseq/mosaicatcher-pipeline

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 180b6591078f2a1bd28b2cdd757c6ee31c6299fc, 26 August 2026
Languages: R (111), Python (106), Shell (29), Jupyter (19), Perl (8)
Size: 517 files, 273 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, environment (.gitpod.Dockerfile, pixi.lock, pixi.toml, github-actions-runner/Dockerfile, github-actions-runner/Dockerfile-1.8.3-T2T.dockerfile, github-actions-runner/Dockerfile-1.8.3.dockerfile, github-actions-runner/Dockerfile-1.8.4.dockerfile, github-actions-runner/Dockerfile-1.8.5.dockerfile, github-actions-runner/Dockerfile-1.8.6.dockerfile, github-actions-runner/Dockerfile-1.9.0.dockerfile, github-actions-runner/Dockerfile-1.9.1.dockerfile, github-actions-runner/Dockerfile-2.0.0.dockerfile), continuous integration, documentation, 19 notebooks
Not found: CITATION.cff, tests
Tools: pandas (80 files), data.table (36 files), tidyverse (32 files), ggplot2 (26 files), NumPy (24 files), reshape2 (14 files), pysam (10 files), Matplotlib (9 files), SAMtools (7 files), pheatmap (6 files), TensorFlow (6 files), ggpubr (5 files), seaborn (5 files), cowplot (3 files), Keras (3 files), Plotly (3 files), SciPy (3 files), Snakemake (3 files), BCFtools (2 files), BEDTools (2 files), DESeq2 (2 files), UMAP (2 files), ComplexHeatmap (1 file), edgeR (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
275 files
At the source: github.com/friendsofstrandseq/mosaicatcher-pipeline

Zenodo 10843397

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (13 files), ggplot2 (7 files), data.table (6 files), circlize (2 files), reshape2 (2 files), SAMtools (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
35 files
At the source:

Zenodo 10838367

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (2 files), pandas (2 files), BEDTools (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
12 files
At the source:

Zenodo 10838365

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (4 files), reshape2 (3 files)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
5 files
At the source:

Code availability

In this manuscript, three singularity containers were applied for HTGTS (10.5281/zenodo.10843397), GRO-seq (10.5281/zenodo.10838367), and Repli-seq (10.5281/zenodo.10838365) analyses. In addition, an automated AlignmentAndQC workflow (https://github.com/DKFZ-ODCF/AlignmentAndQCWorkflows/) was applied to WGS analysis. Bowtie 2 2.5.0 used for LAM-HTGTS, Repli-seq, and GRO-seq. Samtools 1.7 was used for LAM-HTGTS, and 1.8 and above was used for Repli-seq and GRO-seq. These pipelines and software dependencies were containerized using Docker, ensuring reproducible and consistent version control across analyses. Delly 1.2.6 is a publicly available package (https://github.com/dellytools/delly). MosaicCatcher V2 (https://github.com/friendsofstrandseq/mosaicatcher-pipeline) and Strandtools (https://git.embl.org/cosenza/strandtools) are publicly available.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 8 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 484 scripts, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

The whole genome sequencing data generated in this study have been deposited in the European Genome Archive database under accession code PRJEB95922 (https://www.ebi.ac.uk/ena/browser/view/PRJEB95922). The two-fraction Repli-Seq and LAM-HTGTS generated in this study have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE305347 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE305347). The GRO-seq data generated in this study have been deposited in NCBI’s GEO database under accession code GSE305346 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE305346). Processed data were aligned to mm10 genome build. The Strand-seq data have been deposited in the European Genome Archive database under accession code PRJEB105360 [https://www.ebi.ac.uk/ena/browser/view/ PRJEB105360]. We retrieved high-resolution Repli-Seq data generated in this study from GSE257765 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE257765) (related to Fig. 2), LAM-HTGTS data from GSE233842 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE233842) [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE233842 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE233842)] (related to Fig. 2) and GSE74356 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74356) (wild-type Chr12 data, Fig. 6; re-aligned to mm10) for replotting. Source data are provided with this paper.

In this manuscript, three singularity containers were applied for HTGTS (10.5281/zenodo.10843397), GRO-seq (10.5281/zenodo.10838367), and Repli-seq (10.5281/zenodo.10838365) analyses. In addition, an automated AlignmentAndQC workflow (https://github.com/DKFZ-ODCF/AlignmentAndQCWorkflows/) was applied to WGS analysis. Bowtie 2 2.5.0 used for LAM-HTGTS, Repli-seq, and GRO-seq. Samtools 1.7 was used for LAM-HTGTS, and 1.8 and above was used for Repli-seq and GRO-seq. These pipelines and software dependencies were containerized using Docker, ensuring reproducible and consistent version control across analyses. Delly 1.2.6 is a publicly available package (https://github.com/dellytools/delly). MosaicCatcher V2 (https://github.com/friendsofstrandseq/mosaicatcher-pipeline) and Strandtools (https://git.embl.org/cosenza/strandtools) are publicly available.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 18 authors, 4 keywords, 12 MeSH terms, 2 funders, 62 references.

Cite

This paper

Corazzi, L., Ing, A., Benito, E., Cosenza, M. R., Hasenfeld, P., Weber, T., Marx, A. J. M., Ionasz, V. S., Trausch, N., Benedetto, S., Di Muzio, G., Ding, B., Berlanda, J., Giaisi, M., Claudino, N., Höfer, T., Korbel, J. O., & Wei, P.-C. (2026). Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification. Nature communications, 17(1), 3627. https://doi.org/10.1038/s41467-026-71790-5

BibTeX

@article{corazzi2026recurrent,
author = {Corazzi, Lorenzo and Ing, Alex and Benito, Eva and Cosenza, Marco Raffaele and Hasenfeld, Patrick and Weber, Thomas and Marx, Anna J M and Ionasz, Vivien S and Trausch, Nathan and Benedetto, Sarah and Di Muzio, Giulia and Ding, Boyu and Berlanda, Jana and Giaisi, Marco and Claudino, Nina and Höfer, Thomas and Korbel, Jan O and Wei, Pei-Chi},
title = {{Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {3627},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71790-5},
url = {https://doi.org/10.1038/s41467-026-71790-5},
pmid = {42009664},
pmcid = {PMC13096501}
}

RIS

TY - JOUR
AU - Corazzi, Lorenzo
AU - Ing, Alex
AU - Benito, Eva
AU - Cosenza, Marco Raffaele
AU - Hasenfeld, Patrick
AU - Weber, Thomas
AU - Marx, Anna J M
AU - Ionasz, Vivien S
AU - Trausch, Nathan
AU - Benedetto, Sarah
AU - Di Muzio, Giulia
AU - Ding, Boyu
AU - Berlanda, Jana
AU - Giaisi, Marco
AU - Claudino, Nina
AU - Höfer, Thomas
AU - Korbel, Jan O
AU - Wei, Pei-Chi
TI - Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/20
VL - 17
IS - 1
SP - 3627
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71790-5
UR - https://doi.org/10.1038/s41467-026-71790-5
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71790-5",
"type": "article-journal",
"title": "Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification",
"container-title": "Nature communications",
"author": [
{
"family": "Corazzi",
"given": "Lorenzo"
},
{
"family": "Ing",
"given": "Alex"
},
{
"family": "Benito",
"given": "Eva"
},
{
"family": "Cosenza",
"given": "Marco Raffaele"
},
{
"family": "Hasenfeld",
"given": "Patrick"
},
{
"family": "Weber",
"given": "Thomas"
},
{
"family": "Marx",
"given": "Anna J M"
},
{
"family": "Ionasz",
"given": "Vivien S"
},
{
"family": "Trausch",
"given": "Nathan"
},
{
"family": "Benedetto",
"given": "Sarah"
},
{
"family": "Di Muzio",
"given": "Giulia"
},
{
"family": "Ding",
"given": "Boyu"
},
{
"family": "Berlanda",
"given": "Jana"
},
{
"family": "Giaisi",
"given": "Marco"
},
{
"family": "Claudino",
"given": "Nina"
},
{
"family": "Höfer",
"given": "Thomas"
},
{
"family": "Korbel",
"given": "Jan O"
},
{
"family": "Wei",
"given": "Pei-Chi"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "3627",
"DOI": "10.1038/s41467-026-71790-5",
"PMID": "42009664",
"PMCID": "PMC13096501",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71790-5",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
20
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: BCFtools, SAMtools, edgeR, 18 other tools, mouse, cellular / molecular
[2] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: Snakemake, pysam, BEDTools, 17 other tools, genetics / omics, mouse, cellular / molecular
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pysam, edgeR, Keras, 18 other tools, genetics / omics, cellular / molecular
[4] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: BCFtools, SAMtools, edgeR, 17 other tools, genetics / omics
[5] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Snakemake, pysam, BEDTools, 15 other tools, genetics / omics, mouse, cellular / molecular
[6] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: pysam, BEDTools, SAMtools, 15 other tools, mouse, cellular / molecular
[7] doi:10.1038/s41586-026-10612-6 [code]
Acquired genetic and cell-state changes in IDH-mutant glioma progression.
Journal: Nature
In common: Snakemake, BCFtools, pysam, 13 other tools, cellular / molecular
[8] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: pysam, BEDTools, SAMtools, 15 other tools, genetics / omics
[9] doi:10.3389/fnmol.2026.1844705 [code]
Risperidone regulates the expression of schizophrenia-related genes in the forebrain of adult male mice.
Journal: Frontiers in molecular neuroscience
In common: pysam, BEDTools, SAMtools, 14 other tools, genetics / omics, mouse, cellular / molecular
[10] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: BCFtools, BEDTools, SAMtools, 12 other tools, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.