kashiff007
commented
3 years ago I solve the this problem by just assigning the genome to reference_file variable:
reference_file <- ("Genome_CC7_new/CC7_R6_locked.softmasked.scaffolds.fa")
But while running:
snp.list <- extract_bams(bam_files, vcf_files, sample_names, cores=10, reference_file, coverage = 2)
I got the following error:
Reading reference file
Running sample CC7_25_1
Reading VCF file
Scanning file to determine attributes.
File attributes:
meta lines: 289
header_line: 290
variant count: 36668556
column count: 10
Meta line 289 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
Character matrix gt rows: 36668556
Character matrix gt cols: 10
skip: 0
nrows: 36668556
row_num: 0
Processed variant 30978000
*** caught segfault ***
address 0x110, cause 'memory not mapped'
Traceback:
1: .read_body_gz(file, stats = stats, nrows = nrows, skip = skip, cols = cols, convertNA = as.numeric(convertNA), verbose = as.numeric(verbose))
2: vcfR::read.vcfR(vcfFiles[samp], verbose = verbose)
3: vcfR::getFIX(vcfR::read.vcfR(vcfFiles[samp], verbose = verbose))
4: FUN(X[[i]], ...)
5: lapply(sample_list, function(samp) { message(sprintf("Running sample %s", sampleNames[samp])) if (verbose) message("Reading VCF file", appendLF = TRUE) vcf <- vcfR::getFIX(vcfR::read.vcfR(vcfFiles[samp], verbose = verbose)) bam.file <- bamFiles[samp] if (verbose) message("Extracting methylation per SNP", appendLF = TRUE) snp.table <- parallel::mcmapply(function(t, u, v) { snp <- GRanges(seqnames = t, IRanges(start = as.integer(u), width = 1)) if (length(grep(t, c(as.character(1:22), "X", "Y"))) == 0) { if (verbose) message("Bad chrom", appendLF = TRUE) return(NULL) } chrom <- t params <- Rsamtools::ScanBamParam(tag = c("MD", "XM", "XR", "XG"), which = snp) alns.pairs <- readGAlignmentPairs(bam.file, use.names = TRUE, param = params) alns <- unlist(alns.pairs) split <- splitReads(alns, v, snp) alt.reads <- split$alt.reads ref.reads <- split$ref.reads if (length(ref.reads) <= 1 | length(alt.reads) <= 1) { if (verbose) message("0 or 1 read in one or both alleles", appendLF = TRUE) return(NULL) } left <- min(start(alns)) right <- max(end(alns)) window <- GRanges(GenomeInfoDb::seqnames(snp), IRanges(left, right)) dna <- Rsamtools::scanFa(fa, param = window) cgsite <- stringr::str_locate_all(dna, "CG")[[1]][, 1] if (length(cgsite) < 1) { if (verbose) message("No CpG sites associated to this SNP", appendLF = TRUE) return(NULL) } mepos <- cgsite/Biostrings::nchar(dna) conversion <- vapply(alns.pairs, function(x) { XGcondition <- mcols(x)$XG[1] if (XGcondition == "GA") { cgsite <- cgsite + 1 } conversion <- matrix(NA, 1, length(cgsite)) read.start <- start(x) - left + 1 read.end <- end(x) - left + 1 for (pair in 1:2) { something <- mcols(x[pair])$MD tag <- getMD(something) MDtag <- tag$MDtag nucl.num <- tag$nucl.num count <- 0 this.mepos <- cgsite - read.start[pair] + 1 for (p in seq_along(this.mepos)) { if (is.na(conversion[1, p])) { if (cgsite[p] > read.end[pair] || cgsite[p] < read.start[pair]) { conversion[1, p] <- NA next } else { for (i in seq_along(nucl.num)) { count <- count + nucl.num[i] if (count >= this.mepos[p]) { count <- 0 break } } if (MDtag[i] %in% c("C", "G")) { conversion[1, p] <- 1 } else { conversion[1, p] <- 2 } } } else { next } } } return(conversion) }, double(length(cgsite))) if (is.null(dim(conversion))) { ref.conv <- conversion[names(conversion) %in% ref.reads] alt.conv <- conversion[names(conversion) %in% alt.reads] ref.cov <- sum(!is.na(ref.conv)) alt.cov <- sum(!is.na(alt.conv)) ref.meth <- sum(!is.na(ref.conv) & ref.conv == 2) alt.meth <- sum(!is.na(alt.conv) & alt.conv == 2) } else { ref.conv <- conversion[, colnames(conversion) %in% ref.reads] alt.conv <- conversion[, colnames(conversion) %in% alt.reads] ref.cov <- BiocGenerics::rowSums(!is.na(ref.conv)) alt.cov <- BiocGenerics::rowSums(!is.na(alt.conv)) ref.meth <- BiocGenerics::rowSums(!is.na(ref.conv) & ref.conv == 2) alt.meth <- BiocGenerics::rowSums(!is.na(alt.conv) & alt.conv == 2) } GR <- GRanges(gsub("chr", "", chrom), IRanges(start = left + cgsite, width = 1)) mcols(GR)$cov.ref <- ref.cov mcols(GR)$cov.alt <- alt.cov mcols(GR)$meth.ref <- ref.meth mcols(GR)$meth.alt <- alt.meth mcols(GR)$snp <- paste0(GenomeInfoDb::seqnames(snp), ".", start(snp)) filt <- BiocGenerics::rowSums(cbind(GR$cov.ref, GR$cov.alt) >= coverage) >= 2 if (sum(filt) < 2) { if (verbose) message("No CpG sites sufficiently covered at this SNP", appendLF = TRUE) return(NULL) } gr <- GR[filt] return(gr) }, t = vcf[, 1], u = vcf[, 2], v = vcf[, 4], SIMPLIFY = F, USE.NAMES = F, mc.cores = cores) message(sprintf("Done with sample %s", sampleNames[samp])) return(snp.table)})
6: extract_bams(bam_files, vcf_files, sample_names, reference_file, coverage = 2)
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:
What is the possible mistake I am doing?
Thanks in advance, Kashif
sorjuela
Hi,
I know from vignette how to upload one chr/scaffold into reference_file
reference_file <- DNAStringSet(genome[[19]], use.names = TRUE)But having problem in uploading whole genome fasta file containing multiple chr/scaffolds. My genome looks like:
I also looked into manual of extract_bams, but couldn't found the solution.
Can you provide a solution/command by which I can perform this operation.
Thank you, Kashif