SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bcftools consensus from BCF/VCF rwhet052 Bioinformatics 1 02-04-2015 06:03 AM
error when running samtools mpileup and bcftools linp5 Bioinformatics 5 08-22-2014 05:41 AM
Samtools output errors (ref is all 'N' in vcf) rajj Bioinformatics 4 12-03-2013 10:22 AM
samtools/bcftools: How to inspect SNPs not made it into vcf file muenalan Bioinformatics 2 09-06-2012 06:13 AM
Read Depth in vcf (samtools / bcftools) Marie_Noir Bioinformatics 1 04-17-2012 07:48 AM

Reply
 
Thread Tools
Old 06-28-2017, 06:29 AM   #1
MGCBrown
Member
 
Location: Canada

Join Date: May 2017
Posts: 15
Angry Very strange .vcf output after running samtools and bcftools

Hello,

After I run the samtools variant calling pipeline, I'm not getting the conventional .vcf output that I was expecting. Any advice or insight would be greatly appreciated.

Here is my command:

$ samtools mpileup -gf /scratch/$USER/Trinity.adj.fasta NFS_6525_adj_bowtie2.coordSorted.bam NFS_50254_adj_bowtie2.coordSorted.bam SFS_CC1_adj_bowtie2.coordSorted.bam SFS_25428_adj_bowtie2.coordSorted.bam | bcftools view -bvcg - > NFS_6525_50254_SFS_CC1_25428_adj_var.bcf

$ bcftools view NFS_6525_50254_SFS_CC1_25428_adj_var.bcf | vcfutils.pl varFilter -D 100 - > NFS_6525_50254_SFS_CC1_25428_adj_var.flt.vcf

I used VIM to view my NFS_6525_50254_SFS_CC1_25428_adj_var.flt.vcf file and it does not look like any examples that I've seen online. It's a combination of alpha-numeric characters and a list of my Trinity contigs. I've copied a sample of my .vcf file below:

(there is a full list of my Trinity contigs)
##contig=<ID=TRINITY_DN311005_c0_g1_i1,length=296,IDX=667121>
##contig=<ID=TRINITY_DN311035_c0_g1_i1,length=368,IDX=667122>
##contig=<ID=TRINITY_DN311050_c0_g1_i1,length=302,IDX=667123>
##contig=<ID=TRINITY_DN311084_c0_g1_i1,length=319,IDX=667124>
##contig=<ID=TRINITY_DN311076_c0_g1_i1,length=250,IDX=667125>
##contig=<ID=TRINITY_DN311012_c0_g1_i1,length=316,IDX=667126>
##contig=<ID=TRINITY_DN311030_c0_g1_i1,length=269,IDX=667127>
##contig=<ID=TRINITY_DN311065_c0_g1_i1,length=201,IDX=667128>
##contig=<ID=TRINITY_DN311069_c0_g1_i1,length=285,IDX=667129>
##contig=<ID=TRINITY_DN311034_c0_g1_i1,length=283,IDX=667130>

(some of the alphanumeric characters were changed to symbols when I copied and pasted here)
    T7<*>   @ B PE B rE HB @D  % @   1 , {  .
H     G7<*>   @ B E B rE HB @D  % @   1 , "{  .
I     G7<*>   @ B E B rE HB @D  % @   1 , &{  .
J     G7<*>   @ B E B rE HB @D  % @   1 , &{  .
K     G7<*>   @ B E B rE HB @D  % @   1 , &{  .
L     G7<*>   @ B E B rE HB @D  % @   1 , &{  .
M     G7<*>   @ B E B rE HB @D  % @   1 , &{  .
N     A7<*>   @ B E B rE HB @D  % @   1 , ${  .
O     G7<*>   @ B E B rE HB @D  % @   1 , #{  .
P     T7<*>   @ B E B rE HB @D  % @   1 , #{  .
Q     C7<*>   @ B E B rE HB @D  % @   1 , &{  .
R     T7<*>   @ B @E B rE HB @D  % @   1 , &{  .
S     A7<*>   @ B E B rE HB

However, my output looks nothing like this, which is an example I found online:

##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4
20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
20 1234567 microsat1 GTCT G,GTACT 50 PASS NS=3;DP=9;AA=G GT:GQP 0/1:35:4 0/2:17:2 1/1:40:3



Thank you.
MGCBrown is offline   Reply With Quote
Old 06-28-2017, 11:51 AM   #2
SamCurt
Member
 
Location: Iowa

Join Date: May 2010
Posts: 40
Default

Hi,

Which version of bcftools are you using? bcftools view are two different programs before and after 1.0. My guess is the -b argument in bcftools that's the problem, since it makes a BCF output. Try using bcftools without that argument, i.e. bcftools view -vcg [...]
SamCurt is offline   Reply With Quote
Old 06-28-2017, 12:00 PM   #3
MGCBrown
Member
 
Location: Canada

Join Date: May 2017
Posts: 15
Default

Thanks for the reply, I appreciate it. I'm using samtools and bcftools -1.4.1. I'll try running it now. Thanks again.
MGCBrown is offline   Reply With Quote
Old 06-28-2017, 12:13 PM   #4
SamCurt
Member
 
Location: Iowa

Join Date: May 2010
Posts: 40
Default

In that case, you used the wrong program. Before bcftools 1.0 the output of mpileup is called with bcftools view, but after that version the correct program should be bcftools call. So the command line should be:

samtools mpileup -gf /scratch/$USER/Trinity.adj.fasta NFS_6525_adj_bowtie2.coordSorted.bam NFS_50254_adj_bowtie2.coordSorted.bam SFS_CC1_adj_bowtie2.coordSorted.bam SFS_25428_adj_bowtie2.coordSorted.bam | bcftools call -vc -O b -o NFS_6525_50254_SFS_CC1_25428_adj_var.bcf

The definition of -g is also changed in latter versions to mean a GVCF output, and I have yet found a modern equivalent to -c.
SamCurt 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 05:31 AM.


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