Biostar Beta. Not for public use.
Calling SNPs using bcftools mpileup
0
Entering edit mode
23 months ago
jaafari.omid • 40
@jaafari.omid42947

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

snp genome bcftools-mpileup WGS SNP calling • 994 views
0
Entering edit mode

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

0
Entering edit mode

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

cheers,

0
Entering edit mode

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

1
Entering edit mode

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

0
Entering edit mode

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

1
Entering edit mode
23 months ago
@finswimmer37605

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