view sequenza_to_hrdtools_input.R @ 0:b77d7a0a45e8 draft

"planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/snvtocnv commit 10ad3a0ca7cd23ad1e0940844147e1d1b3d069f0"
author artbio
date Sun, 07 Mar 2021 23:19:59 +0000
parents
children
line wrap: on
line source

options(show.error.messages = F, error = function() {
        cat(geterrmessage(), file = stderr()); q("no", 1, F) })

# load packages that are provided in the conda env

library(optparse)
library(tidyverse)

option_list <- list(
  make_option(
    c("-i", "--input"),
    default = NA,
    type = "character",
    help = "Path to Sequenza output segments file"
  ),
  make_option(
    c("-o", "--output"),
    default = NA,
    type = "character",
    help = "output file, to be used as input for HRDetect"
  ),
  make_option(
    c("-s", "--solutions"),
    default = NA,
    type = "character",
    help = "Path to Sequenza list of alternative solutions"
  )
)

opt <- parse_args(OptionParser(option_list = option_list),
                 args = commandArgs(trailingOnly = TRUE))

sequenza_data <- as_tibble(read.delim(opt$input, header = TRUE))
solutions_data <- as_tibble(read.delim(opt$solutions, header = TRUE))


ploidy <- round(solutions_data$ploidy[1])
cellularity <- solutions_data$cellularity[1]

reformatted <- sequenza_data %>%
  select(
    chr = chromosome,
    start = start.pos,
    end = end.pos,
    copynumber = CNt,
    A, B
  ) %>%
  mutate(
    ploidy = ploidy,
    cellularity = cellularity,
    lohtype = case_when(
      copynumber == 0 ~ "HOMD",
      B == 0 & A == ploidy ~ "NLOH",
      B == 0 & A < ploidy & A > 0 ~ "DLOH",
      copynumber > ploidy & A > B ~ "ASCNA",
      copynumber > ploidy & A == B ~ "BCNA",
      TRUE ~ "HET"
    )
  )

message("Preview of output:")
print(reformatted)

reformatted %>%
  write.table(opt$output, quote = F, row.names = F, sep = "\t")

message(sprintf("Output written to %s", opt$output))