I am using this pipeline to make a vcf file. but m having problem at step 4. bcftools view -bvcg my-raw.bcf > my-var.bcf
it says -b is not a valid parameter.??? what to do ?
Convert SAM to BAM for sorting
samtools view -S -b my.sam > my.bam
Sort BAM for SNP calling
samtools sort my.bam my-sorted
Run 'mpileup' to generate VCF format
samtools mpileup -g -f my.fasta my-sorted-1.bam my-sorted-2.bam my-sorted-n.bam > myraw.bcf
Call SNPs...
bcftools view my.var.bcf | vcfutils.pl varFilter - > my.var-final.vcf
•
written
by
hafiz.talhamalik
For calling variants from assembled genomes, relative to a reference, you need an entirely different approach.
Take for example a look here: https://github.com/lh3/minimap2/blob/master/cookbook.md#asm-var Using minimap2 for alignment and paftools for variant calling.
oh okz.. Thanks for this. let me go through this. so in a nutshell you are suggesting samtool's aprroach is not good for this kind of analysis ?
You are aligning an assembly, so the aligner you used before is probably not the best for this job. Furthermore, you will have a coverage of 1 on most locations in the genome, which is of course not enough to call variants using pileup as you would for fastq reads.
Log in to answer this question.
Step 1: check if you have a recent bcftools version
Step 2: if it is the recent version, check if the parameter used are still valid by typing
bcftools view --helpIt looks like hafiz.talhamalik is following an outdated tutorial in variant calling using samtools/bcftools. In the recent versions this looks like this:
fin swimmer
yes my tutorial is old but the command you gave is giving an empty vcf file.
Hello again,
what do you mean by "empty". Just the vcf header or really a blank file?
fin swimmer
completely blank one. no header, nothing
Have you tried running them separately??
set the
--ploidyif required in the variant calling step.yes i tried them separately but still empty file.
Hi, have you checked the sorted bam file? The command you mentioned for sorting is not correct. You can either redirect the ouput using
>or can use the-ooption. Also index the refernece fasta file with samtools in the same directory.yeah I tried that. sorted bam file is fine.
Could you please post the output of these commands:
$ samtools --version$ bcftools --version$ samtools view -H my-sorted.bam$ samtools flagstat my-sorted.bam$ cat ref.fasta.faifin swimmer
argh, it should be
samtools flagstat(without s). I've edited it.This looks like you're trying to align a fasta sequence and not fastq. AFAIK
bwacan do this. But I'm not sure ifbcftoolscan call variants from this file, because the base quality values are missing and ifsequence.fastacontains only one sequence you only have one read per position in the alignment file.So could you please tell us more about what you are trying to do?
fin swimmer
I have downloaded all the fasta sequence of salmonella from ncbi and now i want to compare snp's from those fasta files with my sequenced data's vcf file.
So you have downloaded assembled genomes rather than separate reads?
You can do variant calling from those but it would have been absolutely important if you mention things like this in your first post.
yes assembled genomes.
You have wasted a lot of your and our time by not telling us that earlier. I have added an answer to this thread.
oho. M really sorry for that..!