Read Counting¶

Pipeline step: count_reads (Step 2)

Module: grid.utils.count_reads

Step 2 counts the number of reads that originate within the VNTR region of interest for each sample. The resulting per-sample counts serve as a locus-specific proxy for sequencing depth and feed directly into the diploid copy number calculation in Step 4 (see Diploid Copy Number Estimation).

Read Filter¶

For each sample, GRiD opens the CRAM/BAM file with pysam and fetches all reads whose alignment overlaps the region \([s,\, e)\) on chromosome \(c\). A read \(r\) is counted if and only if it passes every one of the following conditions simultaneously:

\[\text{flag}(r) \in \mathcal{F} \quad \wedge \quad \text{MAPQ}(r) \geq q_{\min} \quad \wedge \quad \text{chr_mate}(r) = c \quad \wedge \quad \neg\,\text{dup}(r) \quad \wedge \quad \neg\,\text{sec}(r) \quad \wedge \quad s \leq \text{pos}(r) < e\]

where:

Symbol	Meaning
\(\mathcal{F}\)	User-specified set of SAM bitwise flags (`count_reads.flags` in config). Each flag must be an exact match — not a bitmask test.
\(q_{\min}\)	Minimum mapping quality (`count_reads.min_mapq`).
\(\text{chr_mate}(r) = c\)	Both reads in the pair align to the same chromosome (proper pair).
\(\neg\,\text{dup}(r)\)	Read is not marked as a PCR/optical duplicate.
\(\neg\,\text{sec}(r)\)	Read is not a secondary alignment.
\(s \leq \text{pos}(r) < e\)	The read’s start position falls within the target region boundaries. Reads that overlap the region but start outside it are excluded.

The per-sample read count is:

\[n_i = \bigl|\bigl\{\, r \in \mathcal{R}_i \;\mid\; \text{all conditions above}\,\bigr\}\bigr|\]

where \(\mathcal{R}_i\) is the set of reads fetched for sample \(i\).

Flag Selection¶

The flags list controls which read orientations are counted. The filter uses an exact flag comparison (not a bitwise AND), so each value in the list targets a specific orientation. The recommended defaults capture properly paired reads on both strands:

Flag	Orientation
`83`	Proper pair, read on reverse strand, mate on forward strand (R1 reverse)
`147`	Proper pair, mate on reverse strand, read on forward strand (R2 reverse)
`81`	Read on reverse strand, mate on forward strand (not requiring proper pair flag)
`145`	Mate on reverse strand (not requiring proper pair flag)

Flags can be added or removed to match the read orientation profile of a specific library preparation. Use the Broad Picard flag explainer to decode or construct flag values.

Implementation¶

Samples are processed in parallel using a ThreadPoolExecutor with the thread count set by threads in the config. Results are written to the output file incrementally as each sample completes, using a threading lock for safe concurrent writes.

The output is a tab-separated file with one row per sample:

HG00096.hg38.cram    1842
HG00097.hg38.cram    1956
NA12878.hg38.cram    1703
...

Downstream Use¶

The count \(n_i\) is used in Step 4 (Diploid Copy Number Estimation) as the numerator of the diploid CN estimator. For each sample \(i\), its locus read count is scaled relative to the mean count of its nearest neighbors \(\mathcal{N}_i\):

\[\widehat{\text{dipCN}}_i = 2 \cdot \frac{n_i / \bar{d}_i}{\frac{1}{|\mathcal{N}_i|}\sum_{j \in \mathcal{N}_i} n_j / \bar{d}_j}\]

where \(\bar{d}_i\) is the genome-wide depth scale for sample \(i\) derived from the mosdepth normalization step. See Diploid Copy Number Estimation for the full derivation.