Having Issues with Variant Calling

Question

[shola_aleru]

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

[Pierre Lindenbaum]

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

[shola_aleru]

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

[Pierre Lindenbaum]

I almost forgot : don't use excel.

[Pierre Lindenbaum]

and what is the version of samtools ?

[shola_aleru]

I am using samtools 1.5

Answer 1

[Pierre Lindenbaum]

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>

[Pierre Lindenbaum]

and do not use excel. Never.

[shola_aleru]

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

[Pierre Lindenbaum]

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

[WouterDeCoster]

Our less can read compressed files without needing a gunzip

[shola_aleru]

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

[finswimmer]

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