change calls

uzh · Dec 5, 2021 · 7163d04 · 7163d04
1 parent ec99455
commit 7163d04
Showing 1 changed file with 7 additions and 3 deletions.
diff --git a/R/app-gatkRnaHaplotyper.R b/R/app-gatkRnaHaplotyper.R
@@ -42,16 +42,18 @@ ezMethodGatkRnaHaplotyper = function(input=NA, output=NA, param=NA){
     ezSystem('mv local.bam withRg.bam')
   }
 
+  ezSystem("samtools index withRg.bam")
   cmd = paste(gatkCall, "SplitNCigarReads", 
               "-R", genomeSeq,
               "-I", "withRg.bam",
               "-L", exomeInterVals,
               ## this is the default! "-rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS",
               "-O splitNtrim.bam")
   ezSystem(cmd)
-
+  ezSystem("samtools index splitNtrim.bam")
+
   #BaseRecalibration is done only if known sites are available
-  if(param$knownSitesAvailable){
+  if(file.exists(dbsnpFile)){
     cmd = paste(gatkCall, "BaseRecalibrator", 
                 "-R", genomeSeq,
                 "-I splitNtrim.bam",
@@ -68,13 +70,15 @@ ezMethodGatkRnaHaplotyper = function(input=NA, output=NA, param=NA){
                 "--add-output-sam-program-record",
                 "--bqsr-recal-file recal.table",
                 "--output recal.bam")
+    ezSystem(cmd)
   } else {
     ezSystem('mv splitNtrim.bam recal.bam')
   }
+  ezSystem("samtools index recal.bam")
 
   ########### haplotyping
   outputFile = basename(output$getColumn("GVCF"))#  paste0(sampleName, ".g.vcf.gz")
-  cmd = paste(gatkCal,'HaplotypeCaller',
+  cmd = paste(gatkCall,'HaplotypeCaller',
               "-R", genomeSeq,
               "-I recal.bam",
               "-L", exomeInterVals,