I am attempting to get the Allele Depth (AD) from my vcf file.
I got this file from bcftools file using the function:
bcftools view -c -g scaffold_filt_PV003_q10_variant.bcf > scaffold_filt_PV003_q10_variant.vcf
but when I see the vcf generated in the column called format does not appear AD, only appear GT:PL:GQ I do not another alternative way to see Allele Depth.
Thanks in advance
Try this:
bcftools mpileup \
--redo-BAQ \
--min-BQ 30 \
--per-sample-mF \
--annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
-f reference.txt \
input.bam | \
bcftools call \
--multiallelic-caller \
--variants-only \
-Ov \
> output.vcf ;
I am not sure why your reference is TXT format - usually it is FASTA (and indexed with samtools faidx). Also, why the -g option?
well, I have follow your code:
samtools faidx cmv3.fasta to index the fasta file
bcftools mpileup --redo-BAQ --min-BQ 30 --per-sample-mF --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR -f cmv3.fasta scaffold_nofilt_PV158.bam bcftools call --multiallelic-caller --variants-only -Ov > scaffold_nofilt_PV158_q10_variant.vcf to create the vcf file
But I get only this message: [main] Unrecognized command.
What did I do wrong? Thanks
your version of bcftools is just too old.
You are missing the pipe symbol (|) just before bcftools call; however, this is my fault because I did not include it in my earlier comment (I have now added it).
If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted.
To be clear, to which parameter are you referring? - AF?
No, I mean Coveage allele-fraction https://bioinformatics-core-shared-training.github.io/intro-to-IGV/InspectingVariantsInIGV.html when I change the threshold to 0.05 in IGV I see more SNP that in VCF file I can not see. How can I change the bcftools paramethers to see them? Thanks in advantage.
If you are expecting variants at those low allele fractions and these are supposed to be germline variants, then that alludes to some problem with the sequencing run. If you wish to try to call these using BCFtools, then add the following parameter to bcftools call:
--pval-threshold 1.0
Note that this will result in a lot of junk being called, too.
BCFtools mpileup is excellent at calling germline variants that are expected at frequencies of 50% or 100%. Something else like GATK may better if you are expecting other types of variants.
Without seeing all information in your output VCF and having access to your BAM files, I cannot really comment further - sorry. What other information is in your VCF? - take a look at the QUAL score, for example. Also, how were your BAM files produced?
My BAM files were produced following these steps:
bowtie2 -f -x /cmv3.fasta /scaffold_nofilt_PV002.fa -S /scaffold_nofilt_PV002.sam
samtools view -bS /scaffold_nofilt_PV002.sam > /scaffold_nofilt_PV002.bam
samtools view -F 0x04 -b /scaffold_nofilt_PV003.bam >/scaffold_nofilt_PV003.bam
samtools view -q 10 -b /scaffold_nofilt_PV003.bam > /scaffold_nofilt_PV003.bam
samtools sort /scaffold_nofilt_PV003.bam /scaffold_nofilt_PV003
samtools index /scaffold_nofilt_PV003.bam
I am not sure that the removing of un-mapped reads is necessary, as these will not be included in the variant calling stage. I don't see any major issue by removing reads with MAPQ < 10.
The results can perhaps be explained by the experimental set-up: From where are the reads derived?; are you aligning reads to cytomegalovirus?
Log in to answer this question.
how did you call the original bcf file ?
I called it scaffold_filt_PV003_q10_variant.bcf, is it answer your question?
what is the tool used to call the variants ? samtools ? bcftools ?gatk ? what was the command line ?
He means: how did you perform the variant calling? - GATK?; SAMtools?
I used samtools to get that file, the command line is
the question was:
and not
begining with a input file:
I got the bcf then I convert a bcf to vcf with this command:
And the
output.vcfhas not a record of Allele Depth. Is there another alternative way to get it? How?