SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
SAMtools missing bcftools khadeeja Bioinformatics 8 06-19-2013 02:06 AM
Samtools mpileup/bcftools cristae8 Bioinformatics 3 05-02-2012 12:43 PM
Samtools/bcftools problem abelhj Bioinformatics 5 04-26-2011 06:02 AM
variant calling using samtools -v- bcftools ksc Bioinformatics 2 04-13-2011 06:44 AM
samtools/bcftools problem? nimmi Bioinformatics 0 01-12-2011 11:59 AM

Reply
 
Thread Tools
Old 12-22-2011, 05:14 AM   #1
alma
Junior Member
 
Location: UK

Join Date: Jun 2011
Posts: 2
Unhappy Problem with samtools/bcftools

Hi

I am trying to call SNPs and indels from some bam files and I tried the following commands as given in the samtools/bcftools manual:

samtools mpileup -uf reference.fa aln1.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

However after running the commands I see that var.raw.bcf and var.flt.vcf files have results from only the first chromsome that is present in the bam file. I also get the following message:

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
[afs] 0:16349.158 1:6949.156 2:48149.686

In my bam file Chr10 appears first followed by Chr11, 12, ...19, 1, 20, 21..., 2 etc. The reads, however, are sorted within the chromosomes. Just to make sure that the problem is not due to such arrangement of the chromsomes, I extracted Chr1, 2 and 3 and created a new bam file and sorted it using following commands.

samtools view -bh aln1.bam Chr1 Chr2 Chr3 -o Chr1to3.bam
samtools sort Chr1to3.bam Chr1to3sorted

Then ran the same command as before for SNP/Indel analysis. But the problem exists.

Can any of you please suggest what might be the possible problems/solutions?

Thanks in advance

ALM
alma is offline   Reply With Quote
Old 12-22-2011, 10:28 AM   #2
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Try running it without chromosome 1 and see what happens. Make sure the headers for all files used are complete and the reference sequence file is complete as well.
Heisman is offline   Reply With Quote
Old 01-05-2012, 03:11 AM   #3
amruta.bn
Junior Member
 
Location: India

Join Date: Jan 2012
Posts: 6
Default

Hi,

I have generated a VCF file using the command

samtools faidx input.reference.fa
samtools index sortedhg19.bam
samtools mpileup -d 8000 -uf input.reference.fa my.sortedhg19.bam | bcftools view -bvcg -> my.raw.bcf
bcftools view my.raw.bcf | perl vcfutils.pl varFilter -D8000 > my.flt.vcf

But when i validated my VCF file using vcf-validator of VCFTools and i am getting the warning message as given below

The header tag 'reference' not present. (Not required but highly recommended.)
The header tag 'contig' not present for CHROM=chr1. (Not required but highly recommended.)
INFO field at chr1:10007 .. Could not validate the float [-nan]

Can any one please help me to fix this.I think only the third one is critical.

Regards,
Amruta
amruta.bn is offline   Reply With Quote
Old 01-05-2012, 03:48 AM   #4
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Could you post that entire line in the file?
Heisman is offline   Reply With Quote
Old 01-06-2012, 12:40 AM   #5
amruta.bn
Junior Member
 
Location: India

Join Date: Jan 2012
Posts: 6
Default

##fileformat=VCFv4.1
##samtoolsVersion=0.1.18 (r982:295)
##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=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=VDB,Number=1,Type=Float,Description="Variant Distance Bias">
##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=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">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 29dec/sorted.29dec.hg19.bam.bam
chr1 10007 . T C 6.98 . DP=127;VDB=-nan;AF1=0.5001;AC1=1;DP4=1,3,0,7;MQ=20;FQ=5.64;PV4=0.36,1,1,1 GT:PL:GQ 0/1:36,0,33:34
chr1 10186 . T A 6.98 . DP=140;VDB=-nan;AF1=0.4999;AC1=1;DP4=2,8,5,7;MQ=20;FQ=9.53;PV4=0.38,1,1,1 GT:PL:GQ 0/1:36,0,76:37
chr1 10549 . T A 9.52 . DP=14;VDB=-nan;AF1=0.5;AC1=1;DP4=2,1,1,2;MQ=20;FQ=9.52;PV4=1,1,1,1 GT:PL:GQ 0/1:39,0,39:39
chr1 12613 . G C 6.02 . DP=6;VDB=-nan;AF1=1;AC1=2;DP4=0,0,0,2;MQ=20;FQ=-33 GT:PL:GQ 1/1:36,6,0:6
chr1 13275 . G A,C 3.41 . DP=9;VDB=-nan;AF1=1;AC1=2;DP4=0,0,1,2;MQ=20;FQ=-33 GT:PL:GQ 1/1:43,17,11,31,0,28:4
chr1 13407 . T A 6.02 . DP=10;VDB=-nan;AF1=1;AC1=2;DP4=0,0,0,2;MQ=20;FQ=-33 GT:PL:GQ 1/1:36,6,0:6
chr1 13411 . C A 5.29 . DP=10;VDB=-nan;AF1=1;AC1=2;DP4=0,0,1,2;MQ=20;FQ=-33 GT:PL:GQ 1/1:35,6,0:6
chr1 16259 . T A 6.02 . DP=8;VDB=-nan;AF1=1;AC1=2;DP4=0,0,2,0;MQ=20;FQ=-33 GT:PL:GQ 1/1:36,6,0:6
chr1 17216 . T G 3.64 . DP=6;VDB=-nan;AF1=0.5205;AC1=1;DP4=1,0,1,1;MQ=20;FQ=-16.1;PV4=1,0.03,1,0.44 GT:PL:GQ 0/1:31,0,11:18
chr1 18712 . T A 3.54 . DP=15;VDB=-nan;AF1=0.4998;AC1=1;DP4=3,1,2,3;MQ=20;FQ=5.35;PV4=0.52,0.042,1,1 GT:PL:GQ 0/1:31,0,45:30
chr1 20700 . C G 3.64 . DP=9;VDB=-nan;AF1=0.5205;AC1=1;DP4=1,0,1,1;MQ=20;FQ=-16.1;PV4=1,1,1,0.28 GT:PL:GQ 0/1:31,0,11:18
chr1 20794 . C A 3.64 . DP=17;VDB=-nan;AF1=0.5205;AC1=1;DP4=1,0,1,1;MQ=20;FQ=-16.1;PV4=1,0.056,1,0.2 GT:PL:GQ 0/1:31,0,11:18
chr1 21113 . T A 3.64 . DP=9;VDB=-nan;AF1=0.5205;AC1=1;DP4=0,1,1,1;MQ=20;FQ=-16.1;PV4=1,0.33,1,1 GT:PL:GQ 0/1:31,0,11:18
chr1 21130 . C A 7.39 . DP=11;VDB=-nan;AF1=0.5718;AC1=1;DP4=0,1,2,3;MQ=20;FQ=-21;PV4=1,0.023,1,0.44 GT:PL:GQ 0/1:36,0,6:10
chr1 21186 . G A 6.02 . DP=9;VDB=-nan;AF1=1;AC1=2;DP4=0,0,1,1;MQ=20;FQ=-33 GT:PL:GQ 1/1:36,6,0:6
chr1 21253 . G A 8.64 . DP=24;VDB=-nan;AF1=0.5;AC1=1;DP4=5,2,3,1;MQ=20;FQ=11.3;PV4=1,0.1,1,0.1 GT:PL:GQ 0/1:38,0,76:40
chr1 21300 . C A 3.54 . DP=19;VDB=-nan;AF1=0.4998;AC1=1;DP4=4,2,1,2;MQ=20;FQ=5.47;PV4=0.52,0.34,1,1 GT:PL:GQ 0/1:31,0,70:30
chr1 21415 . C A 5.46 . DP=16;VDB=-nan;AF1=0.4999;AC1=1;DP4=3,2,1,2;MQ=20;FQ=7.77;PV4=1,0.41,1,0.26 GT:PL:GQ 0/1:34,0,55:34
chr1 21564 . C A 5.29 . DP=4;VDB=-nan;AF1=1;AC1=2;DP4=0,0,2,1;MQ=20;FQ=-33 GT:P



This is the VCF file.Hope you asked for vcf file .
Is there any way to format reference fasta sequence?
i have downloaded individual chromosomes from 1,2...22,X and y and joined into a single file in two orders like 1,11,12,13 ... and 1,2,3,4,.. .But in both cases this error is coming after validation
amruta.bn is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:39 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO