0
$\begingroup$

Maybe I'm having a brain-freeze so excuse me if this is a waste of everyone's time...

I am working with an organism with ~700Mb genome.

I have bam files of two individuals that are whole genome sequenced with the average depth of ~150.

I would like to count the occurrences when both of the individuals are homozygote for the reference allele. Preferably without any hassle or space, I have no need to create a vcf file for the task or save the information in any form. Something like 'gatk or mpileup... argument that only retains sites when all are homozygote reference... delete tmp file...'. So the question is: what would be the "cheapest"/fastest way to count where these individuals are homozygote reference?

Come to think of it... I'm asking for a line count value on an "inverted variant file"/"gVCF minus variation"

Improve this question
$\endgroup$
5
  • $\begingroup$ There is no genome reference for your organism ? What you want to do is very simple, you want to do variant calling using the BAM files (there are a number of ways for that) and then count the sites with variation, then you just do (number_of_all_sites_in_genome minus number_of_variant_sites) $\endgroup$
    – JRodrigoF
    Sep 14, 2021 at 7:00
  • $\begingroup$ I was thinking of all "callable" sites according to my bams (that meet some filtering criteria). Not just reference size minus variant sites. Maybe I will just do some genotype file and parse through it. I just thought I was forgetting some obvious command in the known tools that would be fast and easy. $\endgroup$
    – AWE
    Sep 14, 2021 at 9:20
  • $\begingroup$ Like bcftools call --skip-variants $\endgroup$
    – AWE
    Sep 14, 2021 at 9:44
  • $\begingroup$ gatk has an option to output all sites and not only variants, you could do your calculation on that VCF and you will be able to filter reference sites based on reference allele quality, on the other hand you will have to deal with a massive VCF and slow computation, as opposed to the size minus #_variants approach, $\endgroup$
    – JRodrigoF
    Sep 14, 2021 at 9:59
  • $\begingroup$ For now it's gvcf via bcftools and related filtering and then something like awk if substr$x,1 == substr$x,3 and $x==$y.... Heavy duty for minor reward, the antonym of judo. $\endgroup$
    – AWE
    Sep 14, 2021 at 17:02

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.