From bafad0d186a338a419d3eed3d043ac9ef796a591 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Fri, 30 Sep 2022 09:11:33 -0400 Subject: [PATCH 01/11] update --- R/HRD_wrapper.R | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 R/HRD_wrapper.R diff --git a/R/HRD_wrapper.R b/R/HRD_wrapper.R new file mode 100644 index 0000000..9d02ac5 --- /dev/null +++ b/R/HRD_wrapper.R @@ -0,0 +1,30 @@ +#!/usr/bin/env Rscript +library(HRDex) +library(argparse) + +args = commandArgs(trailingOnly = TRUE) + +if( length(args) < 2) +{ + stop("need to provide 1 argument: seq data file, reference genome (grch37 or grch38)") +} + +seq.dat <- read.table(args[1], header = TRUE) +ref <- args[2] +sample.id <- args[3] + +#column names for facets +#csv text or csv + +seq.dat <- preprocessHRD( seq.dat, ref )#needs required fields +#Doies not check for A or B NA values, which facet may have. +CN.dat <- getCNt(seq.dat ) +HRD <- getHRD.Score( seq.dat, CN.dat, scaleTotal = FALSE) + +#HRD.NTAIr <- getNTAI.raw( seq.dat ) +HRD.LST <- getLST( seq.dat ) +HRD.LOH <- getLOH( seq.dat ) +HRD.NTAIm <- getNTAI.norm( seq.dat, CN.dat ) + + +cat(paste(sample.id,HRD.LOH,HRD.NTAIm,HRD.LST,HRD,sep=",")) From 4fcf374d70ab6ed5d941ff09344abdb252a5fbae Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Thu, 16 Feb 2023 12:28:23 -0500 Subject: [PATCH 02/11] Update rename --- R/HRD_wrapper.R | 30 ------------------------------ R/run_HRDex.R | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 30 deletions(-) delete mode 100644 R/HRD_wrapper.R create mode 100644 R/run_HRDex.R diff --git a/R/HRD_wrapper.R b/R/HRD_wrapper.R deleted file mode 100644 index 9d02ac5..0000000 --- a/R/HRD_wrapper.R +++ /dev/null @@ -1,30 +0,0 @@ -#!/usr/bin/env Rscript -library(HRDex) -library(argparse) - -args = commandArgs(trailingOnly = TRUE) - -if( length(args) < 2) -{ - stop("need to provide 1 argument: seq data file, reference genome (grch37 or grch38)") -} - -seq.dat <- read.table(args[1], header = TRUE) -ref <- args[2] -sample.id <- args[3] - -#column names for facets -#csv text or csv - -seq.dat <- preprocessHRD( seq.dat, ref )#needs required fields -#Doies not check for A or B NA values, which facet may have. -CN.dat <- getCNt(seq.dat ) -HRD <- getHRD.Score( seq.dat, CN.dat, scaleTotal = FALSE) - -#HRD.NTAIr <- getNTAI.raw( seq.dat ) -HRD.LST <- getLST( seq.dat ) -HRD.LOH <- getLOH( seq.dat ) -HRD.NTAIm <- getNTAI.norm( seq.dat, CN.dat ) - - -cat(paste(sample.id,HRD.LOH,HRD.NTAIm,HRD.LST,HRD,sep=",")) diff --git a/R/run_HRDex.R b/R/run_HRDex.R new file mode 100644 index 0000000..8745e9d --- /dev/null +++ b/R/run_HRDex.R @@ -0,0 +1,33 @@ +#!/usr/bin/env Rscript +library(HRDex) +library(argparse) + +p<-ArgumentParser() +p$add_argument('-i','--infile',help='Path to input file.') +p$add_argument('-o','--outfile',help="Path to output file.") +p$add_argument('--tumor',help='Tumor id.') +p$add_argument('--build',help='Genome build; {grch37,grch38}.') + +args<-p$parse_args() + +seq.dat <- read.table(args$infile, header = TRUE) +ref <- tolower(args$build) +sample.id <- args$tumor + +#column names for facets +#csv text or csv + +seq.dat <- preprocessHRD( seq.dat, ref )#needs required fields +#Does not check for A or B NA values, which facet may have. + +CN.dat <- getCNt(seq.dat ) + +#outputting all values for now. +HRD.score <- getHRD.Score( seq.dat, CN.dat, scaleTotal = FALSE) +HRD.LST <- getLST( seq.dat ) +HRD.LOH <- getLOH( seq.dat ) +HRD.NTAIr <- getNTAI.raw( seq.dat ) +HRD.NTAIm <- getNTAI.norm( seq.dat, CN.dat ) + +output=data.frame("TumorID"=args$tumor,"HRD.NTAIr"=HRD.NTAIr,"HRD.NTAIm"=HRD.NTAIm,"HRD.LOH"=HRD.LOH,"HRD.LST"=HRD.LST,"HRD.score"=HRD.score) +write.csv(output,file=args$outfile,row.names=FALSE,col.names=TRUE) From 29117e3761c68e9603907d85f5461cd0ca9f6793 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Tue, 28 Feb 2023 20:56:47 -0500 Subject: [PATCH 03/11] Update increase NA checking. Chromosome name changes --- R/preprocessHRD.R | 53 ++++++++++++++++++++++++++++------------------- R/run_HRDex.R | 2 +- 2 files changed, 33 insertions(+), 22 deletions(-) diff --git a/R/preprocessHRD.R b/R/preprocessHRD.R index 93ccd2a..58e95e4 100644 --- a/R/preprocessHRD.R +++ b/R/preprocessHRD.R @@ -1,12 +1,12 @@ #' Define additional properties of the sequencing data necessary to compute HRD -#' +#' #' @param seq.dat the raw sequencing data #' @param ref the reference genome to use (grch37, grch38) -#' @details define chromosome size, arms, telomere position, and centromere position from reference data; +#' @details define chromosome size, arms, telomere position, and centromere position from reference data; #' compute allelic imbalance and segment size -#' +#' #' @return seq.dat -#' +#' #'@examples #'seq.dat <- sub01.segments[ sub01.segments$chromosome == "chr1",] #' @@ -36,59 +36,70 @@ preprocessHRD <- function( seq.dat, ref ) print(paste(ref, "is not a valid reference genome.", sep = " ")) stop("select one of: grch37, grch38") } - + # check the sequencing data input to make sure it has the data we need if( any(!(seq.cols.needed %in% colnames(seq.dat)))) { - print(paste("column", + print(paste("column", seq.cols.needed[which(!(seq.cols.needed %in% colnames(seq.dat)))], "missing in the seq data.", sep=" ")) print(paste("Looking for columns:", seq.cols.needed, sep=" ")) print(paste("Found:", seq.cols.needed[seq.cols.needed %in% colnames(seq.dat)])) stop("Column mismatch. Exiting.") - } - + } + + # fix chromosome numbers + seq.dat$chromosome[seq.dat$chromosome==23,]<-"X" + seq.dat<-seq.dat[!seq.dat$chromosome %in% c(24,"Y"),] + # remove rows w/ missing entries if( any( is.na(seq.dat$CNt) )) { print(paste("Row: ", which(is.na(seq.dat$CNt)), " contains NA; removing.", sep = "")) seq.dat <- subset(seq.dat, !is.na(seq.dat$CNt)) } - + + if( any( is.na(seq.dat$A) | is.na(seq.dat$B))) + { + seq.dat <- seq.dat[-which(is.na(seq.dat$A) | is.na(seq.dat$B)),] + } + levels(seq.dat$chromosome) <- levels(ref.dat$chromosome) - + seq.dat$frac.chr <- (seq.dat$end.pos - seq.dat$start.pos) / ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] - + + + # segment length seq.dat$brk.len <- 0 seq.dat <- combineSeg(seq.dat, 3e06) - + # matches reference data to the correct chromosome in the subject data key = match(seq.dat$chromosome, ref.dat$chromosome) - + seq.dat$seg.len <- seq.dat$end.pos - seq.dat$start.pos seq.dat$chr.size <- ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] - + seq.dat$AI <- (seq.dat$A > seq.dat$B) & (seq.dat$A != 1) & (seq.dat$B != 0) - + seq.dat$start.arm <- seq.dat$end.arm <- rep("NA", dim(seq.dat)[1]) - + ind = seq.dat$start.pos - ref.dat$centromere.start[key] > 1 seq.dat$start.arm[ind] = "p" seq.dat$start.arm[!ind] = "q" - + ind = seq.dat$end.pos - ref.dat$centromere.end[key] > 1 seq.dat$end.arm[ind] = "p" seq.dat$end.arm[!ind] = "q" - + seq.dat$cross.arm <- seq.dat$start.arm != seq.dat$end.arm - - + + # define end telomeres seq.dat$post.telomere <- (seq.dat$start.pos - ref.dat$p.telomere.end[key] <= 1000) seq.dat$pre.telomere <- (ref.dat$q.telomere.start[key] - seq.dat$end.pos <= 1000) # --- - + return(seq.dat) } # ------------------------------------------------------------------------------------- # diff --git a/R/run_HRDex.R b/R/run_HRDex.R index 8745e9d..9209e19 100644 --- a/R/run_HRDex.R +++ b/R/run_HRDex.R @@ -30,4 +30,4 @@ HRD.NTAIr <- getNTAI.raw( seq.dat ) HRD.NTAIm <- getNTAI.norm( seq.dat, CN.dat ) output=data.frame("TumorID"=args$tumor,"HRD.NTAIr"=HRD.NTAIr,"HRD.NTAIm"=HRD.NTAIm,"HRD.LOH"=HRD.LOH,"HRD.LST"=HRD.LST,"HRD.score"=HRD.score) -write.csv(output,file=args$outfile,row.names=FALSE,col.names=TRUE) +write.table(output,file=args$outfile,sep=",",row.names=FALSE,col.names=TRUE) From 5b774dbe05a935e9bcd366f78242e019b8e2dc8c Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 1 Mar 2023 12:26:54 -0500 Subject: [PATCH 04/11] Delete run_HRDex.R --- R/run_HRDex.R | 33 --------------------------------- 1 file changed, 33 deletions(-) delete mode 100644 R/run_HRDex.R diff --git a/R/run_HRDex.R b/R/run_HRDex.R deleted file mode 100644 index 9209e19..0000000 --- a/R/run_HRDex.R +++ /dev/null @@ -1,33 +0,0 @@ -#!/usr/bin/env Rscript -library(HRDex) -library(argparse) - -p<-ArgumentParser() -p$add_argument('-i','--infile',help='Path to input file.') -p$add_argument('-o','--outfile',help="Path to output file.") -p$add_argument('--tumor',help='Tumor id.') -p$add_argument('--build',help='Genome build; {grch37,grch38}.') - -args<-p$parse_args() - -seq.dat <- read.table(args$infile, header = TRUE) -ref <- tolower(args$build) -sample.id <- args$tumor - -#column names for facets -#csv text or csv - -seq.dat <- preprocessHRD( seq.dat, ref )#needs required fields -#Does not check for A or B NA values, which facet may have. - -CN.dat <- getCNt(seq.dat ) - -#outputting all values for now. -HRD.score <- getHRD.Score( seq.dat, CN.dat, scaleTotal = FALSE) -HRD.LST <- getLST( seq.dat ) -HRD.LOH <- getLOH( seq.dat ) -HRD.NTAIr <- getNTAI.raw( seq.dat ) -HRD.NTAIm <- getNTAI.norm( seq.dat, CN.dat ) - -output=data.frame("TumorID"=args$tumor,"HRD.NTAIr"=HRD.NTAIr,"HRD.NTAIm"=HRD.NTAIm,"HRD.LOH"=HRD.LOH,"HRD.LST"=HRD.LST,"HRD.score"=HRD.score) -write.table(output,file=args$outfile,sep=",",row.names=FALSE,col.names=TRUE) From d26fe7edcd1ee4ea5807cd2a814a7ef12b41a299 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Mon, 6 Mar 2023 21:44:42 -0500 Subject: [PATCH 05/11] Update preprocessHRD.R fix --- R/preprocessHRD.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/preprocessHRD.R b/R/preprocessHRD.R index 58e95e4..53dda0a 100644 --- a/R/preprocessHRD.R +++ b/R/preprocessHRD.R @@ -49,7 +49,7 @@ preprocessHRD <- function( seq.dat, ref ) } # fix chromosome numbers - seq.dat$chromosome[seq.dat$chromosome==23,]<-"X" + seq.dat$chromosome[seq.dat$chromosome==23]<-"X" seq.dat<-seq.dat[!seq.dat$chromosome %in% c(24,"Y"),] # remove rows w/ missing entries From 60d7b849b22027443a1c0d49cd8cc0bba520f57a Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 14 Feb 2024 13:12:55 -0500 Subject: [PATCH 06/11] Create runHRDex.R --- R/runHRDex.R | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 R/runHRDex.R diff --git a/R/runHRDex.R b/R/runHRDex.R new file mode 100644 index 0000000..096d818 --- /dev/null +++ b/R/runHRDex.R @@ -0,0 +1,43 @@ +#!/usr/bin/env Rscript +library(HRDex) +library(argparse) + +p <- ArgumentParser() +p$add_argument('-i', '--infile', help = 'Path to input BED file.') +p$add_argument('-o', '--outfile', help = "Path to output file.") +p$add_argument('--tumor', help = 'Tumor id.') +p$add_argument('--build', help = 'Genome build; {grch37,grch38}.') + +args <- p$parse_args() + +# Reading the BED file and specifying column names +seq.dat <- read.table(args$infile, header = FALSE, sep = "\t", stringsAsFactors = FALSE) +colnames(seq.dat) <- c("chromosome", "start.pos", "end.pos", "name", "score", "strand") + +# Extracting A and B values from the 'name' column and calculating CNt +seq.dat$A <- as.numeric(sub(".*;A(\\d+);.*", "\\1", seq.dat$name)) +seq.dat$B <- as.numeric(sub(".*;B(\\d+).*", "\\1", seq.dat$name)) +seq.dat$CNt <- seq.dat$A + seq.dat$B + +# Ensuring chromosome format is consistent for analysis +seq.dat$chromosome <- as.character(seq.dat$chromosome) +seq.dat$chromosome <- ifelse(seq.dat$chromosome == "23", "X", ifelse(seq.dat$chromosome == "24", "Y", seq.dat$chromosome)) + +ref <- tolower(args$build) +sample.id <- args$tumor + +# Now seq.dat includes the required fields: +seq.dat <- preprocessHRD(seq.dat, ref) # Make sure preprocessHRD is ready to accept this structure + +CN.dat <- getCNt(seq.dat) + +# Computing HRD scores +HRD.score <- getHRD.Score(seq.dat, CN.dat, scaleTotal = FALSE) +HRD.LST <- getLST(seq.dat) +HRD.LOH <- getLOH(seq.dat) +HRD.NTAIr <- getNTAI.raw(seq.dat) +HRD.NTAIm <- getNTAI.norm(seq.dat, CN.dat) + +# Preparing output +output <- data.frame(TumorID = args$tumor, HRD.NTAIr = HRD.NTAIr, HRD.NTAIm = HRD.NTAIm, HRD.LOH = HRD.LOH, HRD.LST = HRD.LST, HRD.score = HRD.score) +write.csv(output, file = args$outfile, row.names = FALSE, col.names = TRUE) From 1881fa77a8c45a3062753f988150f5c6ca72994b Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 14 Feb 2024 13:23:31 -0500 Subject: [PATCH 07/11] Update runHRDex.R --- R/runHRDex.R | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/R/runHRDex.R b/R/runHRDex.R index 096d818..ad02260 100644 --- a/R/runHRDex.R +++ b/R/runHRDex.R @@ -6,10 +6,18 @@ p <- ArgumentParser() p$add_argument('-i', '--infile', help = 'Path to input BED file.') p$add_argument('-o', '--outfile', help = "Path to output file.") p$add_argument('--tumor', help = 'Tumor id.') -p$add_argument('--build', help = 'Genome build; {grch37,grch38}.') +p$add_argument('--build', help = 'Genome build; {grch37, grch38, hg19, hg38}.') args <- p$parse_args() +# Mapping hg19 to grch37 and hg38 to grch38 +genome_build_map <- list(grch37 = "grch37", grch38 = "grch38", hg19 = "grch37", hg38 = "grch38") +ref <- tolower(genome_build_map[[tolower(args$build)]]) + +if (is.null(ref)) { + stop("Invalid genome build. Please use one of: grch37, grch38, hg19, hg38") +} + # Reading the BED file and specifying column names seq.dat <- read.table(args$infile, header = FALSE, sep = "\t", stringsAsFactors = FALSE) colnames(seq.dat) <- c("chromosome", "start.pos", "end.pos", "name", "score", "strand") From 5e4c1513ea64efcfe7982dd80e63e5ea94c63ca6 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 14 Feb 2024 16:43:25 -0500 Subject: [PATCH 08/11] Delete runHRDex.R --- R/runHRDex.R | 51 --------------------------------------------------- 1 file changed, 51 deletions(-) delete mode 100644 R/runHRDex.R diff --git a/R/runHRDex.R b/R/runHRDex.R deleted file mode 100644 index ad02260..0000000 --- a/R/runHRDex.R +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/env Rscript -library(HRDex) -library(argparse) - -p <- ArgumentParser() -p$add_argument('-i', '--infile', help = 'Path to input BED file.') -p$add_argument('-o', '--outfile', help = "Path to output file.") -p$add_argument('--tumor', help = 'Tumor id.') -p$add_argument('--build', help = 'Genome build; {grch37, grch38, hg19, hg38}.') - -args <- p$parse_args() - -# Mapping hg19 to grch37 and hg38 to grch38 -genome_build_map <- list(grch37 = "grch37", grch38 = "grch38", hg19 = "grch37", hg38 = "grch38") -ref <- tolower(genome_build_map[[tolower(args$build)]]) - -if (is.null(ref)) { - stop("Invalid genome build. Please use one of: grch37, grch38, hg19, hg38") -} - -# Reading the BED file and specifying column names -seq.dat <- read.table(args$infile, header = FALSE, sep = "\t", stringsAsFactors = FALSE) -colnames(seq.dat) <- c("chromosome", "start.pos", "end.pos", "name", "score", "strand") - -# Extracting A and B values from the 'name' column and calculating CNt -seq.dat$A <- as.numeric(sub(".*;A(\\d+);.*", "\\1", seq.dat$name)) -seq.dat$B <- as.numeric(sub(".*;B(\\d+).*", "\\1", seq.dat$name)) -seq.dat$CNt <- seq.dat$A + seq.dat$B - -# Ensuring chromosome format is consistent for analysis -seq.dat$chromosome <- as.character(seq.dat$chromosome) -seq.dat$chromosome <- ifelse(seq.dat$chromosome == "23", "X", ifelse(seq.dat$chromosome == "24", "Y", seq.dat$chromosome)) - -ref <- tolower(args$build) -sample.id <- args$tumor - -# Now seq.dat includes the required fields: -seq.dat <- preprocessHRD(seq.dat, ref) # Make sure preprocessHRD is ready to accept this structure - -CN.dat <- getCNt(seq.dat) - -# Computing HRD scores -HRD.score <- getHRD.Score(seq.dat, CN.dat, scaleTotal = FALSE) -HRD.LST <- getLST(seq.dat) -HRD.LOH <- getLOH(seq.dat) -HRD.NTAIr <- getNTAI.raw(seq.dat) -HRD.NTAIm <- getNTAI.norm(seq.dat, CN.dat) - -# Preparing output -output <- data.frame(TumorID = args$tumor, HRD.NTAIr = HRD.NTAIr, HRD.NTAIm = HRD.NTAIm, HRD.LOH = HRD.LOH, HRD.LST = HRD.LST, HRD.score = HRD.score) -write.csv(output, file = args$outfile, row.names = FALSE, col.names = TRUE) From 530677a8a41632f2d37d716e9282fb5d97ee353a Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 14 Feb 2024 19:11:38 -0500 Subject: [PATCH 09/11] Update getCNT.R --- R/getCNT.R | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/R/getCNT.R b/R/getCNT.R index 7dab4b6..5c67f9e 100644 --- a/R/getCNT.R +++ b/R/getCNT.R @@ -1,17 +1,17 @@ #' get the main copynumber state for each chromosome -#' +#' #' @param seq.dat the data.frame of sequencing data -#' +#' #' @details Get the main copynumber state for each chromosome. Calculate the total proportion of #' the chromosome associated with each unique copynumber state; return the copynumber state that is most #' represented. -#' +#' #' @return a data.frame containing the main copynumber state for each input chromosome -#' -#' @examples +#' +#' @examples #' seq.dat <- preprocessHRD( seq.dat ) #' CN.dat <- getCNt( seq.dat ) -#' +#' #' @export @@ -24,35 +24,44 @@ getCNt <- function( seq.dat ) # CN.dat (data.frame), then number { main.CN <- rep(0, length(unique(seq.dat$chromosome))) - + ct1 <- 0 for( chr in unique(seq.dat$chromosome)) { ct1 <- ct1 + 1 - + # reduce data to one chromosome at a time dat <- seq.dat[seq.dat$chromosome == chr,] - + # remove any NAs or this will crash dat <- dat[!is.na(dat$CNt),] CNt.frac <- rep(0, length(unique(dat$CNt))) ct2 <- 0 - + # get the fraction of the chromosome association with each copynumber state for( CNt in unique(dat$CNt) ) { ct2 <- ct2 + 1 CNt.frac[ct2] <- sum(dat$frac.chr[dat$CNt == CNt]) } - + # record values + associated CNt out <- data.frame(CNt=unique(dat$CNt), CNt.frac=CNt.frac) - + # get copynumber state most prominent in the chromosome - main.CN[ct1] <- out$CNt[which(CNt.frac == max(CNt.frac))] - + if(length(CNt.frac) > 0 && any(!is.na(CNt.frac))) { + max_frac_indices <- which(CNt.frac == max(CNt.frac, na.rm = TRUE)) + if(length(max_frac_indices) > 0) { + main.CN[ct1] <- out$CNt[max_frac_indices[1]] # Use the first if there are multiple maxima + } else { + main.CN[ct1] <- NA # Or handle differently if no maximum fraction is found + } + } else { + main.CN[ct1] <- NA # Or set a default value or handle as appropriate + } + } - + CN.dat <- data.frame(chromosome=unique(seq.dat$chromosome), main.CN=main.CN) return(CN.dat) } From db5a26b4b1365f04352f14903d1131be1113e3b3 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Wed, 14 Feb 2024 19:24:04 -0500 Subject: [PATCH 10/11] Update preprocessHRD.R --- R/preprocessHRD.R | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/R/preprocessHRD.R b/R/preprocessHRD.R index 53dda0a..217db5a 100644 --- a/R/preprocessHRD.R +++ b/R/preprocessHRD.R @@ -66,9 +66,12 @@ preprocessHRD <- function( seq.dat, ref ) levels(seq.dat$chromosome) <- levels(ref.dat$chromosome) - seq.dat$frac.chr <- (seq.dat$end.pos - seq.dat$start.pos) / ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] - - + if(nrow(seq.dat) > 0) { + seq.dat$frac.chr <- (seq.dat$end.pos - seq.dat$start.pos) / ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] + } else { + print("No data available for calculating frac.chr.") + return(seq.dat) + } # segment length seq.dat$brk.len <- 0 From 434dff76fe1835bb825fd3fc6d9fad71a63357f8 Mon Sep 17 00:00:00 2001 From: Brad Wubbenhorst Date: Thu, 15 Feb 2024 10:52:18 -0500 Subject: [PATCH 11/11] Update added something to conditionally remove prepending chr. Hopefully didnt break other stuff in the process. --- R/getCNT.R | 4 ++-- R/preprocessHRD.R | 14 +++----------- 2 files changed, 5 insertions(+), 13 deletions(-) diff --git a/R/getCNT.R b/R/getCNT.R index 5c67f9e..0b6a999 100644 --- a/R/getCNT.R +++ b/R/getCNT.R @@ -54,10 +54,10 @@ getCNt <- function( seq.dat ) if(length(max_frac_indices) > 0) { main.CN[ct1] <- out$CNt[max_frac_indices[1]] # Use the first if there are multiple maxima } else { - main.CN[ct1] <- NA # Or handle differently if no maximum fraction is found + main.CN[ct1] <- NA } } else { - main.CN[ct1] <- NA # Or set a default value or handle as appropriate + main.CN[ct1] <- NA } } diff --git a/R/preprocessHRD.R b/R/preprocessHRD.R index 217db5a..a2d7c5e 100644 --- a/R/preprocessHRD.R +++ b/R/preprocessHRD.R @@ -51,17 +51,9 @@ preprocessHRD <- function( seq.dat, ref ) # fix chromosome numbers seq.dat$chromosome[seq.dat$chromosome==23]<-"X" seq.dat<-seq.dat[!seq.dat$chromosome %in% c(24,"Y"),] - - # remove rows w/ missing entries - if( any( is.na(seq.dat$CNt) )) - { - print(paste("Row: ", which(is.na(seq.dat$CNt)), " contains NA; removing.", sep = "")) - seq.dat <- subset(seq.dat, !is.na(seq.dat$CNt)) - } - - if( any( is.na(seq.dat$A) | is.na(seq.dat$B))) - { - seq.dat <- seq.dat[-which(is.na(seq.dat$A) | is.na(seq.dat$B)),] + # Check if any values in seq.dat$chromosome do NOT start with "chr" + if(any(!grepl("^chr", seq.dat$chromosome))) { + ref.dat$chromosome <- sub("^chr", "", ref.dat$chromosome) } levels(seq.dat$chromosome) <- levels(ref.dat$chromosome)