<p>I have a little function that can report this by using samtools to look at the read pileup. Um. It isn't the most elegant piece of code, but it does sorta do the job:</p>
<pre class="lang-sh prettyprint-override"><code>variantSupport () 
{ 
  genome=$1
  bam=$2;
  range=$3;
  samtools mpileup -f &quot;$genome&quot; -r &quot;$range&quot; &quot;$bam&quot; | cut -f 5 |
    tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c |
    awk 'BEGIN{tot=1}$2~/,/{ $2=&quot;.&quot; }
         {
           tot+=$1;
           c[$2]+=$1
         }
         $2~/[a-zA-Z]/{ varReads+=$1; }
         END{
           for(i in c){
             print i,c[i]
           }
           printf &quot;%s of %s reads (%.1f%%), support the variant.\n&quot;,varReads,tot,varReads*100/tot
   }'
}
</code></pre>
<p>Add the function to your <code>~/.bashrc</code> or just paste it directly into a terminal and you can then do:</p>
<pre><code>$ variantSupport /path/to/genome/hg19.fa foo.bam chr13:28607916-28607916
[mpileup] 1 samples in 1 input files
A 915
* 1
. 93
915 of 1010 reads (90.6%), support the variant.
</code></pre>
<p>This is telling me that at position 28607916 of chr13, I have 1 read supporting a deletion (<code>*</code>, see <a href="https://www.htslib.org/doc/samtools-mpileup.html" rel="nofollow noreferrer">https://www.htslib.org/doc/samtools-mpileup.html</a> for an explanation of the symbols), 93 reads supporting the reference residue (T, in this case) and 915 supporting an A, so a T=&gt; A variant.</p>
<p>However, a better approach is to use a variant caller and force it to call the variant you want. You can then investigate the support for it. With <code>mutect2</code>, you can use the <code>--alleles</code> option to force the calling and then use the <code>-L</code> option to limit calling to only the specific region. If you give the same vcf file for both options, that will return results for the variants you want and only those. For example, given this file <code>foo.vcf.gz</code> (which needs to be compressed using <code>bgzip</code> and indexed using <code>tabix</code>):</p>
<pre><code>$ zcat 1.vcf.gz 
##fileformat=VCFv4.3
##reference=/genomes/hg19/hg19.fa
##FORMAT=&lt;ID=AB,Number=A,Type=Float,Description=&quot;Allelic Balance&quot;&gt;
##FORMAT=&lt;ID=DP,Number=1,Type=Integer,Description=&quot;Read depth&quot;&gt;
##FORMAT=&lt;ID=GQ,Number=1,Type=Integer,Description=&quot;Net Genotype quality across all datasets, calculated from GQ scores of callsets supporting the consensus GT, using only one callset from each dataset&quot;&gt;
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description=&quot;Genotype&quot;&gt;
##contig=&lt;ID=chr13,assembly=hg19,length=115169878&gt;
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG003
chr13   28607917    .   G   C   .   .   .   GT:AB:DP:GQ 1/1:0.6:395:318
</code></pre>
<p>You can then do:</p>
<pre><code>gatk Mutect2 \
  -R /genomes/hg19/hg19.fa \
  --alleles 1.vcf.gz  \
  -L 1.vcf.gz \
  -I foo.bam \
  -O out.vcf
</code></pre>
<p>This will produce an <code>out.vcf</code> like this (I am excluding most headers for clarity):</p>
<pre><code>
##FORMAT=&lt;ID=AD,Number=R,Type=Integer,Description=&quot;Allelic depths for the ref and alt alleles in the order listed&quot;&gt;
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  331970
chr13   28607917    .   G   C   .   .   DP=1107;ECNT=1;MBQ=30,0;MFRL=353,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-2.716e+00 GT:AD:AF:DP:F1R2:F2R1:SB    0/1:1043,0:9.616e-04:1043:518,0:496,0:989,54,0,0
</code></pre>
<p>I can see this is not a real variant since the value for AD, the number of reads supporting the reference and the variant is <code>1043,0</code>, meaning that 1043 reads support the reference and none support the variant.</p>
<p>This way, you can get results for your variants of interest and see if there is any evidence supporting them.</p>

CC BY-SA 4.0
2024-03-04T14:47:53.423Z
Thanks! Your answer does what I needed. However, I have noticed that the Allele Depth numbers in the VCF (out.vcf, in your code) are lower than expected. For instance, I got REF=105 and ALT=104, whereas the numbers are  REF = 436 and ALT = 380 if I extract this information with bcftools and set -d to a high number. I need the total read numbers, so how can I edit your command to achieve this?
Nereus
CC BY-SA 4.0
2024-03-04T14:53:34.077Z
@Nereus I don't think you need the total, actually. The ones that are not counted by my function are the the ones  samtools skips by default (see https://www.htslib.org/doc/samtools-mpileup.html; for example, it won't count bases with a fastq quality under 13), and the ones not included in the `AD` value are those the variant caller filtered out. You can't just use all reads, you _must_ only take into account the ones that pass some basic QC filtering otherwise, you just get nonsense. If you see such huge differences, I would investigate the quality (both mapping and fastq) in those regions.
terdon
CC BY-SA 4.0
2024-03-04T15:18:30.087Z
Sounds good, thank you for your advice
Nereus
CC BY-SA 4.0
2024-03-06T14:59:29.750Z
may I ask how you produced the 1.vcf.gz file in your example?
Nereus
CC BY-SA 4.0
2024-03-06T15:02:06.697Z
@Nereus from one I had. I just looked at a bam file I had lying around and made the vcf point to a variant in that file and ran the command I showed you.
terdon
CC BY-SA 4.0
2024-03-06T16:50:31.330Z
I see. Would it also be possible to create a custom VCF if I haven't got one available, e.g. with data from Ensembl? Thank you for your time!
Nereus
CC BY-SA 4.0
2024-03-06T17:07:52.217Z
@Nereus um, sure, why not? I mean, a vcf is just a text file, you can make one manually if you want.
terdon
CC BY-SA 4.0
2024-03-08T16:27:01.847Z
I have noticed that your approach produces the desired result if the variant is actually present; however, when the result is not present an empty VCF file is created (I'm adding an Edit above)
Nereus
CC BY-SA 4.0
2024-03-03T17:45:10.007Z
2024-03-03T17:45:10.007Z
terdon
<p>I'm working with targeted Illumina sequencing data generated with DNA from diseased and healthy tissue (this is an age-related disease and is not cancer/neoplastic). My hypothesis is that the diseased tissue might contain specific somatic mutations in one particular gene (let's call it GENE), while the healthy tissue will not. The reasoning behind this: The same mutations cause autosomal-dominant versions of my disease of interest when present in the germline (my diseased tissue is from sporadic disease). My coverage varies, but is in the range of 1000x - 50'000x. My sequencing protocol employed UMIs.</p>
<p>I have run Mutect2 in default settings and with lower thresholds (<code>--tumor-lod-to-emit</code>) and haven't found any calls that correspond to these pathogenic mutations, so my results are negative. Nevertheless, I have reason to believe that my pathogenic mutations might be present at <strong>very</strong> low VAF (around 0.05-0.1 %), based on other experiments.</p>
<p><strong>My question</strong>: How do I quantify the number of reads that correspond to a particular variant of interest? Let's say that base 2000 in GENE is usually an A, whereas the mutation A&gt;T would be pathogenic. How do I count the number of Ts in position 2000?</p>
<p>As I'm somewhat new to bioinformatics, I'm also very open to any other ideas. Thanks!</p>
<hr />
<p><strong>Edit</strong>: @terdon's approach worked perfectly in cases where I used a VCF for <code>--alleles</code> that had been generated as an output of Mutect2, such as the following one:</p>
<p><em>VCF produced by Mutect</em>:</p>
<pre><code>#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847
</code></pre>
<p>Curiously, if I used a manually produced file such as the following, no output was produced, unless the variant of interest was heterozygous:</p>
<p><em>Manually produced VCF</em>:</p>
<pre><code>#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 rs1799990   A   G   .   .   .
</code></pre>
<p><em>Output using the manual VCF, heterozygous sample</em>:</p>
<pre><code>#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  heterogyzous_sample
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847
</code></pre>
<p><em>Output using the manual VCF, non-heterozygous sample</em>:</p>
<pre><code>#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  novariant_sample
</code></pre>
<p>Manual inspection of this non-heterozygous sample revealed that it contained 1 read with the ALT variant, so this doesn't seem to depend on the variant being absent.</p>

<variant-calling><illumina><gatk>
How do I quantify a specific somatic variant?
# How do I quantify a specific somatic variant?

(Asked by: username_0 on 2024-03-03)

I'm working with targeted Illumina sequencing data generated with DNA from diseased and healthy tissue (this is an age-related disease and is not cancer/neoplastic). My hypothesis is that the diseased tissue might contain specific somatic mutations in one particular gene (let's call it GENE), while the healthy tissue will not. The reasoning behind this: The same mutations cause autosomal-dominant versions of my disease of interest when present in the germline (my diseased tissue is from sporadic disease). My coverage varies, but is in the range of 1000x - 50'000x. My sequencing protocol employed UMIs.
I have run Mutect2 in default settings and with lower thresholds (--tumor-lod-to-emit) and haven't found any calls that correspond to these pathogenic mutations, so my results are negative. Nevertheless, I have reason to believe that my pathogenic mutations might be present at very low VAF (around 0.05-0.1 %), based on other experiments.
My question: How do I quantify the number of reads that correspond to a particular variant of interest? Let's say that base 2000 in GENE is usually an A, whereas the mutation A>T would be pathogenic. How do I count the number of Ts in position 2000?
As I'm somewhat new to bioinformatics, I'm also very open to any other ideas. Thanks!

Edit: @terdon's approach worked perfectly in cases where I used a VCF for --alleles that had been generated as an output of Mutect2, such as the following one:
VCF produced by Mutect:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Curiously, if I used a manually produced file such as the following, no output was produced, unless the variant of interest was heterozygous:
Manually produced VCF:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 rs1799990   A   G   .   .   .

Output using the manual VCF, heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  heterogyzous_sample
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Output using the manual VCF, non-heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  novariant_sample

Manual inspection of this non-heterozygous sample revealed that it contained 1 read with the ALT variant, so this doesn't seem to depend on the variant being absent.

**Tags:** variant-calling, illumina, gatk

---

## 1 answers

### Answer by username_1 on 2024-03-03. Score: 4 | Accepted answer

I have a little function that can report this by using samtools to look at the read pileup. Um. It isn't the most elegant piece of code, but it does sorta do the job:
variantSupport () 
{ 
  genome=$1
  bam=$2;
  range=$3;
  samtools mpileup -f "$genome" -r "$range" "$bam" | cut -f 5 |
    tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c |
    awk 'BEGIN{tot=1}$2~/,/{ $2="." }
         {
           tot+=$1;
           c[$2]+=$1
         }
         $2~/[a-zA-Z]/{ varReads+=$1; }
         END{
           for(i in c){
             print i,c[i]
           }
           printf "%s of %s reads (%.1f%%), support the variant.\n",varReads,tot,varReads*100/tot
   }'
}

Add the function to your ~/.bashrc or just paste it directly into a terminal and you can then do:
$ variantSupport /path/to/genome/hg19.fa foo.bam chr13:28607916-28607916
[mpileup] 1 samples in 1 input files
A 915
* 1
. 93
915 of 1010 reads (90.6%), support the variant.

This is telling me that at position 28607916 of chr13, I have 1 read supporting a deletion (*, see https://www.htslib.org/doc/samtools-mpileup.html for an explanation of the symbols), 93 reads supporting the reference residue (T, in this case) and 915 supporting an A, so a T=> A variant.
However, a better approach is to use a variant caller and force it to call the variant you want. You can then investigate the support for it. With mutect2, you can use the --alleles option to force the calling and then use the -L option to limit calling to only the specific region. If you give the same vcf file for both options, that will return results for the variants you want and only those. For example, given this file foo.vcf.gz (which needs to be compressed using bgzip and indexed using tabix):
$ zcat 1.vcf.gz 
##fileformat=VCFv4.3
##reference=/genomes/hg19/hg19.fa
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##contig=
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG003
chr13   28607917    .   G   C   .   .   .   GT:AB:DP:GQ 1/1:0.6:395:318

You can then do:
gatk Mutect2 \
  -R /genomes/hg19/hg19.fa \
  --alleles 1.vcf.gz  \
  -L 1.vcf.gz \
  -I foo.bam \
  -O out.vcf

This will produce an out.vcf like this (I am excluding most headers for clarity):

##FORMAT=
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  331970
chr13   28607917    .   G   C   .   .   DP=1107;ECNT=1;MBQ=30,0;MFRL=353,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-2.716e+00 GT:AD:AF:DP:F1R2:F2R1:SB    0/1:1043,0:9.616e-04:1043:518,0:496,0:989,54,0,0

I can see this is not a real variant since the value for AD, the number of reads supporting the reference and the variant is 1043,0, meaning that 1043 reads support the reference and none support the variant.
This way, you can get results for your variants of interest and see if there is any evidence supporting them.

> **Comments:**
> **username_0**: Thanks! Your answer does what I needed. However, I have noticed that the Allele Depth numbers in the VCF (out.vcf, in your code) are lower than expected. For instance, I got REF=105 and ALT=104, whereas the numbers are  REF = 436 and ALT = 380 if I extract this information with bcftools and set -d to a high number. I need the total read numbers, so how can I edit your command to achieve this?  
> **username_1**: @Nereus I don't think you need the total, actually. The ones that are not counted by my function are the the ones  samtools skips by default (see https://www.htslib.org/doc/samtools-mpileup.html; for example, it won't count bases with a fastq quality under 13), and the ones not included in the `AD` value are those the variant caller filtered out. You can't just use all reads, you _must_ only take into account the ones that pass some basic QC filtering otherwise, you just get nonsense. If you see such huge differences, I would investigate the quality (both mapping and fastq) in those regions.  
> **username_0**: Sounds good, thank you for your advice  
> **username_0**: may I ask how you produced the 1.vcf.gz file in your example?  
> **username_1**: @Nereus from one I had. I just looked at a bam file I had lying around and made the vcf point to a variant in that file and ran the command I showed you.  
> **username_0**: I see. Would it also be possible to create a custom VCF if I haven't got one available, e.g. with data from Ensembl? Thank you for your time!  
> **username_1**: @Nereus um, sure, why not? I mean, a vcf is just a text file, you can make one manually if you want.  
> **username_0**: I have noticed that your approach produces the desired result if the variant is actually present; however, when the result is not present an empty VCF file is created (I'm adding an Edit above)  

---

