VCF file without SNPs
1
0
Entering edit mode
8.4 years ago
Aiah ▴ 10

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 • 5.3k views
ADD COMMENT
1
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
1
Entering edit mode

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 REPLY
0
Entering edit mode

Is your problem solved?

ADD REPLY
0
Entering edit mode
8.4 years ago
Aiah ▴ 10

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 COMMENT
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

Just try

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 1452 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6