Calling SNPs using bcftools mpileup

Question

[jaafari.omid]

Hi dears all, I am going to do SNP calling using bcftools mpileup for WGS data on fish. Before SNP calling I did this preprocessing as follow:

samtools sort -n F0-F06-10.bam F0-F06-10ns.bam

Add ms and MC tags for markdup to use later

samtools fixmate -m F0-F06-10ns.bam F0-F06-10-fixmate.bam

Markdup needs position order

samtools sort -o F0-F06-10PS.bam F0-F06-10-fixmate.bam

Finally mark duplicates

samtools markdup F0-F06-10PS.bam F0-F06-10-markdup.bam

Then using samtools faidx the fasta reference file got indexed and bcftools mpileup was used for SNP calling. For the last step (bcftools mpileup) I got involved with an error.

samtools faidx /home2/omid/alignment/ref/GCF_001858045.2_O_niloticus_UMD_NMBU_genomic.fa

bcftools mpileup -Ou -f GCF_001858045.2_O_niloticus_UMD_NMBU_genomic.fa.fai F0-F06-10-markdup.bam

error: [E::fai_build_core] Format error, unexpected "N" at line 1

I will be so grateful if any one can help me to solve this problem.

Cheers

[finswimmer]

Hello,

the filenames used in the bcftools command doesn't match the filenames you present before. So please show us the exact command you've used.

Thanks!

fin swimmer

[jaafari.omid]

Hello and Thanks for the comment. I have updated the commands to the exact ones that I used.

cheers,

[jaafari.omid]

You are right, When I changed the reference file to the original reference fasta file I did not see that error but the terminal got unnormal after running the command. I have attached a screen shot of my terminal. can you please let me know what can be the problem?

Cheers, Omid

[finswimmer]

Hello jaafari.omid ,

please use the ADD REPLYbutton below a post you like to reply to.

In your bcftools you define that the output type is uncompressed bcf. But you neither define an output file (-o output) nor redirect it to a file (>output).

Also notice that the results of mpileup are not your final variants. This step just collect metrics about each covered site ,which are then used by bcftools call to find sites which are investigated for variants. So if you goal is finding variants, you have to do something like this:

$ bcftools mpileup -Ou -f ref.fa aln.bam | bcftools call -Ov -mv > output.vcf

fin swimmer

[jaafari.omid]

I am really grateful for your great help. Now it worked fine. Cheers,

Answer 1

[finswimmer]

I'm quite sure that the reference file you are using in your bcftools command is the index file of your reference and not the reference file itself.

The index file must be located in the same location like the reference file. The file you specify with the -f parameter must be the reference genome file.

fin swimmer