needlestack icon indicating copy to clipboard operation
needlestack copied to clipboard

`--output_vcf` parameter doesn't work when including a relative path

Open Fer020707 opened this issue 2 years ago • 6 comments

Hi, i am trying to run the program using:

nextflow run iarcbioinfo/needlestack -with-docker --bed bed/SOPHIA_targetregions_hg19.bed --input_bams bam/SOPHIA/ --ref ref/ucsc.hg19.fasta --output_vcf output/SOPHIA_needlestack.vcf

The bed file is sorted, the bam files are in the corresponding folder with their respective .bai, and in the reference folder are the .fasta, .fai, .dict, .amb, .ann, .bwt, . pac and .sa.

But I get an error in the execution process. I would appreciate your response, thank you.

Error message:

executor >  local (3)
[76/9e2dc1] process > bed                                         [100%] 1 of 1 ✔
[18/a87fe0] process > split_bed                                   [100%] 1 of 1 ✔
[91/4ac861] process > mpileup2vcf (chr1_17345370-chrX_14883636... [  0%] 0 of 1
[-        ] process > collect_vcf_result                          -
Error executing process > 'mpileup2vcf (chr1_17345370-chrX_14883636_regions)'

Caused by:
  Process `mpileup2vcf (chr1_17345370-chrX_14883636_regions)` terminated with an error exit status (1)

Command executed:

  for cur_bam in BAM/*.bam
  do
      if [ "BAM" == "FILE" ]; then
          # use bam file name as sample name
          bam_file_name=$(basename "${cur_bam%.*}")
          # remove whitespaces from name
          SM1="$(echo -e "${bam_file_name}" | tr -d '[[:space:]]')"
          SM2="$(echo -e "${bam_file_name}" | tr -d '[[:space:]]')"
      else
          # get bam file names
   bam_file_name=$(basename "${cur_bam%.*}")
   # remove whitespaces from name
   SM1="$(echo -e "${bam_file_name}" | tr -d '[[:space:]]')"
   # extract sample name from bam file read group info field
   SM2=$(samtools view -H $cur_bam | grep "^@RG" | tail -n1 | sed "s/.*SM:\([^	]*\).*/\1/" | tr -d '[:space:]')
      fi
  
  
      printf "$SM1	$SM2\n" >> names.txt
  done
  
  if [ "FALSE" != "FALSE" ]; then
      abs_pairs_file=$(readlink -f FALSE)
  else
      abs_pairs_file="FALSE"
  fi
  set -o pipefail
  i=1

  ...
  
  { while read bed_line; do
      samtools mpileup --fasta-ref ucsc.hg19.fasta --region $bed_line --ignore-RG --min-BQ 13 --min-MQ 0 --max-idepth 1000000 --max-depth 50000 BAM/*.bam | sed 's/		/	*	*/g'
      i=$((i+1))
  done < chr1_17345370-chrX_14883636_regions
  } | mpileup2readcounts 0 -5 false 3 0 | Rscript /home/fer/.nextflow/assets/iarcbioinfo/needlestack/bin/needlestack.r --pairs_file=${abs_pairs_file} --source_path=/home/fer/.nextflow/assets/iarcbioinfo/needlestack/bin/ --out_file=chr1_17345370-chrX_14883636_regions.vcf --fasta_ref=ucsc.hg19.fasta --input_bams=BAM/ --ref_genome=null --GQ_threshold=50 --min_coverage=30 --min_reads=3 --min_af=0 --SB_type=SOR --SB_threshold_SNV=100 --SB_threshold_indel=100 --output_all_SNVs=false --do_plots=ALL --do_alignments=false --plot_labels=true --add_contours=true --extra_rob=false --min_af_extra_rob=0.2 --min_prop_extra_rob=0.1 --max_prop_extra_rob=0.5 --afmin_power=-1 --sigma=0.1

Command exit status:
  1

...

Command error:
  [W::bam_hdr_read] EOF marker is absent. The input is probably truncated
  [W::bam_hdr_read] EOF marker is absent. The input is probably truncated
  ...
  [warning] samtools mpileup option `L` is functional, but deprecated. Please switch to using bcftools mpileup in future.
  [W::bam_hdr_read] EOF marker is absent. The input is probably truncated
  [W::bam_hdr_read] EOF marker is absent. The input is probably truncated
..
  [mpileup] fail to load index for BAM/ISEM114.recal.reads.bam
  Error : Pileup file is empty

Work dir:
  /home/fer/Documents/work/91/4ac861e494a55ec631f77ea31b0ce6

Tip: you can try to figure out what's wrong by changing to the process work dir and showing the script file named `.command.sh`

Fer020707 avatar Jan 23 '22 00:01 Fer020707