Body: <p>Maybe I'm having a brain-freeze so excuse me if this is a waste of everyone's time...</p>
  <p>I am working with an organism with ~700Mb genome.</p>
  <p>I have bam files of two individuals that are whole genome sequenced with the average depth of ~150.</p>
  <p>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 '<code>gatk or mpileup... argument that only retains sites when all are homozygote reference... delete tmp file...'</code>. So the question is: what would be the &quot;cheapest&quot;/fastest way to count where these individuals are homozygote reference?</p>
  <p>Come to think of it... I'm asking for a line count value on an &quot;inverted variant file&quot;/&quot;gVCF minus variation&quot;</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-09-14T07:00:06.243Z
Comments.Text: 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)
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-09-14T09:20:35.643Z
Comments.Text: 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.
Comments.UserDisplayName: AWE
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-09-14T09:44:05.957Z
Comments.Text: Like `bcftools call --skip-variants`
Comments.UserDisplayName: AWE
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-09-14T09:59:37.257Z
Comments.Text: 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,
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-09-14T17:02:11.230Z
Comments.Text: 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.
Comments.UserDisplayName: AWE
Tags: <sam><variant-calling><reference-genome><wgs><multiple-sequence-alignment>
Title: Fastest way to count where two bam files are homozygote reference (inverted variant calling)
ThreadText: # Fastest way to count where two bam files are homozygote reference (inverted variant calling)
  
  (Asked by: username_0 on 2021-09-13)
  
  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"
  
  **Tags:** sam, variant-calling, reference-genome, wgs, multiple-sequence-alignment
  
  > **Comments:**
  > **username_1**: 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)  
  > **username_0**: 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.  
  > **username_0**: Like `bcftools call --skip-variants`  
  > **username_1**: 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,  
  > **username_0**: 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.  
  
  ---
  
  ## 0 answers
  
  