Variant calling from comparison of assembled genomes where there is only one read depth.

Question

[Hans]

Hello I have generated two sorted bam files using minimap2 with the flag -ax asm5. Each of the files is a comparison of an assembled genome from a specific individual to a reference genome. I used bcftools mpileup and bcftools call to call variants:

bcftools mpileup  -Ou   -f ref.fasta ind1.bam ind2.bam  | bcftools call -Ov    -m  >var.vcf

In the vcf file I get GT 0/1 0/1 on positions where I see in IGV browser that one of the genotypes have SNP while the other does not. I guess that bcftools pipeline in not designed to call variants with 1 read depth. Is there a way to call variants with 1 read depth using bcftools? Is there a better alternative? I do not have the raw reads from which these individual assemblies were generated. Thank you

[samuel.a.odonnell]

To call variants from whole genome alignment you may want to look at other tools like minimap2-paftools or MUMmer-dnadiff

Answer 1

[Istvan Albert]

When ran your snp caller various cutoffs are applied, for example a p-value cutoff etc.

Run the tool with help flag and see what these parameters are. You may be able to turn off all cutoffs and generate all calls.

Since what you are trying to do is a bit "against" how SNP calling is supposed to operate it might not quite work out. In that case, you could try artificially increasing the reads ( basically duplicate/triplicate the reads you already have to "coerce" the algorithm in believing you have more coverage). Obviously this would be quite an artificial construct, thus keep in mind it might introduce other problems.

Answer 2

[swbarnes2]

If only one read shows a different letter, it's almost certainly a sequencing error, and not genuine.

Answer 3

[Hans]

It appears that bcftools takes into consideration the alignments of all genotypes when it calls variants. When I call the variants separately for each gnotype I get the correct results. I can then merge the bcf files using bcftools merge.