Question: VCF file without SNPs
0
gravatar for Aiah
4.2 years ago by
Aiah10
Aiah10 wrote:

Hi I have tried to generate a VCF file using this command

samtools mpileup -ugf A1163.fa.fai aP1.bam 2>stderr.out | bcftools view -bvcg - > aP1.raw.bcf 
bcftools view aP1.raw.bcf | perl vcfutils.pl varFilter -D2 > aP1.flt.vcf

and it generated a small file without SNPs I tiered changing the depth to 2 but still nothing any suggestions ?

here is the generated VCF file

##fileformat=VCFv4.1
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=DP4,Number=4,Type=Integer,Description="# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads">
##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same">
##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">
##INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=IS,Number=2,Type=Float,Description="Maximum number of reads supporting an indel and fraction of indel reads">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes for each ALT allele, in the same order as listed">
##INFO=<ID=G3,Number=3,Type=Float,Description="ML estimate of genotype frequencies">
##INFO=<ID=HWE,Number=1,Type=Float,Description="Chi^2 based HWE test P-value based on G3">
##INFO=<ID=CLR,Number=1,Type=Integer,Description="Log ratio of genotype likelihoods with and without the constraint">
##INFO=<ID=UGT,Number=1,Type=String,Description="The most probable unconstrained genotype configuration in the trio">
##INFO=<ID=CGT,Number=1,Type=String,Description="The most probable constrained genotype configuration in the trio">
##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=PC2,Number=2,Type=Integer,Description="Phred probability of the nonRef allele frequency in group1 samples being larger (,smaller) than in group2.">
##INFO=<ID=PCHI2,Number=1,Type=Float,Description="Posterior weighted chi^2 P-value for testing the association between group1 and group2 samples.">
##INFO=<ID=QCHI2,Number=1,Type=Integer,Description="Phred scaled PCHI2.">
##INFO=<ID=PR,Number=1,Type=Integer,Description="# permutations yielding a smaller PCHI2.">
##INFO=<ID=QBD,Number=1,Type=Float,Description="Quality by Depth: QUAL/#reads">
##INFO=<ID=RPB,Number=1,Type=Float,Description="Read Position Bias">
##INFO=<ID=MDV,Number=1,Type=Integer,Description="Maximum number of high-quality nonRef reads in samples">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias (v2) for filtering splice-site artefacts in RNA-seq data. Note: this version may be broken.">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# high-quality non-reference bases">
##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##INFO=<ID=PR,Number=1,Type=Integer,Description="# permutations yielding a smaller PCHI2.">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
snp • 3.1k views
ADD COMMENTlink modified 4.2 years ago • written 4.2 years ago by Aiah10
1

Can you first try to save the output of each step (samtools mpileup, bcftools view, bcftools view...>vcf)? Does the problem solely occur in the final step?

ADD REPLYlink written 4.2 years ago by H.Hasani920

I am not very expert, but why did you provide A1163.fa.fai instead of A1163.fa? It is true that in the manual it says "indexed", but this means that you have to provide the FASTA and also have it indexed, but NOT provide the index name. Also, posting the stderr.out error file would be useful.

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by Fabio Marroni2.6k

Thank you for your reply -I used the A1163.fa.fai file for the reason you mentioned. Are you saying that I should give the file a different name after indexing ?

-If I ran the command with A1163.fa I get the fallowing message: [bcf_sync] incorrect number of fields (0 != 5) at 0:0

  • I used the stderr.outin the command to avoid these errors: [mpileup] 1 samples in 1 input files [bcf_sync] incorrect number of fields (0 != 5) at 0:0 [afs] 0:0.000 <mpileup> Set max per-file depth to 8000 [bam_pileup_core] the input is not sorted (reads out of order).

I have tried sorteing my bam file but still no SNPs.

ADD REPLYlink written 4.2 years ago by Aiah10

I think you need to address this first error in order to get output. First, can you check your BAM file isn't truncated using SAMtools view? Second, can you make sure that the FASTA file you're using for variant calling is exactly the same you used for alignment?

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by Matt Miossec350
1

Thank you all I have checked BAM file isn't truncated and Yes I am using the same ref fasta file after indexing and tried it without .

ADD REPLYlink written 4.2 years ago by Aiah10

Is your problem solved?

ADD REPLYlink written 4.2 years ago by H.Hasani920
0
gravatar for Aiah
4.2 years ago by
Aiah10
Aiah10 wrote:

Ok so I have been using different versions of samtools and bcftools . Dlownloaded the comaptable ones . had to update the filters in the commands : samtools mpileup -uf A1163.fa P1.bam | bcftools call -c -v --output-type b | bcftools view -m 2 --output-type b > P1.raw.bcf and .

bcftools view P1.raw.bcf | perl vcfutils.pl varFilter -D2 > r1P.flt.vcf

still empty vcf file. Any suggestions ?

ADD COMMENTlink modified 4.2 years ago • written 4.2 years ago by Aiah10

According to bcftools manual call is intended to call variants, therefore, I'm expecting the output to be vcf. If so, you should have a valid file right there and your problems start in the steps following it.

ADD REPLYlink written 4.2 years ago by H.Hasani920

Thank you H. Hasani I did try: amtools mpileup -uf A1163.fa P1.bam | bcftools call -c -v --output-type b | bcftools view -m 2 --output-type b > P1.vcf.

and the result was abit unusual see below: \8B\00\00\00\00\00\FF\00BC\00L\B5\98ݒ\DA6\C73\B9\A3\D7}\00 {\B3\99bbI6\E0N\E9 \BB,f\F6\83\85\CD&W\ED\A0\A9\BF\D6$\FBL}\C9J\96lQ\A6Mo\9Fst~\E7K\BE\B8\BD~\FD\F7ϯ^\9D\9D\ADXLWY\91\D1\BA\ED\826\FA\E9\ECl4\BEy\BC\9A\F6\FB\93\C1l\D6R\BE(X.X\96\F6\9B\838\F2.A r\C29]6[\E3\E1Ǿ\FF\BB\BC\93\93DdY̟\E4\BFJ \DB藍\E01\9B{\F2\B2 \F7$\97Y\92\90tٯ\BE\83$\97\B1ls\E0mW\00a\B7WL{Ny[AW\B4\A0\E9\82\F6Uȿ\BE}[I\E4\8B,l]\C6;\9CQ\D4 z\AD\98\A6k\B1\E9C?lx_ \BBV\D8\F5K!t ;\B5\E2N)Dnah\85\A1\B6\88\DD\C2\C0 ;\B0n!\B6®^L\E8\A2=aX ;n!\AC\85\C8\D71v\DDBO\A8]\F7\9CBY\A1IO\E4Z2H&\BFL\B8 \B6hP\A4}C7l\D9\DF\D8t\C3\C1F:L覃-,\F3Z*\DDx\B0Ń{:\EB\D0\CD[>A\84\B5\D2 [@!4\DE݄\B0%V+r#BQ\E9z\83nF\C82\EA\98\EEAnF\C82\EA\9A\D2DnF\C82\EA!\BD"\E4f\84jF(4YBnF\A8f\84\BA\91\CE<r3B5#\EC\FBƦ\9B\AAa\EAZBnFZeO7r3B5#@\A7\9B\8C\AC\D2\D0DnF\B0f\84\C3H\D3\C4nF\B0f\84\A3P3\C2nF\B0c\95\86&v3\825\A3\00B3\E4܌\CD(\E8\EE\D8\CD֌B\DFt1v3\825\A30\EA\F5\B4\D2\CD֌z!6q\BAA\DBG>\EA\EA\D4c7$\DF6\ECDƪ\9B\92o;)\C4=3\E4ݘ|\DBJ\A1\FC\A6\A5nN\BE\ED\A5\AE\E8rܠ|\DBL\BD\A0 \C0Mʷ\A4z\BAN7+\BFf\D5 ;]l\B4nZ\BE\9Dza\A2P\E76p\F3\ED\E0\EB\81oz p3\F3-\B3u\83* '\B50\AA\A9!\E4G\9D\C8Xvr\93eU\8Be]\FA\E6\89:Ʌ\91m\B0ne\EA\B4\D8\C9N\A2\AD\B3\8C\E4\A3C\D6R)v\D2 #\DBf]\F9%е:\F9\85\D1?\95 ;,nK\E5\C7íܔ\E6\E54\908\A61=\E7o@&6\B4\00bCR\90\CD9-vt\D9n#\E3\BB\D1}ie|7\BC\BAi\DDm\939-\FA~\EB\F1%\A7\FDQLև\C6\C7\E9\92-\88\A0\D9\F2\83\82)I\D2_ J3\ED\A6\DDtYç\CA<\D4\E6ǩ\A0kZz\B8%\9FY\B2M@ZJA\B6%K\F86ϳB\B0t\ADܰtI\E3\A6ݱY/\B7\A3#/\A38#\C2\EDcU\90\85\FA\E5/\F8\C0\CBpr\CAR\A6\E4Si,i.6M\BB\AC\ED< /\BE\ED\93I\F0\90qA\E4^\0\DCߛ\BB \96\E71[P\8F3A)]Ʌq?\98\DE <N\9F\C1\92\CE\E7l-\A3T\A8\E6T\C8[\DF4[զ\BE\89\9BvZ8\9D\9C\E0-IS\EFÆ\89\94\BE\80\F7@V\87P\9D\AA\A5O2\93Jq\B4+\80zG[;\BD}\F8\F7NoI\9E\AB\84<lI\CC\C4˷\DDv\DC^\FC\B7\84\D3\D3|\F6\8E\96:\FBqk\DDq0\85<~};\82\E8 \82\D9\F5 \CC躠kR\82\9C˥.ABE\C1U\A3\FBG\8B\F2Oh\C2\D1^\F3\C9;L\9E\F3DM\AD/\A3\D6DŽ\D1\FD\F4v\A0\C7ޤ\9EV\D7\DF\E8\C1\A6s5\D9t\E9\F1\89e\F0k\9AfB\DEb\F6\8D\D9&˖\BC\B9w\D8\F3r\FDx\B4\99a\99\F3C'\D7\C6^s\EF\A0Pgc0\82'M$\CF$&O\EFTE\AEf\EC\8Arr\E0\9B\91.'}\DEʃ\F3 8\97\A7\F5m\A2\CA\E0݇\ABW\9D\A3\D0\8A@\D5\A7\F2)%sXd\F2XzH\87\BD6\B8\FC\92\B2ȶrL\9E\A7\99\F2 \CA JKUǭw\D2 \CF2\E1%\94\A4ޒ\82\82\C44ݳi:\CC"\DB\E9)\96ps\EFV;;s,\B7,Q\90ٜ\CCYe[. p\92\E41U\FD\A0\9C\A8\A5\CB_L\B1\A1\C3֛<\95\A3\E0뎼\89\B7T?C\B8\9Es91Zeo?\98k\B9P}YR\84\C5`Y=\81\D4\CFͽ\93a\ED\FFW\EE\F1ױ\DE\A0\AC\9B\B1* F+ۇ\8Fv \F5\FBI\BCܰ?\90\99Q\AA \CAYil~\96#\E7\DA<\EC\D0\F1c\FD8{Κ\B8\AB\B7&\B6\DExU%t\E5Ɍ~"\C5\B4\CAo\95\85\C1iKR\F5_*\9B\EA\BB\F9\B3 \AAZp9}\E6\8BU\F9\EE\ECO9\ACb\FB\CE \B7\A1}\EB\86˷n\CA\EA՛\BA\DEx;\E0y\D9V\E4[\E1\95\F9\9D\EF߱c\F4\D3i\B6\95\B2\B2\AD\AE\81\97\00\F4\A5\E9\CBw\D3\FB\DB\C6\E4~\D6ӫQC6f\E3\E1\FDড_<6T\96z\B26\CC \C0W\FF\00\CB\CEwG\B4\00\00\8B\00\00\00\00\00\FF\00BC\00\00\00\00\00\00\00\00\00\00\00 .

Any suggestions ?

ADD REPLYlink written 4.2 years ago by Aiah10

Just try

samtools mpileup -uf A1163.fa P1.bam | bcftools call -c -v --output-type b > output

ADD REPLYlink written 4.2 years ago by H.Hasani920

Thank you I tiered aligning with botwie2 instead of Bwa and it worked smoothly !!! not sure why but its working

ADD REPLYlink written 4.2 years ago by Aiah10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 690 users visited in the last hour