samtools mpileup: error reading from input file

General
1

Dear all,

I am writing a code to analyse yeast sequencing data using R. I was able to do the alignment using BWA and thus obtain a .bam file and its .bai index using samtools. My next step is to perform variant calling using bcftools. Here is my command line and the associated error:

vcf_file <- "pathway/to/stock/variants.vcf"
system(paste("samtools mpileup -f", paste(fasta_dir, "chr*.fa", sep = ""), merged_bam_file,
"| bcftools call --ploidy 1 -mv >", vcf_file))

[mpileup] 17 samples in 17 input files
samtools mpileup: error reading from input file
Failed to read from standard input: unknown file type

  • paste(fasta_dir, "chr*.fa", sep = "") gives the pathway to the directory containing the 17 fasta files of the genome of reference (one for each chromosome), named chrI.fa, chrII.fa, ...

  • merged_bam_file is pathway to the alignment file I obtained and with which I want to perform the variant calling

I understand that the problem comes from the type of my files, but since the alignment file is BAM and my genome of reference is FASTA, I have absolutely no clue what to change to make it work.
Any advices?

Thank you for your help :slight_smile:

2

The sharp end of the script is system(), which sends a command to the external "shell" program to execute a non-R program. For example

> system("date")
#> [1]Sun Jul 16 18:00:38 PDT 2023

Where things go wrong can be anywhere that system() gets a command it can't execute. That's what happens when the first paste sends the name of the directory and the wild card smooshed together. This is corrected in the second.

vcf_file <- "~/projects/demo/variants.vcf"
fasta_dir <- "~/projects/demo/fasta_dir"
merged_bam_file <- "/project/demo/merger.xxx"

paste("samtools mpileup -f", 
      paste(fasta_dir, "chr*.fa", sep = ""), 
      merged_bam_file,
      "| bcftools call --ploidy 1 -mv >", 
      vcf_file)
#> [1] "samtools mpileup -f ~/projects/demo/fasta_dirchr*.fa /project/demo/merger.xxx | bcftools call --ploidy 1 -mv > ~/projects/demo/variants.vcf"
  

paste("samtools mpileup -f", 
      paste(fasta_dir, "/chr*.fa", sep = ""), 
      merged_bam_file,
      "| bcftools call --ploidy 1 -mv >", 
      vcf_file)
#> [1] "samtools mpileup -f ~/projects/demo/fasta_dir/chr*.fa /project/demo/merger.xxx | bcftools call --ploidy 1 -mv > ~/projects/demo/variants.vcf"

Created on 2023-07-16 with reprex v2.0.2

But there are at least a couple of packages that makes this easier. Rsamtools and {Samtools} ease some of the rough spots by relieving the necessity of having to construct elaborate raw commands