inst/script/listvargeno_transform.R

In Halvee/rvatk: Rare Variant Analysis Toolkit

#!/usr/bin/Rscript

require(optparse)
require(data.table)

require(rvatk)

# GLOBAL VARS
kOutputTypes <- c("qualvar","long","collapsing.matrix")

Main <- function() {
  option.list <- GetOptList()
  opt <- parse_args(OptionParser(option_list=option.list))
  geno <- ReadLargeTable(opt$genotypes.csv, 
                         showProgress=F,
                         header=T, check.names=F)

  # read in sample file
  sampleped <- read.table(opt$sampleped.file, stringsAsFactors=F)
  cases <- sampleped[sampleped[,6]==2, 2]
  ctrls <- sampleped[sampleped[,6]==1, 2]
  samples <- c(cases, ctrls)

  # verify no repeats of sampleID
  if (TRUE %in% duplicated(samples)) {
    stop("duplicated samples in input sample file")
  } 

  # read in qualvar file if defined, else read in genotypes CSV
  geno[[opt$gene.name.col]] <- gsub("'","", geno[[opt$gene.name.col]]) 
  cnd <- ReadConditionFile(opt$condition.file)
  if (nrow(cnd) > 0) {
    if ("loo maf" %in% cnd[,1]) {
      cat("Adding loo maf col...\n")
      geno <- AddLooMafCol(geno, length(samples))
    }
    geno <- TblCndSubset(geno, cnd)
  }

  # make sure that only sample from sampleped are included in geno
  geno <- geno[ geno[[opt$sample.name.col]] %in% samples, ]
  
  allele.count <- GetAlleleCounts(geno[[opt$sample.name.col]],
                                  geno[[opt$variant.id.col]],
                                  geno[[opt$genotype.col]],
                                  sampleped)

  if (opt$output.type == "collapsing.matrix") {
    if (is.null(opt$geneset.file)) {
      cat("to create collapsing matrix output need geneset.\n")
      q()
    }
    # read in geneset file
    geneset <- ReadListFile(opt$geneset.file)

    # verify no duplicated genes in geneset
    if (TRUE %in% duplicated(geneset)) {
      stop("duplicated gene names in input geneset file")
    }

    # build collapsing matrix based on geneset and sample files
    mat <- CollapsingMatrixInit(geneset, samples)
    qualvar <- cbind(geno[, c(opt$gene.name.col, opt$variant.id.col,
                              opt$sample.name.col), drop=F], allele.count)
    mat <- CollapsingMatrixRecode(mat,        
                                  qualvar)                                    
    WriteLargeTable(mat, opt$output.file,
                    col.names=T,
                    row.names=F, 
                    quote=F, sep="\t")
  } else if (opt$output.type == "long") {
    output <- geno[, c(opt$gene.name.col, opt$variant.id.col,
                       opt$sample.name.col), drop=F]
  } else {
    qualvar <- cbind(geno[, c(opt$gene.name.col, opt$variant.id.col,
                              opt$sample.name.col), drop=F], allele.count)
    if (is.null(opt$geneset.file)==F) {
      geneset <- ReadListFile(opt$geneset.file)
      qualvar <- qualvar[qualvar[,1] %in% geneset, ]
    }
    WriteLargeTable(qualvar, opt$output.file,
                    col.names=F,
                    row.names=F, quote=F, sep="\t")
  }
}

GetOptList <- function() {
  # Get option list.
  #
  # Returns:
  #   Option list to parse options.
  option.list <- list(
    make_option(c("--genotypes-csv"), action="store", type="character",
                default=NULL, dest="genotypes.csv",
                help="name of genotypes csv generated by ATAV,
                      if no qualvar file provided, required, [default %default]"),

    make_option(c("--output-file"),
                action="store", type="character",
                default=NULL,
                dest="output.file",
                help="output file name, required, [default %default]"),

    make_option(c("--output-type"),
                action="store", type="character",
                default="collapsing.matrix",
                dest="output.type",
                help="format for output file, [default %default]"),

    make_option(c("--sampleped-file"), action="store", type="character",
                default=NULL, dest="sampleped.file", 
                help="name of required sample file, [default %default]"),

    make_option(c("--geneset-file"), action="store", 
                type="character",      
                default=NULL, 
                dest="geneset.file",                          
                help="name of required geneset file, [default %default]"), 

    make_option(c("--condition-file"), action="store", type="character",
                default=NULL, dest="condition.file",
                help="name of condition file to subset genotype csv on,
                      if no qualvar file is provided [default %default]"),

    make_option(c("--gene-name-col"), action="store", type="character",
                default="Gene Name", dest="gene.name.col",
                help="Gene name column to subset genotype csv on,
                      if no qualvar file is provided [default %default]"),

    make_option(c("--variant-id-col"), action="store", type="character",
                default="Variant ID", dest="variant.id.col",
                help="variant ID column to subset genotype csv on,
                      if no qualvar file is provided [default %default]"),

    make_option(c("--sample-name-col"), action="store", type="character",
                default="Sample Name", dest="sample.name.col",
                help="Sample name column to subset genotype csv on,
                      if no qualvar file is provided [default %default]"),

    make_option(c("--genotype-col"), action="store", type="character",
                default="Genotype", dest="genotype.col",
                help="Sample name column to subset genotype csv on,
                      if no qualvar file is provided [default %default]")

  )

  return(option.list)
}

if (interactive() == F) {
  Main()
}