Variant calling from comparison of assembled genomes where there is only one read depth.
3
0
Entering edit mode
29 days ago
Hans ▴ 120

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

variant • 348 views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode
29 days ago

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.

ADD COMMENT
0
Entering edit mode
29 days ago
swbarnes2 9.9k

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

ADD COMMENT
0
Entering edit mode
28 days ago
Hans ▴ 120

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.

ADD COMMENT

Login before adding your answer.

Traffic: 1774 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6