hrd_score/calculate_HRD.R

args <- commandArgs(TRUE)

if(length(args)<1){
  args <- c("--help")
}

##help section
if("--help"%in%args){
  cat("Arguments:
      --out_dir=out_directory             The output directory
      --Tumor_LogR_file=tumor_LogR        The output tumor LogR file from ascatNgs
      --Tumor_BAF_file=tumor_BAF          The output tumor BAF file from ascatNgs
      --Germline_LogR_file=Normal_LogR    The output normal LogR file from ascatNgs
      --Germline_BAF_file=Normal_BAF      The output normal BAF file from ascatNgs
      --HRD_file_prefix=prefix            The name prefix of HRD score results

      These parameters must be set in orders!
      ")
  q(save = "no")
}

message(paste0("\n\nYou have the Rights to copy, edit, send and distribute this files,",
               " but not as Commercial activities.\nEnjoy!\n"))

args[[1]] <- gsub("--out_dir=","",args[[1]])
args[[2]] <- gsub("--Tumor_LogR_file=","",args[[2]])
args[[3]] <- gsub("--Tumor_BAF_file=","",args[[3]])
args[[4]] <- gsub("--Germline_LogR_file=","",args[[4]])
args[[5]] <- gsub("--Germline_BAF_file=","",args[[5]])
args[[6]] <- gsub("--HRD_file_prefix=","",args[[6]])





#################################################
#################################################
########     Deal with ASCAT results     ########
########        Output of ascatNgs       ########
#################################################
#################################################
library(ASCAT)
source("f:/choppy_app/HRD_score/HRD_functions.R")

#setwd(args[[1]])
rep_1_ascat.bc <- ascat.loadData(Tumor_LogR_file = args[[2]],
                               Tumor_BAF_file = args[[3]],
                               Germline_LogR_file = args[[4]],
                               Germline_BAF_file = args[[5]],
                               gender = "XX",
                               sexchromosomes = c("X","Y"))
#dir.create("../res/")
setwd(args[[1]])
##no GC content correction here
#rep_1_ascat.bc <- ascat.GCcorrect(rep_1_ascat.bc,GCcontentfile = "../SnpGcCorrections.tsv")
ascat.plotRawData(rep_1_ascat.bc)
rep_1_ascat.bc = ascat.aspcf(rep_1_ascat.bc)
ascat.plotSegmentedData(rep_1_ascat.bc)
rep_1_ascat.output = ascat.runAscat(rep_1_ascat.bc)


#################################################
#################################################
########    Calculate Three HRD score    ########
#################################################
#################################################

rep_1_x <- cbind(rep_1_ascat.output$segments[,1:4],
                 number_of_propes=count_total_numbers(segment_data = rep_1_ascat.output$segments,
                                                      SNP_pos_info = rep_1_ascat.bc$SNPpos)[,1],
                 total_copy_number=apply(rep_1_ascat.output$segments[,5:6],1,sum),
                 rep_1_ascat.output$segments[,5:6],
                 ploidy=rep_1_ascat.output$ploidy,contamination=rep_1_ascat.output$aberrantcellfraction,
                 stringsAsFactors = FALSE)


##load chromsome information
load("f:/choppy_app/HRD_score/chrominfo.RData")			 
				
##modify names of chromosomes to match the chromosomes in the chrominfo object
rep_1_x[rep_1_x[,2]=="X",2] <- 23
rep_1_x[rep_1_x[,2]=="Y",2] <- 24

rep_1_ntai <- calc.ai(rep_1_x, chrominfo,check.names = FALSE,return.loc = FALSE)
rep_1_lst <- calc.lst(rep_1_x,chrominfo,check.names = FALSE,return.loc = FALSE)
rep_1_hrd <- calc.hrd(rep_1_x,check.names = FALSE,return.loc = FALSE)

HRD_score <- rep_1_ntai[1,1]+rep_1_lst+rep_1_hrd
names(HRD_score) <- "HRD_score"

write.table(HRD_score,file = paste0(args[[1]],"/",args[[6]],".txt"),col.names = FALSE,quote = FALSE)

q(save="no")