search

rwdavies / STITCH

STITCH - Sequencing To Imputation Through Constructing Haplotypes
http://www.nature.com/ng/journal/v48/n8/abs/ng.3594.html
GNU General Public License v3.0
77 stars 17 forks source link

stitch crashed #86

Open rmott opened 1 year ago

rmott commented 1 year ago

Hi Robbie

I gort a segfault runnning stitch. The data comprise about 550 Arabidopsis bamfiles. I divided the genome into 5Mb segments and STITCH processed the first few segments without problem but suddenly crashed on segment #4. The error was in the part where i is downsampling the bam files:

[2023-09-13 14:44:16] downsample sample A77838_1 - 353 of 122150 reads removed [2023-09-13 14:44:16] downsample sample A77800_1 - 336 of 108490 reads removed [2023-09-13 14:44:16] downsample sample A77912_1 - 358 of 108471 reads removed [2023-09-13 14:44:16] downsample sample A77764_1 - 10 of 97862 reads removed [2023-09-13 14:44:17] downsample sample A77690_1 - 17 of 133077 reads removed [2023-09-13 14:44:17] downsample sample A77949_1 - 285 of 102064 reads removed [2023-09-13 14:44:17] downsample sample A83242_1 - 169 of 264566 reads removed [2023-09-13 14:44:17] downsample sample A77967_1 - 47 of 119207 reads removed [2023-09-13 14:44:18] downsample sample A77634_1 - 531 of 127996 reads removed [2023-09-13 14:44:18] downsample sample A77783_1 - 22 of 105780 reads removed

caught segfault address 0x62ffd004, cause 'memory not mapped'

Traceback: 1: cpp_read_reassign(ord = ord, qnameInteger_ord = qnameInteger_ord, bxtagInteger_ord = bxtagInteger_ord, bxtag_bad_ord = bxtag_bad_ord, qname = qname, bxtag = bxtag, strand = strand,\ sampleReadsRaw = sampleReadsRaw, readStart_ord = readStart_ord, readEnd_ord = readEnd_ord, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpperLimit, bxTagUpperL\ imit = bxTagUpperLimit, use_bx_tag = use_bx_tag, save_sampleReadsInfo = save_sampleReadsInfo) 2: merge_reads_from_sampleReadsRaw(sampleReadsRaw = sampleReadsRaw, qname = qname, bxtag = bxtag, strand = strand, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpp\ erLimit, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit, save_sampleReadsInfo = save_sampleReadsInfo, qname_all = qname_all, readStart_all = readStart_all, readEnd_all = re\ adEnd_all) 3: loadBamAndConvert(iBam = iBam, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter =\ bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, \ tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, chrLength = chrLength, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit) 4: FUN(X[[i]], ...) 5: lapply(X = S, FUN = FUN, ...) 6: doTryCatch(return(expr), name, parentenv, handler) 7: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 8: tryCatchList(expr, classes, parentenv, handlers) 9: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcal\ l <- deparse(call)[1L] prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nc\ har(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefi\ x, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))\ ) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 10: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 11: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 12: FUN(X[[i]], ...) 13: lapply(seq_len(cores), inner.do) 14: mclapply(1:length(sampleRanges), mc.cores = nCores, FUN = loadBamAndConvert_across_a_range, sampleRanges = sampleRanges, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nS\ NPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov,\ sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, chrLen\ gth = chrLength, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit) 15: generate_input(bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLi\ mit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName\ = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, nCores = nCores, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUp\ perLimit) 16: generate_or_refactor_input(regenerateInput = regenerateInput, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, refere\ nce = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, outputdir = outputdir, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inpu\ tdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, generateInputOnly = generateInputOnly, nCo\ res = nCores, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit) 17: STITCH(chr = as.character(chr), posfile = posfile, K = 19, S = 1, nGen = 7, outputdir = outputdir, bamlist = "bamfiles.complete.txt", method = "diploid", nCores = 30, splitReadIter\ ations = NA, tempdir = "/scratch0/ucbtrfm/stitch/", regionStart = start, regionEnd = end, buffer = buffer) 18: eval(ei, envir) 19: eval(ei, envir) 20: withVisible(eval(ei, envir)) 21: source("stitch.complete.R")

Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: [2023-09-13 14:44:18] downsample sample A77561_1 - 1 of 102508 reads removed [2023-09-13 14:44:19] downsample sample A77876_1 - 23 of 113443 reads removed [2023-09-13 14:44:20] downsample sample A77618_1 - 14 of 145622 reads removed [2023-09-13 14:44:20] downsample sample A77544_1 - 7 of 149986 reads removed [2023-09-13 14:44:20] downsample sample A77986_1 - 10 of 122339 reads removed

rmott commented 1 year ago

And here is the script:

library(STITCH) options(bitmapType='cairo') chunksize = 5000000 buffer = 5000 outputdir="./tair10.complete.STITCH/" chrlen = read.delim("./tair10.chrlen.txt") rownames(chrlen) = chrlen$chr for( chr in 1:5 ) { len = chrlen$len[chr] chunks = floor(len/chunksize)+1 cat("chr", chr, "chunks", chunks, "\n") start=1 end = 0 posfile=paste0("chr",chr,".pos.txt") for( k in 1:chunks) { start = end+1 end = start + chunksize-1 if ( end > len ) end = len cat("chr", chr, "chunk", k, "start", start, "end", end, "\n") STITCH( chr=as.character(chr), posfile=posfile, K=19,S=1,nGen=7,outputdir=outputdir, bamlist="bamfiles.complete.txt",method="diploid",nCores=30,splitReadIterations = NA,tempdir="\ /scratch0/ucbtrfm/stitch/", regionStart=start, regionEnd=end, buffer=buffer) } }

rwdavies commented 1 year ago

What version is this? Is this a recent one?

On Wed, 13 Sep 2023 at 10:21, Richard Mott @.***> wrote:

And here is the script:

library(STITCH) options(bitmapType='cairo') chunksize = 5000000 buffer = 5000 outputdir="./tair10.complete.STITCH/" chrlen = read.delim("./tair10.chrlen.txt") rownames(chrlen) = chrlen$chr for( chr in 1:5 ) { len = chrlen$len[chr] chunks = floor(len/chunksize)+1 cat("chr", chr, "chunks", chunks, "\n") start=1 end = 0 posfile=paste0("chr",chr,".pos.txt") for( k in 1:chunks) { start = end+1 end = start + chunksize-1 if ( end > len ) end = len cat("chr", chr, "chunk", k, "start", start, "end", end, "\n") STITCH( chr=as.character(chr), posfile=posfile, K=19,S=1,nGen=7,outputdir=outputdir, bamlist="bamfiles.complete.txt",method="diploid",nCores=30,splitReadIterations = NA,tempdir=" /scratch0/ucbtrfm/stitch/", regionStart=start, regionEnd=end, buffer=buffer) } }

— Reply to this email directly, view it on GitHub https://github.com/rwdavies/STITCH/issues/86#issuecomment-1717732403, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7I36PXHWBNPWXPV3ZNHQLX2G6OLANCNFSM6AAAAAA4WPZGMY . You are receiving this because you are subscribed to this thread.Message ID: @.***>

rmott commented 1 year ago

The version is from this summer – it ran OK on a previous dataset. I added about 50 extra bam files to the data this week and re-ran it when it crashed. I am now re-doing it to see if it was a file system problem

From: rwdavies @.> Date: Wednesday, 13 September 2023 at 18:32 To: rwdavies/STITCH @.> Cc: Mott, Richard @.>, Author @.> Subject: Re: [rwdavies/STITCH] stitch crashed (Issue #86)

⚠ Caution: External sender

What version is this? Is this a recent one?

On Wed, 13 Sep 2023 at 10:21, Richard Mott @.***> wrote:

And here is the script:

library(STITCH) options(bitmapType='cairo') chunksize = 5000000 buffer = 5000 outputdir="./tair10.complete.STITCH/" chrlen = read.delim("./tair10.chrlen.txt") rownames(chrlen) = chrlen$chr for( chr in 1:5 ) { len = chrlen$len[chr] chunks = floor(len/chunksize)+1 cat("chr", chr, "chunks", chunks, "\n") start=1 end = 0 posfile=paste0("chr",chr,".pos.txt") for( k in 1:chunks) { start = end+1 end = start + chunksize-1 if ( end > len ) end = len cat("chr", chr, "chunk", k, "start", start, "end", end, "\n") STITCH( chr=as.character(chr), posfile=posfile, K=19,S=1,nGen=7,outputdir=outputdir, bamlist="bamfiles.complete.txt",method="diploid",nCores=30,splitReadIterations = NA,tempdir=" /scratch0/ucbtrfm/stitch/", regionStart=start, regionEnd=end, buffer=buffer) } }

— Reply to this email directly, view it on GitHub https://github.com/rwdavies/STITCH/issues/86#issuecomment-1717732403, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7I36PXHWBNPWXPV3ZNHQLX2G6OLANCNFSM6AAAAAA4WPZGMY . You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHubhttps://github.com/rwdavies/STITCH/issues/86#issuecomment-1718041720, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AADK6MVIA5TXLPF6ZKJPOGLX2HU4BANCNFSM6AAAAAA4WPZGMY. You are receiving this because you authored the thread.Message ID: @.***>

rmott commented 1 year ago

Crashed again, but earlier on:

source("stitch.complete.R") Loading required package: parallel Loading required package: rrbgen chr 1 chunks 7 chr 1 chunk 1 start 1 end 5e+06 [2023-09-13 19:31:21] Running STITCH(chr = 1, nGen = 7, posfile = chr1.pos.txt, K = 19, S = 1, outputdir = ./tair10.complete.STITCH/, nStarts = , tempdir = /scratch0/ucbtrfm/stitch/, bamlist =\ bamfiles.complete.txt, cramlist = , sampleNames_file = , reference = , genfile = , method = diploid, output_format = bgvcf, B_bit_prob = 16, outputInputInVCFFormat = FALSE, downsampleToCov = \ 50, downsampleFraction = 1, readAware = TRUE, chrStart = NA, chrEnd = NA, regionStart = 1, regionEnd = 5e+06, buffer = 5000, maxDifferenceBetweenReads = 1000, maxEmissionMatrixDifference = 1e+\ 10, alphaMatThreshold = 1e-04, emissionThreshold = 1e-04, iSizeUpperLimit = 600, bqFilter = 17, niterations = 40, shuffleHaplotypeIterations = c(4, 8, 12, 16), splitReadIterations = NA, nCores\ = 30, expRate = 0.5, maxRate = 100, minRate = 0.1, Jmax = 1000, regenerateInput = TRUE, originalRegionName = NA, keepInterimFiles = FALSE, keepTempDir = FALSE, outputHaplotypeProbabilities = \ FALSE, switchModelIteration = NA, generateInputOnly = FALSE, restartIterations = NA, refillIterations = c(6, 10, 14, 18), downsampleSamples = 1, downsampleSamplesKeepList = NA, subsetSNPsfile \ = NA, useSoftClippedBases = FALSE, outputBlockSize = 1000, outputSNPBlockSize = 10000, inputBundleBlockSize = NA, genetic_map_file = , reference_haplotype_file = , reference_legend_file = , re\ ference_sample_file = , reference_populations = NA, reference_phred = 20, reference_iterations = 40, reference_shuffleHaplotypeIterations = c(4, 8, 12, 16), output_filename = NULL, initial_min\ _hapProb = 0.2, initial_max_hapProb = 0.8, regenerateInputWithDefaultValues = FALSE, plotHapSumDuringIterations = FALSE, plot_shuffle_haplotype_attempts = FALSE, plotAfterImputation = TRUE, sa\ ve_sampleReadsInfo = FALSE, gridWindowSize = NA, shuffle_bin_nSNPs = NULL, shuffle_bin_radius = 5000, keepSampleReadsInRAM = FALSE, useTempdirWhileWriting = FALSE, output_haplotype_dosages = F\ ALSE, use_bx_tag = TRUE, bxTagUpperLimit = 50000) [2023-09-13 19:31:21] Program start [2023-09-13 19:31:21] Get and validate pos and gen [2023-09-13 19:31:23] Done get and validate pos and gen [2023-09-13 19:31:23] There are 0 variants in the left buffer region -4999 <= position < 1 [2023-09-13 19:31:23] There are 57234 variants in the central region 1 <= position <= 5000000 [2023-09-13 19:31:23] There are 60 variants in the right buffer region 5000000 < position <= 5005000 [2023-09-13 19:31:23] Get BAM sample names [2023-09-13 19:31:23] Done getting BAM sample names [2023-09-13 19:31:23] Generate inputs [2023-09-13 19:31:26] downsample sample A77890_1 - 76 of 96615 reads removed [2023-09-13 19:31:27] downsample sample A77964_1 - 8 of 117108 reads removed [2023-09-13 19:31:29] downsample sample A77614_1 - 80 of 97312 reads removed [2023-09-13 19:31:29] downsample sample A83239_1 - 34 of 228585 reads removed

…..

[2023-09-13 19:32:20] Load and convert BAM 200 of 548 [2023-09-13 19:32:20] downsample sample A83188_1 - 153 of 284786 reads removed [2023-09-13 19:32:21] downsample sample A83208_1 - 22 of 250047 reads removed [2023-09-13 19:32:21] downsample sample A77556_1 - 31 of 122470 reads removed [2023-09-13 19:32:22] downsample sample A77629_1 - 20 of 170104 reads removed [2023-09-13 19:32:22] downsample sample A77942_1 - 62 of 120575 reads removed [2023-09-13 19:32:22] downsample sample A77536_1 - 10 of 149674 reads removed [2023-09-13 19:32:24] downsample sample A83248_1 - 125 of 306130 reads removed [2023-09-13 19:32:24] downsample sample A77813_1 - 30 of 148545 reads removed [2023-09-13 19:32:27] downsample sample A83209_1 - 42 of 230660 reads removed [2023-09-13 19:32:27] downsample sample A83228_1 - 19 of 240880 reads removed [2023-09-13 19:32:27] downsample sample A83189_1 - 642 of 262446 reads removed [2023-09-13 19:32:28] downsample sample A83176_1 - 27 of 227783 reads removed [2023-09-13 19:32:28] downsample sample A77740_1 - 151 of 139589 reads removed [2023-09-13 19:32:32] downsample sample A83249_1 - 146 of 280145 reads removed [2023-09-13 19:32:33] downsample sample A83210_1 - 46 of 242464 reads removed [2023-09-13 19:32:34] downsample sample A83229_1 - 672 of 264224 reads removed [2023-09-13 19:32:38] downsample sample A83190_1 - 1826 of 342947 reads removed [2023-09-13 19:32:39] downsample sample A83250_1 - 307 of 300085 reads removed [2023-09-13 19:32:41] downsample sample A83230_1 - 623 of 256541 reads removed [2023-09-13 19:32:41] downsample sample A83211_1 - 85 of 278438 reads removed [2023-09-13 19:32:45] downsample sample A83191_1 - 131 of 280657 reads removed [2023-09-13 19:32:46] downsample sample A83212_1 - 28 of 232101 reads removed [2023-09-13 19:32:47] downsample sample A83251_1 - 81 of 256875 reads removed [2023-09-13 19:32:47] downsample sample A83231_1 - 106 of 255269 reads removed [2023-09-13 19:32:51] downsample sample A83192_1 - 44 of 251911 reads removed [2023-09-13 19:32:53] downsample sample A83232_1 - 12 of 249519 reads removed [2023-09-13 19:32:54] downsample sample A83252_1 - 130 of 268270 reads removed [2023-09-13 19:32:54] downsample sample A83213_1 - 138 of 273396 reads removed [2023-09-13 19:32:58] downsample sample A83194_1 - 15 of 245436 reads removed [2023-09-13 19:32:58] downsample sample A83233_1 - 13 of 208481 reads removed [2023-09-13 19:33:00] downsample sample A83253_1 - 32 of 243635 reads removed [2023-09-13 19:33:01] downsample sample A83214_1 - 646 of 292822 reads removed [2023-09-13 19:33:04] downsample sample A83234_1 - 4 of 236761 reads removed [2023-09-13 19:33:04] downsample sample A83195_1 - 194 of 232258 reads removed [2023-09-13 19:33:06] downsample sample A83254_1 - 25 of 215002 reads removed

caught segfault address 0x3cd49000, cause 'memory not mapped'

Traceback: 1: cpp_read_reassign(ord = ord, qnameInteger_ord = qnameInteger_ord, bxtagInteger_ord = bxtagInteger_ord, bxtag_bad_ord = bxtag_bad_ord, qname = qname, bxtag = bxtag, strand = strand,\ sampleReadsRaw = sampleReadsRaw, readStart_ord = readStart_ord, readEnd_ord = readEnd_ord, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpperLimit, bxTagUpperL\ imit = bxTagUpperLimit, use_bx_tag = use_bx_tag, save_sampleReadsInfo = save_sampleReadsInfo) 2: merge_reads_from_sampleReadsRaw(sampleReadsRaw = sampleReadsRaw, qname = qname, bxtag = bxtag, strand = strand, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpp\ erLimit, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit, save_sampleReadsInfo = save_sampleReadsInfo, qname_all = qname_all, readStart_all = readStart_all, readEnd_all = re\ adEnd_all) 3: loadBamAndConvert(iBam = iBam, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter =\ bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, \ 4: FUN(X[[i]], ...) 5: lapply(X = S, FUN = FUN, ...) 6: doTryCatch(return(expr), name, parentenv, handler) 7: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 8: tryCatchList(expr, classes, parentenv, handlers) 9: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcal\ l <- deparse(call)[1L] prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nc\ har(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefi\ x, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))\ ) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 10: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 11: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 12: FUN(X[[i]], ...) 13: lapply(seq_len(cores), inner.do) 14: mclapply(1:length(sampleRanges), mc.cores = nCores, FUN = loadBamAndConvert_across_a_range, sampleRanges = sampleRanges, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nS\ NPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov,\ sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, chrLen\ gth = chrLength, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit) 15: generate_input(bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLi\ mit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName\ = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, nCores = nCores, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUp\ perLimit) 16: generate_or_refactor_input(regenerateInput = regenerateInput, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, refere\ nce = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, outputdir = outputdir, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inpu\ 17: STITCH(chr = as.character(chr), posfile = posfile, K = 19, S = 1, nGen = 7, outputdir = outputdir, bamlist = "bamfiles.complete.txt", method = "diploid", nCores = 30, splitReadIter\ ations = NA, tempdir = "/scratch0/ucbtrfm/stitch/", regionStart = start, regionEnd = end, buffer = buffer) 18: eval(ei, envir) 19: eval(ei, envir) 20: withVisible(eval(ei, envir)) 21: source("stitch.complete.R")

Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: [2023-09-13 19:33:08] downsample sample A83215_1 - 521 of 283482 reads removed [2023-09-13 19:33:12] downsample sample A83196_1 - 154 of 263333 reads removed [2023-09-13 19:33:14] downsample sample A83216_1 - 19 of 258735 reads removed [2023-09-13 19:33:16] downsample sample A83256_1 - 2 of 196172 reads removed

From: rwdavies @.> Date: Wednesday, 13 September 2023 at 18:32 To: rwdavies/STITCH @.> Cc: Mott, Richard @.>, Author @.> Subject: Re: [rwdavies/STITCH] stitch crashed (Issue #86)

⚠ Caution: External sender

What version is this? Is this a recent one?

On Wed, 13 Sep 2023 at 10:21, Richard Mott @.***> wrote:

And here is the script:

library(STITCH) options(bitmapType='cairo') chunksize = 5000000 buffer = 5000 outputdir="./tair10.complete.STITCH/" chrlen = read.delim("./tair10.chrlen.txt") rownames(chrlen) = chrlen$chr for( chr in 1:5 ) { len = chrlen$len[chr] chunks = floor(len/chunksize)+1 cat("chr", chr, "chunks", chunks, "\n") start=1 end = 0 posfile=paste0("chr",chr,".pos.txt") for( k in 1:chunks) { start = end+1 end = start + chunksize-1 if ( end > len ) end = len cat("chr", chr, "chunk", k, "start", start, "end", end, "\n") STITCH( chr=as.character(chr), posfile=posfile, K=19,S=1,nGen=7,outputdir=outputdir, bamlist="bamfiles.complete.txt",method="diploid",nCores=30,splitReadIterations = NA,tempdir=" /scratch0/ucbtrfm/stitch/", regionStart=start, regionEnd=end, buffer=buffer) } }

— Reply to this email directly, view it on GitHub https://github.com/rwdavies/STITCH/issues/86#issuecomment-1717732403, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7I36PXHWBNPWXPV3ZNHQLX2G6OLANCNFSM6AAAAAA4WPZGMY . You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHubhttps://github.com/rwdavies/STITCH/issues/86#issuecomment-1718041720, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AADK6MVIA5TXLPF6ZKJPOGLX2HU4BANCNFSM6AAAAAA4WPZGMY. You are receiving this because you authored the thread.Message ID: @.***>

rwdavies commented 1 year ago

Can you install a recent version and double check? Just in case? I thought I had fixed an issue related to this and am surprised to see it reoccurring.

Robbie

On Wed, 13 Sep 2023 at 14:36, Richard Mott @.***> wrote:

Crashed again, but earlier on:

source("stitch.complete.R") Loading required package: parallel Loading required package: rrbgen chr 1 chunks 7 chr 1 chunk 1 start 1 end 5e+06 [2023-09-13 19:31:21] Running STITCH(chr = 1, nGen = 7, posfile = chr1.pos.txt, K = 19, S = 1, outputdir = ./tair10.complete.STITCH/, nStarts = , tempdir = /scratch0/ucbtrfm/stitch/, bamlist =\ bamfiles.complete.txt, cramlist = , sampleNames_file = , reference = , genfile = , method = diploid, output_format = bgvcf, B_bit_prob = 16, outputInputInVCFFormat = FALSE, downsampleToCov = \ 50, downsampleFraction = 1, readAware = TRUE, chrStart = NA, chrEnd = NA, regionStart = 1, regionEnd = 5e+06, buffer = 5000, maxDifferenceBetweenReads = 1000, maxEmissionMatrixDifference = 1e+\ 10, alphaMatThreshold = 1e-04, emissionThreshold = 1e-04, iSizeUpperLimit = 600, bqFilter = 17, niterations = 40, shuffleHaplotypeIterations = c(4, 8, 12, 16), splitReadIterations = NA, nCores\ = 30, expRate = 0.5, maxRate = 100, minRate = 0.1, Jmax = 1000, regenerateInput = TRUE, originalRegionName = NA, keepInterimFiles = FALSE, keepTempDir = FALSE, outputHaplotypeProbabilities = \ FALSE, switchModelIteration = NA, generateInputOnly = FALSE, restartIterations = NA, refillIterations = c(6, 10, 14, 18), downsampleSamples = 1, downsampleSamplesKeepList = NA, subsetSNPsfile \ = NA, useSoftClippedBases = FALSE, outputBlockSize = 1000, outputSNPBlockSize = 10000, inputBundleBlockSize = NA, genetic_map_file = , reference_haplotype_file = , reference_legend_file = , re\ ference_sample_file = , reference_populations = NA, reference_phred = 20, reference_iterations = 40, reference_shuffleHaplotypeIterations = c(4, 8, 12, 16), output_filename = NULL, initial_min\ _hapProb = 0.2, initial_max_hapProb = 0.8, regenerateInputWithDefaultValues = FALSE, plotHapSumDuringIterations = FALSE, plot_shuffle_haplotype_attempts = FALSE, plotAfterImputation = TRUE, sa\ ve_sampleReadsInfo = FALSE, gridWindowSize = NA, shuffle_bin_nSNPs = NULL, shuffle_bin_radius = 5000, keepSampleReadsInRAM = FALSE, useTempdirWhileWriting = FALSE, output_haplotype_dosages = F\ ALSE, use_bx_tag = TRUE, bxTagUpperLimit = 50000) [2023-09-13 19:31:21] Program start [2023-09-13 19:31:21] Get and validate pos and gen [2023-09-13 19:31:23] Done get and validate pos and gen [2023-09-13 19:31:23] There are 0 variants in the left buffer region -4999 <= position < 1 [2023-09-13 19:31:23] There are 57234 variants in the central region 1 <= position <= 5000000 [2023-09-13 19:31:23] There are 60 variants in the right buffer region 5000000 < position <= 5005000 [2023-09-13 19:31:23] Get BAM sample names [2023-09-13 19:31:23] Done getting BAM sample names [2023-09-13 19:31:23] Generate inputs [2023-09-13 19:31:26] downsample sample A77890_1 - 76 of 96615 reads removed [2023-09-13 19:31:27] downsample sample A77964_1 - 8 of 117108 reads removed [2023-09-13 19:31:29] downsample sample A77614_1 - 80 of 97312 reads removed [2023-09-13 19:31:29] downsample sample A83239_1 - 34 of 228585 reads removed

…..

[2023-09-13 19:32:20] Load and convert BAM 200 of 548 [2023-09-13 19:32:20] downsample sample A83188_1 - 153 of 284786 reads removed [2023-09-13 19:32:21] downsample sample A83208_1 - 22 of 250047 reads removed [2023-09-13 19:32:21] downsample sample A77556_1 - 31 of 122470 reads removed [2023-09-13 19:32:22] downsample sample A77629_1 - 20 of 170104 reads removed [2023-09-13 19:32:22] downsample sample A77942_1 - 62 of 120575 reads removed [2023-09-13 19:32:22] downsample sample A77536_1 - 10 of 149674 reads removed [2023-09-13 19:32:24] downsample sample A83248_1 - 125 of 306130 reads removed [2023-09-13 19:32:24] downsample sample A77813_1 - 30 of 148545 reads removed [2023-09-13 19:32:27] downsample sample A83209_1 - 42 of 230660 reads removed [2023-09-13 19:32:27] downsample sample A83228_1 - 19 of 240880 reads removed [2023-09-13 19:32:27] downsample sample A83189_1 - 642 of 262446 reads removed [2023-09-13 19:32:28] downsample sample A83176_1 - 27 of 227783 reads removed [2023-09-13 19:32:28] downsample sample A77740_1 - 151 of 139589 reads removed [2023-09-13 19:32:32] downsample sample A83249_1 - 146 of 280145 reads removed [2023-09-13 19:32:33] downsample sample A83210_1 - 46 of 242464 reads removed [2023-09-13 19:32:34] downsample sample A83229_1 - 672 of 264224 reads removed [2023-09-13 19:32:38] downsample sample A83190_1 - 1826 of 342947 reads removed [2023-09-13 19:32:39] downsample sample A83250_1 - 307 of 300085 reads removed [2023-09-13 19:32:41] downsample sample A83230_1 - 623 of 256541 reads removed [2023-09-13 19:32:41] downsample sample A83211_1 - 85 of 278438 reads removed [2023-09-13 19:32:45] downsample sample A83191_1 - 131 of 280657 reads removed [2023-09-13 19:32:46] downsample sample A83212_1 - 28 of 232101 reads removed [2023-09-13 19:32:47] downsample sample A83251_1 - 81 of 256875 reads removed [2023-09-13 19:32:47] downsample sample A83231_1 - 106 of 255269 reads removed [2023-09-13 19:32:51] downsample sample A83192_1 - 44 of 251911 reads removed [2023-09-13 19:32:53] downsample sample A83232_1 - 12 of 249519 reads removed [2023-09-13 19:32:54] downsample sample A83252_1 - 130 of 268270 reads removed [2023-09-13 19:32:54] downsample sample A83213_1 - 138 of 273396 reads removed [2023-09-13 19:32:58] downsample sample A83194_1 - 15 of 245436 reads removed [2023-09-13 19:32:58] downsample sample A83233_1 - 13 of 208481 reads removed [2023-09-13 19:33:00] downsample sample A83253_1 - 32 of 243635 reads removed [2023-09-13 19:33:01] downsample sample A83214_1 - 646 of 292822 reads removed [2023-09-13 19:33:04] downsample sample A83234_1 - 4 of 236761 reads removed [2023-09-13 19:33:04] downsample sample A83195_1 - 194 of 232258 reads removed [2023-09-13 19:33:06] downsample sample A83254_1 - 25 of 215002 reads removed

caught segfault address 0x3cd49000, cause 'memory not mapped'

Traceback: 1: cpp_read_reassign(ord = ord, qnameInteger_ord = qnameInteger_ord, bxtagInteger_ord = bxtagInteger_ord, bxtag_bad_ord = bxtag_bad_ord, qname = qname, bxtag = bxtag, strand = strand,\ sampleReadsRaw = sampleReadsRaw, readStart_ord = readStart_ord, readEnd_ord = readEnd_ord, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpperLimit, bxTagUpperL\ imit = bxTagUpperLimit, use_bx_tag = use_bx_tag, save_sampleReadsInfo = save_sampleReadsInfo) 2: merge_reads_from_sampleReadsRaw(sampleReadsRaw = sampleReadsRaw, qname = qname, bxtag = bxtag, strand = strand, readStart = readStart, readEnd = readEnd, iSizeUpperLimit = iSizeUpp\ erLimit, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit, save_sampleReadsInfo = save_sampleReadsInfo, qname_all = qname_all, readStart_all = readStart_all, readEnd_all = re\ adEnd_all) 3: loadBamAndConvert(iBam = iBam, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter =\ bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, \ 4: FUN(X[[i]], ...) 5: lapply(X = S, FUN = FUN, ...) 6: doTryCatch(return(expr), name, parentenv, handler) 7: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 8: tryCatchList(expr, classes, parentenv, handlers) 9: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcal\ l <- deparse(call)[1L] prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nc\ har(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b")

  • nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefi\ x, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))\ ) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 10: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 11: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 12: FUN(X[[i]], ...) 13: lapply(seq_len(cores), inner.do) 14: mclapply(1:length(sampleRanges), mc.cores = nCores, FUN = loadBamAndConvert_across_a_range, sampleRanges = sampleRanges, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nS\ NPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov,\ sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, chrLen\ gth = chrLength, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUpperLimit) 15: generate_input(bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, reference = reference, iSizeUpperLimit = iSizeUpperLi\ mit, bqFilter = bqFilter, chr = chr, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inputdir = inputdir, useSoftClippedBases = useSoftClippedBases, regionName\ = regionName, tempdir = tempdir, chrStart = chrStart, chrEnd = chrEnd, nCores = nCores, save_sampleReadsInfo = save_sampleReadsInfo, use_bx_tag = use_bx_tag, bxTagUpperLimit = bxTagUp\ perLimit) 16: generate_or_refactor_input(regenerateInput = regenerateInput, bundling_info = bundling_info, L = L, pos = pos, nSNPs = nSNPs, bam_files = bam_files, cram_files = cram_files, refere\ nce = reference, iSizeUpperLimit = iSizeUpperLimit, bqFilter = bqFilter, chr = chr, outputdir = outputdir, N = N, downsampleToCov = downsampleToCov, sampleNames = sampleNames, inpu\ 17: STITCH(chr = as.character(chr), posfile = posfile, K = 19, S = 1, nGen = 7, outputdir = outputdir, bamlist = "bamfiles.complete.txt", method = "diploid", nCores = 30, splitReadIter\ ations = NA, tempdir = "/scratch0/ucbtrfm/stitch/", regionStart = start, regionEnd = end, buffer = buffer) 18: eval(ei, envir) 19: eval(ei, envir) 20: withVisible(eval(ei, envir)) 21: source("stitch.complete.R")

Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: [2023-09-13 19:33:08] downsample sample A83215_1 - 521 of 283482 reads removed [2023-09-13 19:33:12] downsample sample A83196_1 - 154 of 263333 reads removed [2023-09-13 19:33:14] downsample sample A83216_1 - 19 of 258735 reads removed [2023-09-13 19:33:16] downsample sample A83256_1 - 2 of 196172 reads removed

From: rwdavies @.> Date: Wednesday, 13 September 2023 at 18:32 To: rwdavies/STITCH @.> Cc: Mott, Richard @.>, Author @.> Subject: Re: [rwdavies/STITCH] stitch crashed (Issue #86)

⚠ Caution: External sender

What version is this? Is this a recent one?

On Wed, 13 Sep 2023 at 10:21, Richard Mott @.***> wrote:

And here is the script:

library(STITCH) options(bitmapType='cairo') chunksize = 5000000 buffer = 5000 outputdir="./tair10.complete.STITCH/" chrlen = read.delim("./tair10.chrlen.txt") rownames(chrlen) = chrlen$chr for( chr in 1:5 ) { len = chrlen$len[chr] chunks = floor(len/chunksize)+1 cat("chr", chr, "chunks", chunks, "\n") start=1 end = 0 posfile=paste0("chr",chr,".pos.txt") for( k in 1:chunks) { start = end+1 end = start + chunksize-1 if ( end > len ) end = len cat("chr", chr, "chunk", k, "start", start, "end", end, "\n") STITCH( chr=as.character(chr), posfile=posfile, K=19,S=1,nGen=7,outputdir=outputdir,

bamlist="bamfiles.complete.txt",method="diploid",nCores=30,splitReadIterations

= NA,tempdir=" /scratch0/ucbtrfm/stitch/", regionStart=start, regionEnd=end, buffer=buffer) } }

— Reply to this email directly, view it on GitHub https://github.com/rwdavies/STITCH/issues/86#issuecomment-1717732403, or unsubscribe < https://github.com/notifications/unsubscribe-auth/AB7I36PXHWBNPWXPV3ZNHQLX2G6OLANCNFSM6AAAAAA4WPZGMY>

. You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub< https://github.com/rwdavies/STITCH/issues/86#issuecomment-1718041720>, or unsubscribe< https://github.com/notifications/unsubscribe-auth/AADK6MVIA5TXLPF6ZKJPOGLX2HU4BANCNFSM6AAAAAA4WPZGMY>.

You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/rwdavies/STITCH/issues/86#issuecomment-1718131157, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7I36MEDWCFD3STSMYYXNDX2H4MHANCNFSM6AAAAAA4WPZGMY . You are receiving this because you commented.Message ID: @.***>