Getting Started

GRiD is driven by a single YAML configuration file. This page walks through every section — what each parameter does, why it is needed, and how to obtain or prepare the required input files.

Copy the bundled example and edit it for your cohort:

cp grid/example_config.yaml my_config.yaml
# edit my_config.yaml, then:
grid wgs my_config.yaml

---

Global Settings

These fields apply to every pipeline step.

samples_file: "path/to/sample_ids.txt"
directory_loc: "path/to/crams/"
reference_genome: "path/to/hg38.fa"
output_dir: "path/to/output/"
threads: 8
file_type: "cram"
chrom: "chr6"
start_bp: 160605062
end_bp: 160647661
output_file_type: "tsv"

Parameter	Type	Description
samples_file	path	Plain-text file with one sample ID per line. IDs must match the file stems of the CRAM/BAM files in directory_loc.
directory_loc	path	Directory containing the CRAM (or BAM) files for your cohort.
reference_genome	path	Reference genome FASTA (hg19 or hg38). Must be indexed (samtools faidx). Download from UCSC, the Broad Institute, or 1000 Genomes.
output_dir	path	Directory where all pipeline outputs are written. Created if it does not exist.
threads	int	Number of CPU threads shared across steps. Match to your SLURM --cpus-per-task allocation.
file_type	str	Sequencing file format: "cram" or "bam".
chrom	str	Target chromosome in the same notation as your CRAM headers (e.g. "chr6" for hg38, "6" for hg19).
start_bp	int	Start coordinate (0-based) of the VNTR region of interest. For LPA KIV-2 in hg38: 160605062.
end_bp	int	End coordinate of the VNTR region of interest. For LPA KIV-2 in hg38: 160647661.
output_file_type	str	Output table format: "tsv" or "csv".

Note

Supported VNTR loci. Mukamel et al. (2021) Science (doi:10.1126/science.abg8289) catalogued 734 coding VNTRs with elevated IBD2 rates across the human genome, identifying loci where copy number is likely to vary and be heritable. Their full code is available at Zenodo 4776804.

The same file is bundled with GRiD at 734_possible_coding_vntr_regions.IBD2R_gt_0.25.uniq.txt. To run GRiD on a different locus, simply update the chrom, start_bp, and end_bp fields in your config file to match the coordinates of your target region. For example, to run on LPA KIV-2.

---

Step 1 — index

Verifies or creates .crai / .bai index files for all CRAM/BAM files in directory_loc. Only needed if your files are not already indexed.

index:
  run: False
  output_file_prefix: "index_file_results"

Parameter	Description
run	Set True to run this step. Skip if files are already indexed.
output_file_prefix	Prefix for the log output file.

Step 2 — count_reads

Counts properly paired reads overlapping the target VNTR region from each CRAM/BAM file. The resulting per-sample read counts are used downstream as a proxy for locus-specific sequencing depth.

count_reads:
  run: True
  output_file_prefix: "read_counts"
  min_mapq: 1
  flags:
    - 83    # proper pair, read reverse strand
    - 147   # proper pair, mate reverse strand
    - 81    # read reverse strand
    - 145   # mate reverse strand

Parameter	Description
run	Set True to run this step.
output_file_prefix	Prefix for the output read-count file.
min_mapq	Minimum mapping quality. Reads below this threshold are ignored. 1 excludes unmapped reads while retaining multi-mappers; higher values (e.g. 20) restrict to uniquely mapped reads.
flags	List of SAM bitwise flags to include. Only reads whose flag exactly matches one of the listed values are counted. The defaults capture properly paired reads on both strands. See SAM flag documentation for flag definitions.

Step 3 — mosdepth

Runs mosdepth in windowed mode to compute binned read depth across the full genome for each sample. The resulting depth profiles are used for normalization and neighbor finding.

mosdepth must be available on your PATH:

conda install -c bioconda mosdepth
# or: module load mosdepth  (HPC)

mosdepth:
  run: True
  output_file_prefix: "mosdepth_results"
  bin_size: 1000
  mode: "fast"
  region_name: "LPA"
  work_dir: "path/to/mosdepth_work/"
  remove_intermediate: True

Parameter	Description
run	Set True to run mosdepth.
output_file_prefix	Prefix for the aggregated depth output file.
bin_size	Genomic bin size in base pairs. 1000 (1 kb) is recommended — smaller bins increase resolution but memory and runtime.
mode	"fast" skips per-base depth (faster, less memory) or "full" for per-base output. "fast" is sufficient for normalization.
region_name	Label applied to the VNTR region in output files (e.g. "LPA").
work_dir	Scratch directory for per-sample mosdepth intermediate files. Should be on fast local or scratch storage.
remove_intermediate	If True, per-sample mosdepth files are deleted after aggregation to save disk space.

Step 3a — mosdepth.normalize

Normalizes the depth matrix: filters low/high-coverage bins, applies a repeat-mask, and z-scores across individuals. See Coverage Normalization for the full method.

mosdepth:
  normalize:
    run: True
    min_depth: 20
    max_depth: 100
    top_frac: 0.1
    output_file_prefix: "mosdepth_results_normalized"
    repeat_mask_file: "path/to/repeat_mask.bed"

Parameter	Description
run	Set True to normalize after mosdepth.
min_depth	Bins with median depth below this value are excluded.
max_depth	Bins with median depth above this value are excluded.
top_frac	Fraction of bins with the highest variance to retain for normalization (e.g. 0.1 keeps the top 10 % most variable bins).
output_file_prefix	Prefix for the normalized depth matrix output.
repeat_mask_file	BED file of repetitive regions to exclude before normalization. Download UCSC RepeatMasker tracks from the UCSC Table Browser (group: Repeats, track: RepeatMasker).

Step 3b — mosdepth.neighbors

Finds read depth-matched nearest neighbors in the normalized depth space using scikit-learn NearestNeighbors. See Nearest-Neighbor Identification.

mosdepth:
  neighbors:
    run: True
    output_file_prefix: "neighbor_coverage"
    num_neighbors: 5
    zmax: 2.0
    sigma2_max: 1000

Parameter	Description
run	Set True to run neighbor finding.
output_file_prefix	Prefix for the neighbor assignment output file.
num_neighbors	Number of nearest neighbors to retain per sample. 5 is a conservative default; increase for small or heterogeneous cohorts.
zmax	Clip normalized depth values to \(\pm\) zmax before distance computation. Reduces the influence of outlier bins. Typical range: 1.5–3.0.
sigma2_max	Maximum allowed neighbor variance. Samples whose nearest-neighbor set has variance above this threshold are flagged. Increase if many samples are being excluded in diverse cohorts.

Step 4 — compute_diploid_genotypes

Estimates diploid copy number for each sample using the neighbor-normalized depth ratio. For each sample, its locus depth is scaled by the mean depth of its nearest neighbors at the same locus. See Diploid Copy Number Estimation.

compute_diploid_genotypes:
  run: True
  output_file_prefix: "diploid_genotypes"
  n_nbr: 300

Parameter	Description
run	Set True to estimate diploid copy numbers.
output_file_prefix	Prefix for the diploid CN output file.
n_nbr	Number of nearest neighbors to use when computing each sample’s diploid CN estimate. Default: 300. Increase for small cohorts; decrease to restrict to the closest matches.

Step 5 — compute_haploid_genotypes

Decomposes diploid CN estimates into haplotype-specific copy numbers using an iterative phasing algorithm. Two neighbor sources are supported via method. See Haplotype Inference and Computing IBS/IBD Neighbors.

compute_haploid_genotypes:
  run: True
  output_file_prefix: "haploid_genotypes"
  method: "ibs"          # "ibs" or "ibd"
  min_neighbors: 1
  max_neighbors: 10
  n_iters: 100

  # IBS method
  ibs_output: "path/to/ibs_neighbors_chr6.tsv.gz"

  # IBD method
  # ibd_output: "path/to/ilash_output.txt"
  # weighted: False
  # weight_scale: 1000000
  # min_length: 0.5
  # min_match: 0.70

Parameter	Description
run	Set True to run haplotype inference.
output_file_prefix	Prefix for the haploid CN output files.
method	Neighbor source: "ibs" (computeIBSpbwt output) or "ibd" (iLASH output).
min_neighbors	Minimum neighbors required per haplotype. Samples below this threshold are excluded from iterative phasing.
max_neighbors	Maximum neighbors used per haplotype per iteration.
n_iters	Iterations for the update loop. 100 is a safe default.
ibs_output	(IBS method) Path to IBS neighbors file from computeIBSpbwt.
ibd_output	(IBD method) Path to iLASH output file.
weighted	(IBD method) If True, weight each IBD segment by distance to the target region and match score (Lorentzian kernel). Default: False.
weight_scale	(IBD method, weighted) Distance half-width in bp for the Lorentzian weight function. Default: 1000000.
min_length	(IBD method) Minimum IBD segment length in cM. Default: 0.5.
min_match	(IBD method) Minimum iLASH match score. Default: 0.70.

Note

Either ibs_output or ibd_output must be provided depending on the chosen method. See Computing IBS/IBD Neighbors for file format details and IBS/IBD Neighbor Computation Example for an annotated example script.

---

Full Example Config

grid/example_config.yaml

samples_file: "path/to/sample_ids.txt"
directory_loc: "path/to/crams/"
reference_genome: "path/to/hg38.fa"
output_dir: "path/to/output/"
threads: 4
file_type: "cram"

# VNTR region of interest
chrom: "chr6"
start_bp: 160605062
end_bp: 160647661

output_file_type: "tsv"

index:
  run: False
  output_file_prefix: "index_file_results"

count_reads:
  run: True
  output_file_prefix: "read_counts"
  min_mapq: 1
  flags:
    - 83   # proper pair, reverse strand
    - 147  # proper pair, mate reverse
    - 81   # reverse, mate forward (unpaired-ish)
    - 145  # mate reverse

mosdepth:
  run: True
  output_file_prefix: "mosdepth_results"
  bin_size: 1000
  mode: "fast"
  region_name: "LPA"
  work_dir: "path/to/mosdepth_work/"
  remove_intermediate: True

  normalize:
    run: True
    min_depth: 20
    max_depth: 100
    top_frac: 0.1
    output_file_prefix: "mosdepth_results_normalized"
    repeat_mask_file: "path/to/repeat_mask.bed"

  neighbors:
    run: True
    output_file_prefix: "neighbor_coverage"
    num_neighbors: 5
    zmax: 2.0
    sigma2_max: 1000

compute_diploid_genotypes:
  run: True
  output_file_prefix: "diploid_genotypes"
  n_nbr: 300

compute_haploid_genotypes:
  run: True
  output_file_prefix: "haploid_genotypes"
  method: "ibs"          # "ibs" or "ibd"
  min_neighbors: 1
  max_neighbors: 10
  n_iters: 100

  # IBS method: provide computeIBSpbwt output
  ibs_output: "path/to/ibs_neighbors_chr6.tsv.gz"

  # IBD method: provide iLASH output and filtering thresholds
  # ibd_output: "path/to/ilash_output.txt"
  # weighted: False         # True = Lorentzian distance+match weighting
  # weight_scale: 1000000   # bp half-width for distance weight
  # min_length: 0.5         # minimum IBD segment length (cM)
  # min_match: 0.70         # minimum IBD match score