Question: Having Issues with Variant Calling
0
gravatar for shola_aleru
16 months ago by
shola_aleru10
shola_aleru10 wrote:

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

ADD COMMENTlink modified 16 months ago by Pierre Lindenbaum121k • written 16 months ago by shola_aleru10

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

ADD REPLYlink modified 16 months ago • written 16 months ago by Pierre Lindenbaum121k
1

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

ADD REPLYlink written 16 months ago by shola_aleru10

I almost forgot : don't use excel.

ADD REPLYlink written 16 months ago by Pierre Lindenbaum121k

and what is the version of samtools ?

ADD REPLYlink written 16 months ago by Pierre Lindenbaum121k

I am using samtools 1.5

ADD REPLYlink written 16 months ago by shola_aleru10
1
gravatar for Pierre Lindenbaum
16 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum121k wrote:

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>
ADD COMMENTlink written 16 months ago by Pierre Lindenbaum121k

and do not use excel. Never.

ADD REPLYlink written 16 months ago by Pierre Lindenbaum121k

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

ADD REPLYlink written 16 months ago by shola_aleru10

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
ADD REPLYlink written 16 months ago by Pierre Lindenbaum121k
1

Our less can read compressed files without needing a gunzip

ADD REPLYlink written 16 months ago by WouterDeCoster40k

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

ADD REPLYlink written 16 months ago by shola_aleru10
1

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

ADD REPLYlink written 16 months ago by finswimmer11k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1493 users visited in the last hour