modify igv batch script to work directly on a MAF file, add example u…

…sage.
agorelick · Sep 1, 2017 · c0ad1db · c0ad1db
1 parent 2d0b697
commit c0ad1db
Showing 1 changed file with 31 additions and 13 deletions.
diff --git a/analysis/igv_batch_screenshots_Ed.R b/analysis/igv_batch_screenshots_Ed.R
@@ -2,6 +2,13 @@
 
 '%!in%' <- function(x,y)!('%in%'(x,y))
 
+# A modified version of Craig's script, with the idea being to parse a MAF file automatically. 
+# Note the options for selecting all non-synonymous mutations, or only LOF mutations.
+
+# Example usage:
+# Rscript igv_batch_screenshots_Ed.r -m MergedMAF_TCGATHCA.maf -f T -b /Users/ereznik/Documents/luna/ifs/e63data/makarovv/HCC/rawbam/ 
+# -o /Users/ereznik/Documents/mtimpact/results/hurthle/IGVscreenshots/IGVscript.txt -d /Users/ereznik/Documents/mtimpact/results/hurthle/IGVscreenshots/ -g b37
+
 catverbose <- function(...) {
   cat(format(Sys.time(), "%Y%m%d %H:%M:%S |"), ..., "\n")
 }
@@ -12,17 +19,17 @@ write_batch_script <- function(samples, out, dir, genome, squish=T) {
   for (id in ids) {
     cat('new')
     # cat(paste0('\ngenome ', genome))
-    dat <- samples[patient == id]
+    dat <- samples[which(samples$patient == id),]
     bams_to_load <- paste(dat$tumor,
                           dat$normal, 
                           sep = ",")
     cat(paste0('\nload ', bams_to_load))
     cat(paste0('\nsnapshotDirectory ', dir))
     cat('\nmaxPanelHeight 600\n')
     for (i in 1:nrow(dat)) {
-      cat(paste0('\ngoto chr', dat$chr[i], ':', dat$start[i]))
+      cat(paste0('\ngoto ', dat$chr[i], ':', dat$start[i]))
       cat(paste0('\nsort base'))
-      cat(paste0('\ngoto chr', dat$chr[i], ':', dat$start[i]-1,'-',dat$end[i]+1))
+      cat(paste0('\ngoto ', dat$chr[i], ':', dat$start[i]-1,'-',dat$end[i]+1))
       if (squish) {
         cat('\nsquish')
       }
@@ -40,27 +47,38 @@ if( ! interactive() ) {
   rm(tmp)
 
   parser=ArgumentParser()
-  parser$add_argument('-s', '--samples', type='character', help='Tab-delimited text file mapping sample IDs to bam files')
+  parser$add_argument('-m', '--maf', type='character',help='MAF-file, in format automatically generated by vcf2maf')
+  parser$add_argument('-b', '--bamdir',type='character',help='path to BAM files')
+  parser$add_argument('-f', '--filter',type='logical',help = 'if TRUE, only look at LOF mutations, otherwise all non-synonymous mutations')
   parser$add_argument('-o', '--out', type='character', help='Output batch script file')
   parser$add_argument('-d', '--dir', type='character', help='Screenshot directory')
   parser$add_argument('-g', '--genome', type='character', default='b37', help='Genome version')
   parser$add_argument('-sq', '--squish', action="store_true", default=T, help='Squish reads')
+
   args=parser$parse_args()
 
-  samples <- fread(args$samples)
+  maf <- read.csv(args$maf, header = TRUE,sep = '\t',stringsAsFactors = FALSE)
+  badmuts = c('Nonsense_Mutation','Missense_Mutation','Nonstop_Mutation','Translation_Start_Site',
+              'In_Frame_Del','Frame_Shift_Ins','Frame_Shift_Del','tRNA')
+  LOFmuts = c('Nonsense_Mutation','Frame_Shift_Ins','Frame_Shift_Del','tRNA')
+  if (args$filter){
+    maf = maf[which(maf$Variant_Classification %in% LOFmuts),]
+  }else{
+    maf = maf[which(maf$Variant_Classification %in% badmuts),]
+  }
+
+  # format of the file is sample_id,gene,chromosome,start_position,end_position,path to tumor bam, path to normal bam
+  samples = maf[,c('Tumor_Sample_Barcode','Hugo_Symbol','Chromosome','Start_Position','End_Position'),]
+  colnames(samples) = c('patient','gene','chr','start','end')
+  samples$tumor = paste(args$bamdir,maf$Tumor_Sample_Barcode,sep = '/')
+  samples$normal = paste(args$bamdir,maf$Matched_Norm_Sample_Barcode,sep = '/')
+
   out <- args$out
   dir <- args$dir
   genome <- args$genome
   #squish <- args$squish # set squish to false for now
-  squish <- FALSE
+  squish <- TRUE
 
-  setnames(samples, c('patient', 
-                      'gene', 
-                      'chr', 
-                      'start',
-		      'end', 
-                      'tumor', 
-                      'normal'))
   write_batch_script(samples, out, dir, 
                      genome, squish)