Multi-Sample Comparison Guide

Ran the pipeline on two or more people? This guide explains how to compare their results for carrier screening, family planning, and inherited disease investigation.

Privacy between the samples of one data directory. The container of every step script gets the whole GENOME_DIR read-only at /genome (the references, the databases, and every sample's folder), and only its own sample's folder writable. So a container is not a confidentiality boundary between the samples of one data directory: any tool run for one person can read the other people's files. To keep two people's data apart, give each their own GENOME_DIR.

---

Family analysis (step 38, opt-in)

When two or more related people are processed, step 38 analyses them together: one joint call, a check that the samples are the people the pedigree says, and for a child with a sequenced parent, which allele came from which parent. It runs only when you give it a pedigree; a single-sample run is not changed.

The pedigree

A tab-separated file with a header line and one row per person (assets/pedigree.example.tsv):

sample  sex father  mother  partner
child_a female  father_a    mother_a    .
father_a    male    .   .   mother_a
mother_a    female  .   .   father_a

sample is the sample name the pipeline uses. One row per sample you sequence, or will. sex is male, female or unknown. father, mother and partner name another person, or . (or empty) for none. A parent or partner named without a row of their own is a person who was not sequenced: a child's unsequenced father can still be named. The step stops, before any work, when a column name is unknown, a person has two rows, a father is female or a mother male, someone is their own parent or ancestor, two people name different partners, the samplesheet's sex contradicts the pedigree's, a gVCF holds another sample name than its own, or fewer than two named people have a gVCF.

Several families in one file

A families file is the same table with a leading family column, so one file says who is related to whom across all your samples (assets/families.example.tsv):

family  sample  sex father  mother  partner
family_a    child_a female  father_a    mother_a    .
family_a    father_a    male    .   .   mother_a
family_a    mother_a    female  .   .   father_a
family_b    partner_x   male    .   .   partner_y
family_b    partner_y   female  .   .   partner_x

A child is a row whose father or mother is filled; partners name each other. Each family gets its own folder, $GENOME_DIR/<family>/ (or <outdir>/<family>/ with --family_name), so a family name follows the sample-name rule and is nobody's name. A file without the column is one family and reads exactly as before. On top of the pedigree checks, a families file stops before any work, naming the file and line, when the column is filled on some rows only, a person is in two families, a family name breaks the sample-name rule or is a person's name, or a row names a parent or partner who has a row in another family. A person therefore belongs to one family: a child who has a family of their own appears in only one of the two.

Running it

With a families file, run-all.sh runs the family steps for you:

GENOME_DIR=/path/to/data FAMILIES=families.tsv ./scripts/run-all.sh partner_x male
GENOME_DIR=/path/to/data FAMILIES=families.tsv ./scripts/run-all.sh partner_y female

Each call processes its own sample only; the other members are never run as part of it. Before the pipeline starts, the sex given must be the one the file gives the sample (exit 2 otherwise). After the pipeline and before the reports, bin/family_report.py plan looks at what the family has under GENOME_DIR (each member's gVCF with its index, step 06's ClinVar hits, step 26b's panel) and writes $GENOME_DIR/<family>/family_plan.tsv, one row per analysis: joint_call, trio_phasing, inheritance_map, shared_segments (step 38) and couple_carrier per couple (step 39), each runs, skipped or waiting, with the reason in words. Step 38 runs when two or more members have a gVCF (a member still without one is left out and named); step 39 runs when both partners have ClinVar hits, with or without a gVCF (it never reads one). Otherwise the step is recorded in logs/run_status.tsv as skipped (waiting: partner_y has no gVCF yet: run run-all.sh partner_y first) for step 38, skipped (waiting: partner_y has no ClinVar hits yet (step 06 writes clinvar/partner_y_clinvar_hits.vcf): run run-all.sh partner_y first) for step 39, or skipped (<reason>), so the first member's run waits and the last member's run does the family. A step that ran is recorded ok or failed, with its log in logs/38-family.log and logs/39-couple.log. Run step 38 again by hand after a later member is processed if it ran without them.

By hand, from the scripts, once each member has a gVCF from step 03 (<sample>/vcf/<sample>.g.vcf.gz) and somalier's sites are installed (setup.sh --sample-qc-data):

GENOME_DIR=/path/to/data ./scripts/38-family.sh my_family pedigree.tsv

The results go to $GENOME_DIR/my_family/, a folder of its own: the family name follows the sample-name rule and may not be a sample's name. Without the second argument the step reads $GENOME_DIR/my_family/pedigree.tsv. The second argument may also be a families file: the step then copies it to $GENOME_DIR/my_family/families.tsv and reads the rows of family my_family, written out as $GENOME_DIR/my_family/pedigree.tsv. Step 39 takes the same two arguments.

In the Nextflow pipeline, give each member a samplesheet row and add --pedigree pedigree.tsv --somalier_sites <sites.hg38.vcf.gz>. A person processed earlier is a VCF row with the gvcf and gvcf_index columns; a row from FASTQ or a BAM is called in the same run and its gVCF is used. The family results go to <outdir>/family/, or to <outdir>/<name>/ with --family_name <name>. --pedigree may be a families file: with one family it runs that family, with more --family_name names the one to run, and without it the run stops before any task and lists the families. FAMILY_PREPARE writes the same family_plan.tsv as run-all.sh does, for the members of the samplesheet.

The shared segments (item 6 below) need the phased reference panel and genetic maps of step 26b (setup.sh --ancestry-fine-panel). Step 38 reads them from ANCESTRY_FINE_DIR (default reference/ancestry_fine/<ANCESTRY_FINE_PANEL>) and skips the segments with one line when they are not there; the Nextflow family workflow runs them when --ancestry_fine_panel is given. To add the segments to a family analysed before the panel was installed, without joint-calling again:

GENOME_DIR=/path/to/data FAMILY_IBD_ONLY=true ./scripts/38-family.sh my_family

It reads the family.vcf.gz and pedigree.tsv the earlier run left in $GENOME_DIR/my_family/.

What it does

1. Joint call. GLnexus (GLNEXUS_IMAGE, its DeepVariantWGS preset) merges the members' DeepVariant gVCFs into one VCF, family.vcf.gz. A site where one person has a variant gets the others' genotypes from their reference blocks, so "not called" and "matches the reference" are told apart. The container gets 16 GB (GLNEXUS_MEM in step 38, the GLNEXUS process in conf/base.config), and GLnexus's own --mem-gbytes budget is half of it. That budget sizes only GLnexus's database cache and write buffers; the alleles it discovers and the sites it merges are held outside it. So a budget near the limit leaves no room for them. Whether 16 GB is enough for three whole genomes under the half rule has not been measured yet; when the task is killed for memory, Nextflow retries it once with 32 GB (capped at --max_memory), and step 38 takes a larger GLNEXUS_MEM.
2. Relationship check. somalier extracts each member's genotypes at its sites from the gVCF and somalier relate measures the relatedness of every pair: about the share of the genome two people inherited from common ancestors (1 the same person, 0.5 parent and child or full siblings, 0.25 second degree, 0 unrelated). The pedigree gives the expected value. A pair must fall in the same degree as declared (bins at 0.707, 0.354, 0.177 and 0.088); a declared parent and child must also be opposite homozygotes at under 1% of the sites, which siblings are not; a pair declared unrelated or distant must not measure first degree. A pair that disagrees stops the step with exit 1, before any phasing, and relationships.tsv says which pair and why: two samples were swapped, or a pedigree line is wrong. A pair compared at fewer than 1,000 sites is reported as not checked.
3. Phasing by inheritance. For each child whose two parents are sequenced, WhatsHap (WHATSHAP_IMAGE, whatshap phase --ped) phases the family by inheritance alone; no reads are given. A site where the child and both parents are heterozygous stays unphased, because no genotype tells it.
4. The report. For each child with a sequenced parent, <child>_inheritance.tsv.gz gives, at every site of the joint VCF where the child and a sequenced parent have a genotype, the allele from the father and the allele from the mother when the genotypes tell (how is genotype), or ambiguous (both parents heterozygous like the child) or inconsistent (a genotype error, a deletion or a new mutation). On chrX and chrY outside the pseudoautosomal regions a son's X comes from his mother and his Y from his father; the child's sex must be known there. chrM is left out. With one parent sequenced, the allele the other parent passed on is known wherever the sequenced parent is homozygous or differs from the child: that is the half of the unsequenced parent's genome the child carries, and family_report.md gives its share of the sites. The half the child did not inherit is not seen.

De novo candidates (de_novo.tsv, a child allele neither parent has) are called only when both parents are sequenced. A candidate is pass with GQ 20 or more and depth 10 or more in all three and, when the child is heterozygous, 25% to 75% of the child's reads on the new allele; otherwise low. Most candidates of a whole genome are genotype errors: confirm one that matters by another method.

When step 23 (the clinical filter) ran for a child, every pair of its candidate variants the child carries heterozygous in one gene is listed in compound_het.tsv: in trans when one came from each parent (both copies of the gene are hit), in cis when both came from one parent, unknown when the genotypes do not tell.

1. The chromosome map. For each child with a sequenced parent, bin/family_segments.py (run by the report) maps the stretches of each chromosome the child got from that parent. It reads the joint VCF's genotypes at the sites where the parent is heterozygous and the genotypes give the allele passed on. When one of that parent's own parents is sequenced too, each such site also says which of the parent's two copies was passed on, the grandfather's or the grandmother's, and a switch from one to the other along a chromosome is a crossover. Without a sequenced grandparent the stretches say only "from this parent". With one parent sequenced, the child's other copy is written over the same stretches as the copy that did not come from the sequenced parent.

Isolated sites that point to the other copy are genotype errors or short gene conversions: a run of them becomes a segment of its own only with 20 sites or more over 500 kb or more; a shorter run is absorbed and counted in the segment's discordant column. A stretch with no informative site for 5 Mb (a centromere, a gap of the reference) is split there, and a stretch with fewer than 20 informative sites is left out. A segment starts and ends at its first and last informative site, so a crossover lies somewhere in the space between two segments. Autosomes are mapped, and chrX outside the pseudoautosomal regions from the mother only (a father passes his single X whole to a daughter, and none to a son). A section of family_report.md counts each child's segments, Mb and crossovers per parent.

<child>_dnapainter.csv holds the same segments in the import format of DNA Painter (Chr,Start,End,cM,SNPs,Match,Confidence,Group,Side,Notes,Colour,Ancestor; Side is paternal or maternal; Group is the grandparent the segment came from, or the parent when no grandparent is sequenced, so each grandparent is painted as a group of its own; Ancestor is 1 when the group names a parent or grandparent, sequenced or not, as DNA Painter's format asks). Positions are GRCh38. DNA Painter's import page names no genome build, and matches imported from testing companies usually carry build 37 positions, so a segment here and a match from such a company are offset by the difference between the builds. cM is left empty: no genetic map is used. Importing needs a DNA Painter subscription.

1. Shared segments. On each autosome the panel covers, the joint call's genotypes at the panel's SNVs (same position and alleles; the panel keeps at most one common SNV per 2 kb) are phased by Beagle 5.5 (BEAGLE_IMAGE; Browning, Zhou and Browning 2018, https://doi.org/10.1016/j.ajhg.2018.07.015) against the panel's phased genotypes. hap-ibd (HAP_IBD_IMAGE; Zhou, Browning and Browning 2020, https://doi.org/10.1016/j.ajhg.2020.02.010) then finds the segments each pair of members shares identical by descent, haplotype by haplotype, with its defaults: a segment of 2 cM or more on Beagle's GRCh38 genetic map. bin/ibd_report.py merges a pair's segments into the stretches the two share, one copy or both, as consumer DNA sites show them. A phasing error breaks hap-ibd's segments, and a piece shorter than its 2 cM seed between two breaks is lost. Two people who share a copy are never opposite homozygotes (one 0/0, the other 1/1) bar a genotype error, so two stretches of a pair are joined when the markers between them hold at most one opposite homozygote, and a pair's first and last stretch reach the chromosome's first and last marker on the same test. The Browning lab, which wrote hap-ibd, joins segments by the same test with its merge-ibd-segments tool, in gaps under 0.6 cM (https://faculty.washington.edu/browning/refined-ibd.html). Here the gap has no length limit: two unrelated people are opposite homozygotes at several percent of the panel's common SNVs, and the panel holds about 200 of them per cM, so a gap of 1 cM or more with at most one is almost always a stretch they share. A shorter gap passes by chance more easily, but it is only ever a gap between stretches hap-ibd already found, or between one of them and the end of the chromosome. Per pair it gives the number of stretches, the total and longest cM, the cM where both copies are shared (IBD2), and the share of the analysed map (the cM between the first and last phased marker of each chromosome). The relatedness is (IBD1 / 4 + IBD2 / 2) / analysed cM, times 2, binned as KING bins it (0.707, 0.354, 0.177, 0.088, plus 0.044 for a fourth degree). A first-degree pair sharing both copies over 10% or more reads as full siblings, otherwise as parent and child. Unlike the relationship check this estimate does not stop the step: compare its degree with the declared column.

Two members marked as partners also get the couple carrier screen of step 39: scripts/39-couple.sh, or in the Nextflow pipeline the COUPLE_CARRIER process of the family workflow (with clinvar in --tools).

Outputs

The Nextflow pipeline writes them to <outdir>/family/ (<outdir>/<name>/ with --family_name), step 38 to $GENOME_DIR/<family_name>/ (with the pedigree.tsv it read, and families.tsv when it was given a families file).

File	Content
family.vcf.gz (+ .tbi)	the joint call
members.tsv, family.ped, trios.ped	the pedigree as read, somalier's pedigree, WhatsHap's trios
family_plan.tsv	family, analysis, members, decision (runs, skipped, waiting) and reason, for each analysis of steps 38 and 39
pedigree.tsv	the family's own rows (from a families file, without the family column)
somalier/family.pairs.tsv, family.samples.tsv, family.html	somalier relate's own tables and plot
relationships.tsv	each pair: declared relationship, expected and measured relatedness, opposite homozygotes, sites, verdict and reason
family.phased.vcf.gz (+ .tbi)	WhatsHap's phasing (only when a child has both parents sequenced)
<child>_inheritance.tsv.gz	allele from each parent, site by site
<child>_segments.tsv	the chromosome map: child, side, parent, chrom, start, end, length_bp, sites, discordant, ancestor (the grandparent, or .), basis (parent and grandparents, parent, or other copy)
<child>_dnapainter.csv	the same segments for DNA Painter's import (GRCh38 positions)
<child>_segments.svg	the map drawn: per chromosome, the copy from the father above the copy from the mother
de_novo.tsv, compound_het.tsv	as above
family_report.md	all of it in words, pair by pair and child by child
ibd_pairs.tsv	each pair: declared relationship, expected relatedness, stretches, shared_cm, ibd2_cm, longest_cm, analysed_cm, shared_fraction, relatedness, degree (with the panel only)
ibd_segments.tsv	each shared stretch: the pair, chrom, start, end (GRCh38), cm, ibd2_cm
ibd_report.md	the shared segments in words
ibd/	step 38 only: the genotypes at the panel's SNVs, Beagle's phased VCF and hap-ibd's own output, per chromosome

On the live site

The live site (bin/owngenome_live.py) lists each family folder under GENOME_DIR as a run of its own, of kind family: any folder holding members.tsv or pedigree.tsv. The plugin bin/owngenome_live_family.py and its page script assets/live/plugins/family.js add a Family card to that run's page:

• The pedigree. bin/pedigree_svg.py draws it from members.tsv, or from pedigree.tsv when step 38 has not run, each time the page asks. A square is a male, a circle a female, a diamond a sex not known; a filled symbol is a sequenced person (with only pedigree.tsv, a person with a row of their own). A line joins partners, and the children hang from it, one generation per row. It draws what the file declares; the relationship check below is what tests it.
• Each child's chromosome map, <child>_segments.svg, under the child's name, or a line saying step 38 has not drawn one.
• The relationship check (relationships.tsv), the shared segments (ibd_pairs.tsv) and the plan (family_plan.tsv, when there is one), row by row as the steps wrote them.
• The couple carrier screen, per couple: how many genes couple_<a>_<b>.tsv lists as shared and as x_linked, and a link to couple_<a>_<b>.md. A count of 0 is not an all-clear; the report says what the screen does not see.

A sample's own page gets a link to the family it is in. The plugin reads the files the steps wrote and nothing else: it calls no tool and computes no genetics. It serves only the pedigree it draws and the result files named above, at /family/<family>/<file>, behind the site's login; any other file of the folder, such as family.vcf.gz, is 404. tests/test_pedigree_svg.py and tests/test_owngenome_live_family.py check the drawing and the card on an invented family (tests/fixtures/families_live/), with a swapped parent as the control that must fail.

python3 bin/pedigree_svg.py "$GENOME_DIR/my_family" > pedigree.svg   # the same drawing, without the site

Limits

• The gVCFs must be DeepVariant's: the GLnexus preset reads DeepVariant's fields.
• run-all.sh decides from the files under GENOME_DIR, not from whether a member's run succeeded: a member whose gVCF is there counts as ready. Step 38 runs on the members ready at that moment; a member processed later needs step 38 run again (the next member's run-all.sh does it, or 38-family.sh by hand).
• In the Nextflow pipeline the folder is family/ unless --family_name names it, also when a families file holds one family.
• No reads are used for phasing, so an all-heterozygous trio site stays unphased, and a duo is phased only where the genotypes tell.
• somalier's relatedness is one number per pair for the whole genome. The chromosome map covers a child and its parents (and grandparents). Stretches shared by any other pair, such as siblings or cousins, come from the shared segments, which need step 26b's panel.
• The shared segments see only stretches of 2 cM or more, and only on the chromosomes the panel covers (chr1 to chr22). Phasing errors can split a long stretch into shorter ones, so the count of stretches can be too high and the longest stretch too short. Pairs from one endogamous population share many short stretches by chance. The degree estimate of a pair beyond the third degree is rough, and it is rough for any pair when fewer than the 22 autosomes are analysed.
• The chromosome map needs informative sites: a crossover is placed only between two of them, and with short reads and no sequenced grandparent nothing tells the parent's two copies apart. No genetic map is used, so lengths are in bases, not centimorgans. Two siblings with no sequenced grandparent are not compared with each other.
• Runtime on whole genomes has not been measured yet; the CI test runs the joint call on 2.45 Mb of the GIAB Ashkenazi trio, and the shared segments on all of chr20 of the same trio (GIAB's benchmark genotypes as the joint call, FAMILY_IBD_ONLY=true).

---

Carrier Screening for Partners

The most immediately useful multi-sample analysis: checking whether both partners carry pathogenic variants in the same recessive gene. If so, each child has a 25% chance of being affected.

Step 39 does this for two partners of a pedigree: the recessive genes where both carry a change, counting ClinVar hits, Delly deletions of 50 bp or more over an exon (a long one only when read depth shows the loss, and none that matches a common deletion of AnnotSV's benign set), SMN1 copy loss and FMR1 and FXN expansions, the X-linked findings, and what the screen does not see. The cross-check below is the same comparison by hand, for ClinVar hits only.

Quick Cross-Check

source versions.env   # from the repository root
PARTNER_A="sample_a"
PARTNER_B="sample_b"

# Gene symbols of each partner's ClinVar hits (step 6 writes them with ClinVar's GENEINFO)
for SAMPLE in $PARTNER_A $PARTNER_B; do
  docker run --rm -v "${GENOME_DIR}:/genome" "${BCFTOOLS_IMAGE}" \
    bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/GENEINFO\t%INFO/CLNSIG\t[%GT]\n' \
      /genome/${SAMPLE}/clinvar/${SAMPLE}_clinvar_hits.vcf \
    > /tmp/${SAMPLE}_clinvar_genes.txt
done

# GENEINFO is SYMBOL:GeneID, several joined by |
for SAMPLE in $PARTNER_A $PARTNER_B; do
  cut -f5 /tmp/${SAMPLE}_clinvar_genes.txt | tr '|' '\n' | cut -d: -f1 | grep -v '^\.$' | sort -u \
    > /tmp/${SAMPLE}_genes.txt
done

# Confirm each per-sample gene file is non-empty before reading the intersection
wc -l /tmp/${PARTNER_A}_genes.txt /tmp/${PARTNER_B}_genes.txt

# Find genes where BOTH partners have pathogenic hits
comm -12 /tmp/${PARTNER_A}_genes.txt /tmp/${PARTNER_B}_genes.txt

Before reading the result, check the wc -l line. A gene file with 0 lines means that partner has no hit with a gene name, or the hits file is missing or from an older run. An empty intersection then says nothing about carrier risk. Each file should list the genes of that partner's hits; rerun step 6 if it does not.

If both gene files have lines and the intersection is empty: No shared recessive carrier risk was found among the small variants in ClinVar. This is the most common result. It does not cover copy-number carriers such as SMA (see below).

If genes appear in both lists: Check whether: 1. Both variants are in the same gene (not just nearby genes) 2. Both are classified as Pathogenic or Likely pathogenic (not VUS) 3. The gene follows autosomal recessive inheritance (check OMIM or ClinVar) 4. Both partners are heterozygous (carriers), not homozygous

Common Carrier Genes (Not Usually a Concern)

These genes frequently show up in ClinVar carrier screens. Being a carrier is very common and only relevant if your partner carries the same gene:

Gene	Condition	Carrier Frequency
CFTR	Cystic fibrosis	1 in 25 (European)
GJB2	Hearing loss (DFNB1)	1 in 30
HFE	Hemochromatosis	1 in 10 (H63D), 1 in 150 (C282Y)
MUTYH	Colorectal cancer risk	1 in 50
HEXA	Tay-Sachs disease	1 in 30 (Ashkenazi), 1 in 300 (general)

Spinal muscular atrophy is not in this table on purpose. Most SMA carriers have one copy of SMN1 instead of two: a copy-number loss, not a small variant. ClinVar screening of a VCF cannot see it, so this check says nothing about SMA carrier status.

The opt-in step 35 (Parascopy) estimates SMN1 and SMN2 copy number from each partner's BAM. What it can and cannot say for a couple:

• Assessed: each partner's SMN1 copy number, with a quality. One partner with one SMN1 copy (quality 20 or more, filter PASS) is a likely carrier; both partners with one copy is the pattern that gives a 1 in 4 chance of an affected child.
• Not assessed: a "2+0" carrier, two SMN1 copies on one chromosome and none on the other, reads as two copies like a non-carrier; small variants inside SMN1; the other paralog genes (GBA, CYP21A2, HBA1/HBA2 and so on).

A result that matters for family planning needs a clinical SMN1 carrier test; ask for one that reports the 2+0 risk.

---

Comparing Pharmacogenomics

PharmCAT results can differ dramatically between partners. Compare the HTML reports side by side for genes that affect commonly prescribed medications. HLA-A, HLA-B and CYP2D6 come from the BAM-based callers through step 36; CYP2D6 has a result only when pypgx and Cyrius (opt-in) agree, so a partner whose CYP2D6 reads indeterminate has no comparable result, not a normal one:

Gene	One Partner is Rapid, Other is Poor?	Clinical Impact
CYP2C19	PPIs, SSRIs, clopidogrel	Different SSRI dosing needed
CYP2D6	Codeine, tramadol, psych meds	Codeine dangerous for poor metabolizers
NAT2	Isoniazid, hydralazine	Slow acetylators get more side effects at standard doses
CYP2C9	Warfarin, NSAIDs	Warfarin dosing differs

---

Parent-Child Analysis

If you have WGS data for a parent and child, you can investigate:

Inherited vs De Novo Variants

A de novo variant is one that appeared for the first time in the child (not present in either parent). These are rare (~50-100 per genome) and occasionally clinically significant.

source versions.env   # from the repository root
PARENT="parent_name"
CHILD="child_name"

# Find variants in the child that are NOT in the parent
docker run --rm -v "${GENOME_DIR}:/genome" "${BCFTOOLS_IMAGE}" \
  bcftools isec -C \
    /genome/${CHILD}/vcf/${CHILD}.vcf.gz \
    /genome/${PARENT}/vcf/${PARENT}.vcf.gz \
    -p /genome/${CHILD}/vcf/de_novo_candidates/

Note: This is a rough screen. With both parents sequenced, step 38 calls de novo candidates from a joint call instead. True de novo detection requires trio analysis (both parents + child) and careful filtering for sequencing errors. Many variants flagged by bcftools isec will be false positives (present in the parent but missed by the variant caller due to low coverage at that position).

Carrier Inheritance Tracing

If the child is a carrier for a recessive condition, you can check which parent contributed the variant:

source versions.env   # from the repository root
GENE_REGION="chr13:20189473-20189473"  # Example: GJB2 position

for SAMPLE in $PARENT $CHILD; do
  echo "--- ${SAMPLE} ---"
  docker run --rm -v "${GENOME_DIR}:/genome" "${BCFTOOLS_IMAGE}" \
    bcftools view -r "$GENE_REGION" /genome/${SAMPLE}/vcf/${SAMPLE}.vcf.gz
done

---

Structural Variant Comparison

SVs called by multiple callers in one person have lower false-positive rates. SVs shared between family members add further confidence:

source versions.env   # from the repository root
# Compare Manta SVs between two samples
# (Simple overlap check using bedtools-style comparison)
for SAMPLE in $PARTNER_A $PARTNER_B; do
  docker run --rm -v "${GENOME_DIR}:/genome" "${BCFTOOLS_IMAGE}" \
    bcftools query -f '%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\n' \
      /genome/${SAMPLE}/manta/results/variants/diploidSV.vcf.gz \
    > /tmp/${SAMPLE}_svs.bed
done

echo "Partner A SVs: $(wc -l < /tmp/${PARTNER_A}_svs.bed)"
echo "Partner B SVs: $(wc -l < /tmp/${PARTNER_B}_svs.bed)"

# Exact position matches (most stringent)
comm -12 <(sort /tmp/${PARTNER_A}_svs.bed) <(sort /tmp/${PARTNER_B}_svs.bed) | wc -l
echo "Shared SVs (exact match)"

Expected: Partners (unrelated) share very few SVs at exact positions. Parent-child pairs share ~50% of SVs.

---

Mitochondrial Haplogroup Comparison

Relationship	Expected mtDNA Result
Partners	Different haplogroups (unless same maternal lineage)
Siblings	Identical haplogroup (same mother)
Mother-child	Identical haplogroup
Father-child	Different haplogroups (mtDNA is maternal only)

If siblings have different mitochondrial haplogroups, it may indicate different biological mothers (adoption, etc.) or a very rare paternal mtDNA inheritance event.

The Y-chromosome haplogroup (step 37, male samples) is the paternal counterpart: a father and his sons, and brothers, share it.

---

Telomere Length Comparison

Telomere length (step 10) is most informative when compared between samples of similar age sequenced on the same platform:

• Partners of similar age: Should be roughly similar. Significant differences (>20%) may reflect lifestyle, stress, or genetic factors.
• Parent-child: Parent will generally have shorter telomeres. The difference correlates loosely with age gap.
• Siblings: Similar telomere lengths expected. Large differences may be worth investigating with a clinical telomere assay.

Caution: TelomereHunter's tel_content is a rough estimate. Only use it for relative comparisons between samples processed identically. Do not compare with values from other labs, papers, or sequencing platforms.

---

What This Pipeline Cannot Do (Yet)

• Read-backed or statistical phasing (which variants are on which chromosome copy in one person): step 38 phases a child only by inheritance, from the parents' genotypes; a single person needs read-based phasing (long reads help) or statistical phasing tools like Eagle2/SHAPEIT
• Comparable polygenic risk scores: step 25 computes raw scores per person, but turning them into percentiles that can be compared between people needs a matched reference population (see step 25)
• Ancestry estimates: one sample at a time. Step 26 projects each sample onto a 1000 Genomes panel and step 26b estimates its admixture fractions and paints its chromosomes; neither compares two samples with each other, so read their outputs side by side

These are planned for future pipeline versions.