1

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 :)

Stibu
  • 15,166
  • 6
  • 57
  • 71
Juliette G.
  • 11
  • 1
  • You're 100% sure that the FASTA files are properly formatted? I personally hate trying to execute system commands from within R, I'd rather do it on the command line to be sure that I haven't screwed something up, but you might also just check that the output of your `paste` command is giving you what you want. Assuming it's not something trivial like that, I would suggest visiting [biostars.org](http://biostars.org), which really specializes in bioinformatics help per se. – C. Murtaugh Jul 15 '23 at 04:06
  • Thanks a lot for your answer, I will try to post my question/find topics on the website you suggested.Regarding my FASTA files, I opened them and check the format and they look fine to me. I already checked also the output of paste and the pathway are correct. – Juliette G. Jul 16 '23 at 16:49
  • 1
    I was able to find the cause of the problem! The results of the samtools mpileup cannot be used by bcftools call, instead, I used bcftools mpileup and it worked well. I leave this here in case other people are interested. – Juliette G. Jul 20 '23 at 09:31

0 Answers0