inst/example/battenberg_snp6_refitting.R

In Wedge-lab/battenberg: Battenberg subclonal copy number caller

args = commandArgs(TRUE)
TUMOURNAME = toString(args[1])
RUN_DIR = toString(args[2])
PRESET_RHO = as.numeric(args[3])
PRESET_PSI = as.numeric(args[4])

library(Battenberg)

###############################################################################
# 2015-05-01
# A pure R Battenberg v2.0.0 SNP6 refitting pipeline implementation.
# sd11@sanger.ac.uk
###############################################################################

# Sample specific
#IS.MALE = F
# TUMOURNAME = "NASCR-0016"
# NORMALCEL = "/nfs/cgpstats1/pvl/ASCAT/NeoAva/CELfiles/NASCR-0016B1.CEL"
# TUMOURCEL = "/nfs/cgpstats1/pvl/ASCAT/NeoAva/CELfiles/NASCR-0016.CEL"
# RUN_DIR = "/lustre/scratch110/sanger/sd11/battenberg_package_test/NASCR-0016_bb_v2.0_singlecore"

# Parallelism parameters
NTHREADS = 6

# General static
IMPUTEINFOFILE = "/lustre/scratch110/sanger/sd11/Documents/GenomeFiles/battenberg_impute/impute_info.txt"
G1000PREFIX = "/lustre/scratch110/sanger/sd11/Documents/GenomeFiles/battenberg_1000genomesloci2012/1000genomesAlleles2012_chr"
IMPUTE_EXE = "impute2"

# General SNP6 specific
PROBLEMLOCI = NA
SNP6_REF_INFO_FILE = "/lustre/scratch110/sanger/sd11/Documents/GenomeFiles/battenberg_snp6/snp6_ref_info_file.txt"
APT_PROBESET_GENOTYPE_EXE = "apt-probeset-genotype"
APT_PROBESET_SUMMARIZE_EXE = "apt-probeset-summarize"
NORM_GENO_CLUST_EXE = "~pvl/PennCNV/gw6/bin/normalize_affy_geno_cluster.pl"
BIRDSEED_REPORT_FILE = "birdseed.report.txt" # No control over the name of this file, as it is automatically generated by APT within cel2baf.logr

# Parameters
PLATFORM_GAMMA = 0.55
PHASING_GAMMA = 1
SEGMENTATION_GAMMA = 10
CLONALITY_DIST_METRIC = 0
ASCAT_DIST_METRIC = 1
MIN_PLOIDY = 1.6
MAX_PLOIDY = 4.8
MIN_RHO = 0.1
MIN_GOODNESS_OF_FIT = 0.63
BALANCED_THRESHOLD = 0.51
MIN_NORMAL_DEPTH = 10

# Change to work directory and load the chromosome information
setwd(RUN_DIR)
chrom_names = get.chrom.names(IMPUTEINFOFILE, TRUE)

# Infer what the gender is
gender = infer_gender_birdseed(BIRDSEED_REPORT_FILE)
is_male = gender == "male"
chrom_names = get.chrom.names(IMPUTEINFOFILE, is_male)

# Fit a clonal copy number profile
fit.copy.number(samplename=TUMOURNAME,
                outputfile.prefix=paste(TUMOURNAME, "_", sep=""),
                inputfile.baf.segmented=paste(TUMOURNAME, ".BAFsegmented.txt", sep=""), 
                inputfile.baf=paste(TUMOURNAME,"_mutantBAF.tab", sep=""), 
                inputfile.logr=paste(TUMOURNAME,"_mutantLogR.tab", sep=""), 
                dist_choice=CLONALITY_DIST_METRIC, 
                ascat_dist_choice=ASCAT_DIST_METRIC, 
                min.ploidy=MIN_PLOIDY, 
                max.ploidy=MAX_PLOIDY, 
                min.rho=MIN_RHO, 
                min.goodness=MIN_GOODNESS_OF_FIT, 
                uninformative_BAF_threshold=BALANCED_THRESHOLD, 
                gamma_param=PLATFORM_GAMMA, 
                use_preset_rho_psi=T, 
                preset_rho=PRESET_RHO, 
                preset_psi=PRESET_PSI, 
                read_depth=30)

# Go over all segments, determine which segements are a mixture of two states and fit a second CN state
callSubclones(sample.name=TUMOURNAME, 
              baf.segmented.file=paste(TUMOURNAME, ".BAFsegmented.txt", sep=""), 
              logr.file=paste(TUMOURNAME,"_mutantLogR.tab", sep=""), 
              rho.psi.file=paste(TUMOURNAME, "_rho_and_psi.txt",sep=""), 
              output.file=paste(TUMOURNAME,"_subclones.txt", sep=""), 
              output.figures.prefix=paste(TUMOURNAME,"_subclones_chr", sep=""), 
              output.gw.figures.prefix=paste(TUMOURNAME,"_BattenbergProfile", sep=""),
              chr_names=chrom_names, 
              gamma=PLATFORM_GAMMA, 
              segmentation.gamma=NA, 
              siglevel=0.05, 
              maxdist=0.01, 
              noperms=1000)