extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/scripts/run_ascat.R

Code
Comments
Other
Rev Date Author Line
7271 29 Jun 23 nicklas 1 #!/usr/bin/env Rscript
7271 29 Jun 23 nicklas 2 #
7271 29 Jun 23 nicklas 3 # Wrapper script for running ASCAT
7271 29 Jun 23 nicklas 4 # 
7271 29 Jun 23 nicklas 5
7278 10 Aug 23 nicklas 6 # Parse arguments
7278 10 Aug 23 nicklas 7 library("argparser")
7278 10 Aug 23 nicklas 8 parser <- arg_parser("Wrapper script for runnint ASCAT", hide.opts=T)
7278 10 Aug 23 nicklas 9 parser <- add_argument(parser, "inputDir", help="Path to directory with LogR and BAF files")
7278 10 Aug 23 nicklas 10 parser <- add_argument(parser, "gcFile", help="Path to GC correction file")
7278 10 Aug 23 nicklas 11 parser <- add_argument(parser, "--outputDir", default=".", help="Path to directory where output files are saved")
7278 10 Aug 23 nicklas 12 parser <- add_argument(parser, "--gender", default="XX", help="Gender: XX or XY")
7278 10 Aug 23 nicklas 13 parser <- add_argument(parser, "--genomeVersion", default="hg38", help="Genome version: hg19 or hg38")
7405 08 Nov 23 nicklas 14 parser <- add_argument(parser, "--penalty", default=70, help="Penalty parameter to ascat.aspcf() method")
7341 11 Sep 23 nicklas 15 parser <- add_argument(parser, "--tau", default=NA, type="numeric", help="Parameter used in the test for allelic imbalance. The meaning depends on the selected test.")
7278 10 Aug 23 nicklas 16 parser <- add_argument(parser, "--rho", default=NA, type="numeric", help="Optional argument to override ASCAT optimization and supply rho parameter for ascat.runAscat() method")
7278 10 Aug 23 nicklas 17 parser <- add_argument(parser, "--psi", default=NA, type="numeric", help="Optional argument to override ASCAT optimization and supply psi parameter for ascat.runAscat() method")
7405 08 Nov 23 nicklas 18 parser <- add_argument(parser, "--imbalanceTest", default="bimodality_coefficient", help="Selects which test for allelic imbalance that is used in ascat.aspcf(): mad_segment, legacy or bimodality_coefficient")
7278 10 Aug 23 nicklas 19 parser <- add_argument(parser, "--progress", default="progress", help="Path to file for reporting progress (Reggie-style)")
7341 11 Sep 23 nicklas 20 parser <- add_argument(parser, "--debug", flag=T, help="If set, ASCAT data objects are also saved to the output directory")
7278 10 Aug 23 nicklas 21 argv <- parse_args(parser)
7271 29 Jun 23 nicklas 22
7278 10 Aug 23 nicklas 23
7271 29 Jun 23 nicklas 24 library("ASCAT")
7271 29 Jun 23 nicklas 25
7271 29 Jun 23 nicklas 26 ## Load LogR and BAF files
7278 10 Aug 23 nicklas 27 cat('50 Loading LogR and BAF', file=argv$progress)
7271 29 Jun 23 nicklas 28 ascat.bc = ascat.loadData(
7278 10 Aug 23 nicklas 29     paste0(argv$inputDir, '/tumor_logr.txt'), 
7278 10 Aug 23 nicklas 30     paste0(argv$inputDir, '/tumor_baf.txt'), 
7278 10 Aug 23 nicklas 31     paste0(argv$inputDir, '/normal_logr.txt'), 
7278 10 Aug 23 nicklas 32     paste0(argv$inputDir, '/normal_baf.txt'), 
7278 10 Aug 23 nicklas 33     gender=argv$gender, 
7278 10 Aug 23 nicklas 34     genomeVersion=argv$genomeVersion
7271 29 Jun 23 nicklas 35   )
7372 06 Oct 23 nicklas 36 ## Save ASCAT data objects when debugging
7372 06 Oct 23 nicklas 37 if (argv$debug) save(ascat.bc, file=paste0(argv$outputDir, '/ascat.bc.Rdata'))
7271 29 Jun 23 nicklas 38
7271 29 Jun 23 nicklas 39 ## Make GC correction of LogR
7278 10 Aug 23 nicklas 40 cat('60 Performing GC correction', file=argv$progress)
7278 10 Aug 23 nicklas 41 ascat.bc = ascat.correctLogR(ascat.bc, GCcontentfile=argv$gcFile)
7372 06 Oct 23 nicklas 42 if (argv$debug) save(ascat.bc, file=paste0(argv$outputDir, '/ascat.bc.Rdata'))
7278 10 Aug 23 nicklas 43 ascat.plotRawData(ascat.bc, img.dir=argv$outputDir)
7372 06 Oct 23 nicklas 44
7278 10 Aug 23 nicklas 45 cat('70 Performing segmentation', file=argv$progress)
7341 11 Sep 23 nicklas 46 ascat.bc = ascat.aspcf(ascat.bc, penalty=argv$penalty, tau=argv$tau, imbalance.test=argv$imbalanceTest)
7372 06 Oct 23 nicklas 47 if (argv$debug) save(ascat.bc, file=paste0(argv$outputDir, '/ascat.bc.Rdata'))
7278 10 Aug 23 nicklas 48 ascat.plotSegmentedData(ascat.bc, img.dir=argv$outputDir)
7271 29 Jun 23 nicklas 49
7341 11 Sep 23 nicklas 50 # Save information about the allelic imbalance test
7341 11 Sep 23 nicklas 51 write.table(ascat.bc$ImbalanceTest, sep='\t', file=paste0(argv$outputDir, '/imbalance.test.txt'), quote=F, row.names=F, col.names=T)
7341 11 Sep 23 nicklas 52
7278 10 Aug 23 nicklas 53 cat('80 Solving ploidy and purity', file=argv$progress)
7278 10 Aug 23 nicklas 54 ascat.output = ascat.runAscat(ascat.bc, img.dir=argv$outputDir, write_segments=T, gamma=1, rho_manual=argv$rho, psi_manual=argv$psi)
7372 06 Oct 23 nicklas 55 if (argv$debug) save(ascat.output, file=paste0(argv$outputDir, '/ascat.output.Rdata'))
7271 29 Jun 23 nicklas 56
7405 08 Nov 23 nicklas 57 ## Generate metrics
7405 08 Nov 23 nicklas 58 ascat.m = ascat.metrics(ascat.bc, ascat.output)
7405 08 Nov 23 nicklas 59
7405 08 Nov 23 nicklas 60 ## Add some more information
7287 15 Aug 23 nicklas 61 purity = ascat.output$aberrantcellfraction[1]
7280 14 Aug 23 nicklas 62 ploidy = ascat.output$ploidy[1]
7421 14 Nov 23 nicklas 63 if (!is.null(ploidy) && !is.null(purity))
7421 14 Nov 23 nicklas 64 {
7421 14 Nov 23 nicklas 65   ascat.m$normal_contamination = 2*(1-purity)/(2*(1-purity)+purity*ploidy)
7421 14 Nov 23 nicklas 66   ascat.m$non_aberrant = ifelse(is.null(ascat.output$nonaberrantarrays), 0, 1)
7421 14 Nov 23 nicklas 67 }
7405 08 Nov 23 nicklas 68 write.table(t(ascat.m), paste0(argv$outputDir, '/stats.txt'), row.names=T, col.names=c('ascat.metrics\tvalue'), quote=F, sep="\t")
7271 29 Jun 23 nicklas 69

Generated on Tue Nov 21 03:33:45 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)