Merge pull request #85 from CuppenResearch/release-v2.4.0

Release v2.4.0

Author: rernst

LICENSE

@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2016 Cuppen Research
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+

README.md

@@ -47,6 +47,7 @@ perl illumina_pipeline.pl /path/to/output_dir/settings.config>
 - [IGVtools](https://www.broadinstitute.org/igv/igvtools)
 - [Contra](http://contra-cnv.sourceforge.net/)
 - [FREEC](http://bioinfo-out.curie.fr/projects/freec/)
+- [QDNAseq](https://github.com/ccagc/QDNAseq) (used tag: v1.9.2-HMF.1)
 - [Varscan](http://varscan.sourceforge.net/)
 - [Strelka](https://sites.google.com/site/strelkasomaticvariantcaller/)
 - [Freebayes](https://github.com/ekg/freebayes)
@@ -338,6 +339,14 @@ CNVCHECK_MEM	maximum_memory
 CNV_MODE	sample_control
 CNV_TARGETS	/path/to/target.bed | Optional, use for targeted data e.g. exome.
 
+## QDNASEQ
+CNV_QDNASEQ	yes
+QDNASEQ_QUEUE	queue_name
+QDNASEQ_TIME	estimated runtime
+QDNASEQ_THREADS	number of threads
+QDNASEQ_MEM	maximum memory
+QDNASEQ_PATH	/hpc/local/CentOS7/cog_bioinf/QDNAseq_v1.9.2-HMF.1
+
 ## Contra
 CNV_CONTRA	yes/no
 CONTRA_THREADS	number_of_threads
@@ -435,5 +444,6 @@ NIPT_MASTER_THREADS	number_of_threads
 CHECKING_QUEUE	queue_name
 CHECKING_TIME	estimated runtime
 CHECKING_THREADS	number_of_threads
+CHECKING_RM	list,of,files,to,remove
 
 ```

illumina_check.pm

@@ -306,9 +306,11 @@ sub runCheck {
     print BASH "else\n";
     print BASH "\techo \"The pipeline completed successfully. The md5sum file will be created.\">>$logFile\n";
 
-    # Remove all tmp folders and empty logs except .done files if pipeline completed successfully
-    print BASH "\trm -r $opt{OUTPUT_DIR}/tmp\n";
-    print BASH "\trm -r $opt{OUTPUT_DIR}/*/tmp\n";
+    # Remove files/folders based on CHECKING_RM and empty logs except .done files if pipeline completed successfully
+    my @removeFiles = split(",", $opt{CHECKING_RM});
+    foreach my $removeFile (@removeFiles){
+	print BASH "\tfind $opt{OUTPUT_DIR} -name $removeFile -print0 | xargs -0 rm -r -- \n";
+    }
     print BASH "\tfind $opt{OUTPUT_DIR}/logs -size 0 -not -name \"*.done\" -delete\n";
     print BASH "\tfind $opt{OUTPUT_DIR}/*/logs -size 0 -not -name \"*.done\" -delete\n";
     print BASH "\tfind $opt{OUTPUT_DIR}/somaticVariants/*/logs -size 0 -not -name \"*.done\" -delete\n";
@@ -321,12 +323,6 @@ sub runCheck {
 	    }
 	}
     }
-    if($opt{SOMATIC_VARIANTS} eq "yes" && $opt{SOMVAR_VARSCAN} eq "yes"){
-	foreach my $sample (@{$opt{SOMATIC_SAMPLES_UNIQ}}){
-		print BASH "\trm $opt{OUTPUT_DIR}/$sample/mapping/$sample*.pileup.gz\n";
-		print BASH "\trm $opt{OUTPUT_DIR}/$sample/mapping/$sample*.pileup.gz.tbi\n";
-	}
-    }
     # Send email.
     print BASH "\tmail -s \"IAP DONE $runName\" \"$opt{MAIL}\" < $logFile\n";
 

illumina_copyNumber.pm

@@ -151,14 +151,27 @@ sub runCopyNumberTools {
 		    my $freec_job = runFreec($sample, $sample_out_dir, $sample_job_dir, $sample_log_dir, $sample_bam, "", \@running_jobs, \%opt);
 		    if($freec_job){push(@cnv_jobs, $freec_job)};
 		}
+		if($opt{CNV_QDNASEQ} eq "yes"){
+		    print "\n###SCHEDULING QDNASEQ####\n";
+		    my $qdnaseq_job = runqDNAseq($sample, $sample_out_dir, $sample_job_dir, $sample_log_dir, $sample_bam, \@running_jobs, \%opt);
+		    if($qdnaseq_job){push(@cnv_jobs, $qdnaseq_job)};
+		}
+		
+		
 	    ## Check copy number analysis
 	    my $job_id = "CHECK_".$sample."_".get_job_id();
 	    my $bash_file = $sample_job_dir."/".$job_id.".sh";
 
 	    open CHECK_SH, ">$bash_file" or die "cannot open file $bash_file \n";
 	    print CHECK_SH "#!/bin/bash\n\n";
 	    print CHECK_SH "echo \"Start Check\t\" `date` `uname -n` >> $sample_log_dir/check.log\n\n";
-	    print CHECK_SH "if [[ -f $sample_log_dir/freec.done ]]\n";
+	    if($opt{CNV_QDNASEQ} eq "yes" && $opt{CNV_FREEC} eq "yes"){
+		print CHECK_SH "if [ -f $sample_log_dir/freec.done -a -f $sample_log_dir/qdnaseq.done ]\n";
+	    } elsif ($opt{CNV_QDNASEQ} eq "yes" && $opt{CNV_FREEC} eq "no") {
+		print CHECK_SH "if [ -f $sample_log_dir/qdnaseq.done ]\n";
+	    } elsif ($opt{CNV_QDNASEQ} eq "no" && $opt{CNV_FREEC} eq "yes") {
+		print CHECK_SH "if [ -f $sample_log_dir/freec.done ]\n";
+	    }
 	    print CHECK_SH "then\n";
 	    print CHECK_SH "\ttouch $sample_log_dir/$sample.done\n";
 	    print CHECK_SH "fi\n\n";
@@ -177,6 +190,62 @@ sub runCopyNumberTools {
 }
 
 ### Copy number analysis tools
+sub runqDNAseq {
+    ###
+    # Run qDNAseq
+    ###
+    my ($sample_name, $out_dir, $job_dir, $log_dir, $sample_bam, $running_jobs, $opt) = (@_);
+    my @running_jobs = @{$running_jobs};
+    my %opt = %{$opt};
+    
+    # Skip qdnseq if done file exists
+    if (-e "$log_dir/qdnaseq.done"){
+	print "WARNING: $log_dir/qdnaseq.done exists, skipping \n";
+	return;
+    }
+    
+    ## Create qdnaseq output directory
+    my $qdnaseq_out_dir = "$out_dir/qdnaseq";
+    if(! -e $qdnaseq_out_dir){
+	make_path($qdnaseq_out_dir) or die "Couldn't create directory: $qdnaseq_out_dir\n";
+    }
+    
+    my $command = "Rscript $opt{IAP_PATH}/scripts/run_QDNAseq.R -qdnaseq_path $opt{QDNASEQ_PATH} ";
+    $command .= "-s $sample_name ";
+    $command .= "-b $sample_bam ";
+
+    ## Create qDNAseq bash script
+    my $job_id = "QDNASEQ_".$sample_name."_".get_job_id();
+    my $bash_file = $job_dir."/".$job_id.".sh";
+
+    open QDNASEQ_SH, ">$bash_file" or die "cannot open file $bash_file \n";
+    print QDNASEQ_SH "#!/bin/bash\n\n";
+    print QDNASEQ_SH "if [ -s $sample_bam ]\n";
+    print QDNASEQ_SH "then\n";
+    print QDNASEQ_SH "\techo \"Start QDNASEQ\t\" `date` \"\t $sample_bam\t\" `uname -n` >> $log_dir/qdnaseq.log\n";
+    print QDNASEQ_SH "\tcd $qdnaseq_out_dir\n";
+    print QDNASEQ_SH "\t$command\n";
+    print QDNASEQ_SH "\tif [ -s $qdnaseq_out_dir/$sample_name*.vcf ]\n";
+    print QDNASEQ_SH "\tthen\n";
+    print QDNASEQ_SH "\t\ttouch $log_dir/qdnaseq.done\n";
+    print QDNASEQ_SH "\tfi\n";
+    print QDNASEQ_SH "\techo \"End QDNASEQ\t\" `date` \"\t $sample_bam\t\" `uname -n` >> $log_dir/qdnaseq.log\n";
+    print QDNASEQ_SH "else\n";
+    print QDNASEQ_SH "\techo \"ERROR:  Input bam files do not exist.\" >> $log_dir/qdnaseq.log\n";
+    print QDNASEQ_SH "fi\n";
+    close QDNASEQ_SH;
+
+    ## Run job
+    my $qsub = &qsubTemplate(\%opt,"QDNASEQ");
+    if ( @running_jobs ){
+	system "$qsub -o $log_dir -e $log_dir -N $job_id -hold_jid ".join(",",@running_jobs)." $bash_file";
+    } else {
+	system "$qsub -o $log_dir -e $log_dir -N $job_id $bash_file";
+    }
+    return $job_id;
+}
+
+
 sub runFreec {
     ###
     # Run freec and plot result
@@ -185,7 +254,7 @@ sub runFreec {
     my @running_jobs = @{$running_jobs};
     my %opt = %{$opt};
 
-    ## Skip Contra if .done file exist
+    ## Skip Freec if .done file exist
     if (-e "$log_dir/freec.done"){
 	print "WARNING: $log_dir/freec.done exists, skipping \n";
 	return;

illumina_pipeline.pl

@@ -594,7 +594,7 @@ sub checkConfig{
 	    if(! $opt{FREEBAYES_TIME}){ print "ERROR: No FREEBAYES_TIME option found in config files.\n"; $checkFailed = 1; }
 	    if(! $opt{FREEBAYES_SETTINGS}){ print "ERROR: No FREEBAYES_SETTINGS option found in config files.\n"; $checkFailed = 1; }
 	    if(! $opt{FREEBAYES_SOMATICFILTER}){ print "ERROR: No FREEBAYES_SOMATICFILTER option found in config files.\n"; $checkFailed = 1; }
-	    if(! $opt{FREEBAYES_GERMLINEFILTER}){ print "ERROR: No FREEBAYES_GERMLINEFILTER option found in config files.\n"; $checkFailed = 1; }
+	    #if(! $opt{FREEBAYES_GERMLINEFILTER}){ print "ERROR: No FREEBAYES_GERMLINEFILTER option found in config files.\n"; $checkFailed = 1; }
 	}
 	if(! $opt{SOMVAR_MUTECT}){ print "ERROR: No SOMVAR_MUTECT option found in config files.\n"; $checkFailed = 1; }
 	if($opt{SOMVAR_MUTECT} eq "yes"){
@@ -656,6 +656,14 @@ sub checkConfig{
 	    if(! $opt{FREEC_WINDOW}){ print "ERROR: No FREEC_WINDOW option found in config files.\n"; $checkFailed = 1; }
 	    if(! $opt{FREEC_TELOCENTROMERIC}){ print "ERROR: No FREEC_TELOCENTROMERIC option found in config files.\n"; $checkFailed = 1; }
 	}
+	if(! $opt{CNV_QDNASEQ}){ print "ERROR: No CNV_QDNASEQ  in config files.\n"; $checkFailed = 1; }
+	if($opt{CNV_QDNASEQ} eq "yes"){
+	    if(! $opt{QDNASEQ_PATH}){ print "ERROR: No QDNASEQ_PATH option found in config files.\n"; $checkFailed = 1; }
+	    if(! $opt{QDNASEQ_QUEUE}){ print "ERROR: No QDNASEQ_QUEUE option found in config files.\n"; $checkFailed = 1; }
+	    if(! $opt{QDNASEQ_THREADS}){ print "ERROR: No QDNASEQ_THREADS option found in config files.\n"; $checkFailed = 1; }
+	    if(! $opt{QDNASEQ_MEM}){ print "ERROR: No QDNASEQ_MEM option found in config files.\n"; $checkFailed = 1; }
+	    if(! $opt{QDNASEQ_TIME}){ print "ERROR: No QDNASEQ_TIME option found in config files.\n"; $checkFailed = 1; }
+	}
     }
     ## SV_CALLING
     if($opt{SV_CALLING} eq "yes"){
@@ -793,6 +801,7 @@ sub checkConfig{
         if(! $opt{CHECKING_THREADS}){ print "ERROR: No CHECKING_THREADS found in .ini file\n"; $checkFailed = 1; }
         if(! $opt{CHECKING_MEM}){ print "ERROR: No CHECKING_MEM found in .ini file\n"; $checkFailed = 1; }
         if(! $opt{CHECKING_TIME}){ print "ERROR: No CHECKING_TIME found in .ini file\n"; $checkFailed = 1; }
+        if(! $opt{CHECKING_RM}){ print "ERROR: No CHECKING_RM found in .ini file\n"; $checkFailed = 1; }
     }
 
     if ($checkFailed) { 

illumina_somaticVariants.pm

@@ -215,10 +215,13 @@ sub runStrelka {
     print STRELKA_SH "\techo \"Start Strelka\t\" `date` \"\t $sample_ref_bam \t $sample_tumor_bam\t\" `uname -n` >> $log_dir/strelka.log\n\n";
 
     # Run Strelka
-    print STRELKA_SH "\t$opt{STRELKA_PATH}/bin/configureStrelkaWorkflow.pl --tumor $sample_tumor_bam --normal $sample_ref_bam --ref $opt{GENOME} --config $opt{STRELKA_INI} --output-dir $strelka_out_dir\n\n";
+    if (-e $strelka_out_dir){
+	print STRELKA_SH "\trm -r $strelka_out_dir \n";
+    }
+    print STRELKA_SH "\t$opt{STRELKA_PATH}/bin/configureStrelkaWorkflow.pl --tumor $sample_tumor_bam --normal $sample_ref_bam --ref $opt{GENOME} --config $opt{IAP_PATH}/$opt{STRELKA_INI} --output-dir $strelka_out_dir\n\n";
 
     print STRELKA_SH "\tcd $strelka_out_dir\n";
-    print STRELKA_SH "\tmake -j 8\n\n";
+    print STRELKA_SH "\tmake -j $opt{STRELKA_THREADS}\n\n";
 
     # Check strelka completed
     print STRELKA_SH "\tif [ -f $strelka_out_dir/task.complete ]\n";
@@ -542,26 +545,13 @@ sub runFreeBayes {
     # Uniqify freebayes output 
     print FREEBAYES_SH "\tuniq $freebayes_out_dir/$sample_tumor_name.vcf > $freebayes_out_dir/$sample_tumor_name.uniq.vcf\n";
     print FREEBAYES_SH "\tmv $freebayes_out_dir/$sample_tumor_name.uniq.vcf $freebayes_out_dir/$sample_tumor_name.vcf\n\n";
-
-    # get sample ids
-    print FREEBAYES_SH "\tsample_R=`grep -P \"^#CHROM\" $freebayes_out_dir/$sample_tumor_name.vcf | cut -f 10`\n";
-    print FREEBAYES_SH "\tsample_T=`grep -P \"^#CHROM\" $freebayes_out_dir/$sample_tumor_name.vcf | cut -f 11`\n\n";
-
-    # annotate somatic and germline scores
-    print FREEBAYES_SH "\t$opt{VCFLIB_PATH}/vcfsamplediff VT \$sample_R \$sample_T $freebayes_out_dir/$sample_tumor_name.vcf > $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf\n";
-    print FREEBAYES_SH "\tsed -i 's/SSC/FB_SSC/' $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf\n"; # to resolve merge conflicts with varscan vcfs
-    print FREEBAYES_SH "\tgrep -P \"^#\" $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf > $freebayes_out_dir/$sample_tumor_name\_germline.vcf\n";
-    print FREEBAYES_SH "\tgrep -P \"^#\" $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf > $freebayes_out_dir/$sample_tumor_name\_somatic.vcf\n";
-    print FREEBAYES_SH "\tgrep -i \"VT=germline\" $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf >> $freebayes_out_dir/$sample_tumor_name\_germline.vcf\n";
-    print FREEBAYES_SH "\tgrep -i \"VT=somatic\" $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf >> $freebayes_out_dir/$sample_tumor_name\_somatic.vcf\n";
-    print FREEBAYES_SH "\trm $freebayes_out_dir/$sample_tumor_name\_VTannot.vcf\n\n";
-
-    # Filter
-    print FREEBAYES_SH "\tcat $freebayes_out_dir/$sample_tumor_name\_somatic.vcf | java -Xmx".$opt{FREEBAYES_MEM}."G -jar $opt{SNPEFF_PATH}/SnpSift.jar filter \"$opt{FREEBAYES_SOMATICFILTER}\" > $freebayes_out_dir/$sample_tumor_name\_somatic_filtered.vcf\n";
-    print FREEBAYES_SH "\tcat $freebayes_out_dir/$sample_tumor_name\_germline.vcf | java -Xmx".$opt{FREEBAYES_MEM}."G -jar $opt{SNPEFF_PATH}/SnpSift.jar filter \"$opt{FREEBAYES_GERMLINEFILTER}\" > $freebayes_out_dir/$sample_tumor_name\_germline_filtered.vcf\n\n";
+    
+    # Filter freebayes vcf
+    print FREEBAYES_SH "python $opt{IAP_PATH}/scripts/filterFreebayes.py -v $freebayes_out_dir/$sample_tumor_name.vcf |";
+    print FREEBAYES_SH "java -Xmx".$opt{FREEBAYES_MEM}."G -jar $opt{SNPEFF_PATH}/SnpSift.jar filter \"$opt{FREEBAYES_SOMATICFILTER}\" > $freebayes_out_dir/$sample_tumor_name\_somatic_filtered.vcf \n";
 
     #Check freebayes completed
-    print FREEBAYES_SH "\tif [ -s $freebayes_out_dir/$sample_tumor_name\_somatic_filtered.vcf -a -s $freebayes_out_dir/$sample_tumor_name\_germline_filtered.vcf ]\n";
+    print FREEBAYES_SH "\tif [ -s $freebayes_out_dir/$sample_tumor_name\_somatic_filtered.vcf ]\n";
     print FREEBAYES_SH "\tthen\n";
     print FREEBAYES_SH "\t\t$rm_command\n";
     print FREEBAYES_SH "\t\ttouch $log_dir/freebayes.done\n\n"; ## Check on complete output!!!

scripts/filterFreebayes.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+
+"""
+filterFreebayes.py
+Applies BCBIO developed filters to the freebayes somatic output.
+2 main filters applied:
+1. LOD Tumor, LOD Normal > 3.5
+2. Freq of Tumor / Freq Normal > 2.7
+"""
+
+import argparse
+import re
+
+TUMOR_PARTS_INDEX = 10
+NORMAL_PARTS_INDEX = 9
+
+# Thresholds are like phred scores, so 3.5 = phred35
+TUMOR_THRESHOLD = 3.5
+NORMAL_THRESHOLD = 3.5
+
+# Code copied from BCBIO Freebayes pipeline.
+def customFilterFreebayes(vcf_file):
+    with open(vcf_file, "r") as f:
+        stripped_lines = (line.strip("\n") for line in f)
+        somatic_lines = (line for line in stripped_lines if _check_line(line))
+        print "\n".join(somatic_lines)
+
+### Filtering
+
+def _check_line(line):
+    parts = line.split("\t")
+    return line.startswith("#") or (_check_lods(parts, TUMOR_THRESHOLD, NORMAL_THRESHOLD) and _check_freqs(parts))
+
+def _check_lods(parts, tumor_thresh, normal_thresh):
+    """Ensure likelihoods for tumor and normal pass thresholds.
+
+    Skipped if no FreeBayes GL annotations available.
+    """
+    try:
+        gl_index = parts[8].split(":").index("GL")
+    except ValueError:
+        return True
+    try:
+        tumor_gls = [float(x) for x in parts[TUMOR_PARTS_INDEX].split(":")[gl_index].split(",") if x != "."]
+        if tumor_gls:
+            tumor_lod = max(tumor_gls[i] - tumor_gls[0] for i in range(1, len(tumor_gls)))
+        else:
+            tumor_lod = -1.0
+    # No GL information, no tumor call (so fail it)
+    except IndexError:
+        tumor_lod = -1.0
+    try:
+        normal_gls = [float(x) for x in parts[NORMAL_PARTS_INDEX].split(":")[gl_index].split(",") if x != "."]
+        if normal_gls:
+            normal_lod = min(normal_gls[0] - normal_gls[i] for i in range(1, len(normal_gls)))
+        else:
+            normal_lod = normal_thresh
+    # No GL inofmration, no normal call (so pass it)
+    except IndexError:
+        normal_lod = normal_thresh
+    return normal_lod >= normal_thresh and tumor_lod >= tumor_thresh
+
+def _check_freqs(parts):
+    """Ensure frequency of tumor to normal passes a reasonable threshold.
+
+    Avoids calling low frequency tumors also present at low frequency in normals,
+    which indicates a contamination or persistent error.
+    """
+    thresh_ratio = 2.7
+    try:  # FreeBayes
+        ao_index = parts[8].split(":").index("AO")
+        ro_index = parts[8].split(":").index("RO")
+    except ValueError:
+        ao_index, ro_index = None, None
+    try:  # VarDict
+        af_index = parts[8].split(":").index("AF")
+    except ValueError:
+        af_index = None
+    if af_index is None and ao_index is None:
+        raise NotImplementedError("Unexpected format annotations: %s" % parts[0])
+    def _calc_freq(item):
+        try:
+            if ao_index is not None and ro_index is not None:
+                ao = sum([int(x) for x in item.split(":")[ao_index].split(",")])
+                ro = int(item.split(":")[ro_index])
+                freq = ao / float(ao + ro)
+            elif af_index is not None:
+                freq = float(item.split(":")[af_index])
+        except (IndexError, ValueError, ZeroDivisionError):
+            freq = 0.0
+        return freq
+    tumor_freq, normal_freq = _calc_freq(parts[TUMOR_PARTS_INDEX]), _calc_freq(parts[NORMAL_PARTS_INDEX])
+    return normal_freq <= 0.001 or normal_freq <= tumor_freq / thresh_ratio
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(formatter_class=lambda prog: argparse.HelpFormatter(prog,max_help_position=100, width=200))
+
+    required_named = parser.add_argument_group("required named arguments")
+    required_named.add_argument("-v", "--vcf_file", help="path/to/file.vcf", required=True)
+
+    args = parser.parse_args()
+    customFilterFreebayes(args.vcf_file)

scripts/makeBAFplot.R

@@ -18,7 +18,6 @@ f2 <- function(x) sum(unlist(x)[2:length(x)])
 f1 <- function(x) sum(unlist(x))
 
 determine_peaks <- function(region) {
-  # ADD CHECK FOR NUMBER OF ROWS >= 50?
   quant <- quantile(region$baf)
   if(quant[2]==quant[4]) {
     # LOW COMPLEXITY
@@ -35,6 +34,10 @@ determine_peaks <- function(region) {
   tp<-turnpoints(ts_y)
   peaks <- data.frame(baf=d$x[tp$tppos], dens=d$y[tp$tppos])
   df <- data.frame(subset(peaks, dens>=2.0))
+  if (nrow(df)<2) {
+    # complex region
+    df <- data.frame(subset(peaks, dens>=1.0))
+  }
   df$CHROM <- region$CHROM[1]
   df$POS <- mean(region$POS)
   df$CALL <- NA