This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Having Issues with Variant Calling

Hi! I am a beginner in bioinformatics so I apologize if I am not too clear on the terminology.
I am trying to variant call some transcriptome data. So far I have:

  • mapped the reads to a refrence genome using Bowtie2
  • converted the .sam files into .bam files using samtools view
  • sorted and indexed the bamfiles using samtools

Next I know I will use the mpileup function but whenever I specify a region, the resulting vcf only gives me this when I open in excel: contig= ID=NW_019351050.1,length=6341

When I do not specify a region, I still get that same contig line but I get other regions that look like this in excel. 

NW_019316579.1  613728  .   C   <*> 0   .   DP=1;I16=0,1,0,0,40,1600,0,0,42,1764,0,0,25,625,0,0;QS=1,0;MQ0F=0   PL  0,3,40

I am a bit confused on what I am actually viewing. Any help would be appreciated! Thank you in advance

The code I used for mpileup: samtools mpileup -v -r NW_019351050.1 -f genomic.fna sorted.bam > variant.vcf.gz

sequencing snp samtools variant calling
• 2,604 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by shola_aleru • shola_aleru

what do you open in excel ? variant.vcf.gz ? the compressed file variant.vcf.gz ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I gunzip the vcf.gz file and then I view the vcf in excel.

• 0 views • link
written by shola_aleru • shola_aleru

1 answer

samtools mpileup produces a raw vcf file that must be called with "bcftools call"

see http://www.htslib.org/workflow/

to convert your BAM file into genomic positions we first use mpileup to produce a BCF file that contains all of the locations in the genome. We use this information to call genotypes and reduce our list of sites to those found to be variant by passing this file into bcftools call.

You can do this using a pipe as shown here:

bcftools mpileup -Ou -f <ref.fa> <sample1.bam> <sample2.bam> <sample3.bam> | bcftools call -vmO z -o <study.vcf.gz>
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

What is the reason why to never use excel? I don't know any other way to see what I actually have.

• 0 views • link
written by shola_aleru • shola_aleru

What is the reason why to never use excel?

excel introduces error in gene names, excel prevents you from having good practices with the command line

I don't know any other way to see what I actually have. 

gunzip -c  variant.vcf.gz   | more
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Our less can read compressed files without needing a gunzip

• 0 views • link
written by WouterDeCoster WouterDeCoster

Thank you!

Maybe I am not understanding it well but I followed that pipeline to get the vcf.gz file and then used the gunzip -c command, sent that to a file and the region I am interested in still looks like this :

contig= ID=NW_019351050.1,length=6341

• 0 views • link
written by shola_aleru • shola_aleru

Hello,

it is better to post every single command you use and also the first line of your resulting file, because "I followed that pipeline" and "looks like this" let so much space for a misinterpretation.

fin swimmer

• 0 views • link
written by finswimmer finswimmer

Log in to answer this question.

More posts like this