Biostar Beta. Not for public use.
Question: How can I get Allele Depth of my vcf file ?
0
Entering edit mode

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

ADD COMMENTlink 9 months ago apl00028 • 20 • updated 9 months ago zx8754 7.5k
Entering edit mode
0

how did you call the original bcf file ?

ADD REPLYlink 9 months ago
Pierre Lindenbaum
120k
Entering edit mode
0

I called it scaffold_filt_PV003_q10_variant.bcf, is it answer your question?

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

what is the tool used to call the variants ? samtools ? bcftools ?gatk ? what was the command line ?

ADD REPLYlink 9 months ago
Pierre Lindenbaum
120k
Entering edit mode
0

He means: how did you perform the variant calling? - GATK?; SAMtools?

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
0

I used samtools to get that file, the command line is

bcftools view -c -g input.bcf > output.vcf
ADD REPLYlink 9 months ago
apl00028
• 20
• updated 9 months ago
h.mon
25k
Entering edit mode
0

the question was:

how did you perform the variant calling

and not

how to you convert a vcf to bcf

ADD REPLYlink 9 months ago
Pierre Lindenbaum
120k
Entering edit mode
0

begining with a input file:

samtools mpileup -g -f reference.txt input.bam > output.bcf

I got the bcf then I convert a bcf to vcf with this command:

bcftools view input.bcf > output.vcf

And the output.vcf has not a record of Allele Depth. Is there another alternative way to get it? How?

ADD REPLYlink 9 months ago
apl00028
• 20
• updated 9 months ago
RamRS
21k
4
Entering edit mode

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?

ADD COMMENTlink 9 months ago Kevin Blighe 43k
Entering edit mode
0

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

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
1

your version of bcftools is just too old.

ADD REPLYlink 9 months ago
Pierre Lindenbaum
120k
Entering edit mode
0

Okey, I updated the bcftools version and I got another issue:

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

What is that mean? Thanks in advantage

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

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

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
1

GREAT! IT WORKS! THANKS!

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted.

Upvote|Bookmark|Accept

ADD REPLYlink 9 months ago
Pierre Lindenbaum
120k
Entering edit mode
0

Can I see those SNP with lower allele frequency? I have seen that vcf file records those SNP with higher frequency.

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

To be clear, to which parameter are you referring? - AF?

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
0

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.

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

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.

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
0

I get the same number of records using --pval-threshold 1.0 and whitout it, is there another way?

Thanks in advantage.

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

so, most probably these SNPs are sequecing errors, are not them?

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

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?

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
0

My BAM files were produced following these steps:

Firstly I did an alignment using bowtie2

  bowtie2 -f -x /cmv3.fasta /scaffold_nofilt_PV002.fa -S  /scaffold_nofilt_PV002.sam

After that I converted sam file in bam file

 samtools view -bS /scaffold_nofilt_PV002.sam > /scaffold_nofilt_PV002.bam

Then, I removed no mapped reads:

 samtools view -F 0x04 -b /scaffold_nofilt_PV003.bam >/scaffold_nofilt_PV003.bam

and I select only those data with a map quality of 10

 samtools view -q 10 -b /scaffold_nofilt_PV003.bam > /scaffold_nofilt_PV003.bam

Finally I sort my bam file and I index them:

 samtools sort /scaffold_nofilt_PV003.bam /scaffold_nofilt_PV003
 samtools index /scaffold_nofilt_PV003.bam
ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

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?

ADD REPLYlink 9 months ago
Kevin Blighe
43k
Entering edit mode
0

My results come from wild plant species of different host, I want to see the variation among cucumber mosaic virus RNA3 of differnt hosts or same host and different habitats.

ADD REPLYlink 9 months ago
apl00028
• 20
Entering edit mode
0

and my bcftools is the 1.9 the latest

ADD REPLYlink 9 months ago
apl00028
• 20

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.0