HLA calls¶

c4 reads HLA typing from HLA-HD, HIBAG and the lab's standardized files. Every loader gives you the same thing: one column per gene, holding HLAgenotype objects. Each object knows its gene, its two alleles and its resolution. It can print itself in any format, count copies of an allele, and compare itself to another call.

Load HLA calls¶

HLA tables are in the catalog under hla, one block per caller. Load them like any other table:

from c4.cohort import cohort

demo = cohort("demo")
hla = demo.load("hla.hlahd", verbose=False)
hla[["patient_id", "seqID", "A", "B", "DRB1"]]

  patient_id    seqID                A                B                   DRB1
0       P001  EX-0001  A*02:01/A*24:02  B*07:02/B*44:02  DRB1*04:01/DRB1*15:01
1       P002  EX-0002  A*01:01/A*03:01  B*08:01/B*35:01  DRB1*01:01/DRB1*03:01
2       P004  EX-0003  A*02:01/A*02:01  B*44:03/B*51:01  DRB1*07:01/DRB1*13:01
3       P003  EX-0004  A*11:01/A*24:02  B*40:01/B*46:01  DRB1*09:01/DRB1*12:02
4       P005  EX-0005  A*03:01/A*68:01  B*07:02/B*15:01  DRB1*11:01/DRB1*15:01

The gene columns look like text, but each cell is an HLAgenotype. The two alleles are always stored in sorted order, so DRB1*15:01/04:01 in the file becomes DRB1*04:01/DRB1*15:01.

How a file is read depends on its catalog entry's loader:

Loader	Used for	Reads
`HLA`	every entry under `hla:` unless it says otherwise	gene columns `A` (`A*02:01/24:02`) or allele pairs `A_1`, `A_2`
`HLA-HD`	HLA-HD output	same as `HLA`
`HIBAG`	raw HIBAG output	`allele1_A` and `allele2_A`, plus `prob_A`

Only genes in config.hla_genes (A, B, C, DRB1, DQA1, DQB1, DPB1) are turned into genotypes. Any other column is left as it is in the file.

HIBAG probabilities¶

HIBAG output has one file per reference panel. Each genotype loaded from it carries HIBAG's posterior probability as .prob. The prob_A column is kept too:

hibag = demo.load("hla.hibag", panel="african", verbose=False)
hibag.loc[0, "A"].prob

0.97

Work with a genotype¶

gt = hla.loc[0, "A"]
gt

HLAgenotype('A*02:01/A*24:02')

str(gt) gives A*02:01/A*24:02. For another layout, use genotype():

Call	Result
`gt.genotype(format="standard2")`	`'A02:01/A24:02'`
`gt.genotype(format="standard")`	`'A*02:01/24:02'`
`gt.genotype(format="standard_no_gene")`	`'02:01/24:02'`
`gt.genotype(format="tuple")`	`(('02', '01'), ('24', '02'))`
`gt.genotype()`	the two alleles, as `HLAallele` objects

Pass digits to get a lower resolution. Two digits make one field, so 4 is 2-field:

g = HLAgenotype("A*11:01:01/24:02:01")
g.genotype(format="standard2", digits=4)

'A*11:01/A*24:02'

Other useful attributes:

Attribute	Meaning
`gt.a1`, `gt.a2`	the two alleles (`HLAallele`), missing one last
`gt.gene`	`"A"`
`gt.fields`	resolution of the lower-resolution allele, in fields
`gt.missing_alleles`	`0`, `1` or `2`

A missing allele is written ., so A*02:01/. has one missing allele. A blank cell becomes ./..

Allele dosages¶

copies_of() counts how many of the two alleles match, which gives an allele dosage for association tests:

hla["A"].map(lambda gt: gt.copies_of("A*02:01"))

0    1
1    0
2    2
3    0
4    0
Name: A, dtype: int64

Save it as a column with hla["A*02:01"] = ....

The match is made at the resolution of the allele you ask about, so copies_of("A*02:01") also counts A*02:01:01, and copies_of("A*02") counts every A*02 allele (the 1-field allele group).

Rows with no HLA call

After a merge, patients with no row in the HLA table have NaN in the gene columns, not a genotype. Check before calling a method:

from c4.hla import HLAgenotype

df = demo.load("pheno", "hla.hlahd", verbose=False)
df["A*02:01"] = df["A"].map(
    lambda gt: gt.copies_of("A*02:01") if isinstance(gt, HLAgenotype) else float("nan")
)

Compare two callers¶

concordance() scores how well two calls for the same gene agree:

Score	Meaning
`2`	both alleles match
`1`	one allele matches
`0`	neither matches
`NaN`	either call is missing an allele

The two calls are compared at the lower of their resolutions, so A*11:01:01 and A*11:01:02 match as 11:01. Pass digits=4 to compare at no more than 2-field.

To compare HLA-HD and HIBAG patient by patient, load both. They share the column A, so the second table's copy is prefixed with its name:

both = demo.load("hla.hlahd", "hla.hibag", how="inner", verbose=False)
both["A_concordance"] = [
    a.concordance(b) for a, b in zip(both["A"], both["hla.hibag.raw.A"])
]
both[["patient_id", "A", "hla.hibag.raw.A", "A_concordance"]]

  patient_id                A  hla.hibag.raw.A  A_concordance
0       P001  A*02:01/A*24:02  A*02:01/A*24:02              2
1       P002  A*01:01/A*03:01  A*01:01/A*03:01              2
2       P004  A*02:01/A*02:01  A*02:01/A*02:01              2
3       P003  A*11:01/A*24:02  A*11:01/A*24:02              2

For two files that aren't in the catalog, use c4 hla_concordance instead.

Change table layout¶

Some tools want one column per allele (A_1, A_2) rather than one per gene. split_genotype and merge_genotype convert between the two, and both return objects:

from c4.hla import split_genotype, merge_genotype

split = split_genotype(hla, gene_list=["A", "B"])    # A -> A_1, A_2 (HLAallele)
merged = merge_genotype(split, gene_list=["A", "B"]) # A_1, A_2 -> A (HLAgenotype)

Save to a file¶

Writing a table with to_csv saves each genotype as A*02:01/A*24:02. For another format, convert the columns first:

out = hla.copy()
for gene in ["A", "B", "DRB1"]:
    out[gene] = out[gene].map(lambda gt: gt.genotype(format="standard_no_gene", digits=4))
out.to_csv("hla_2field.tsv", sep="\t", index=False)

c4 reads a saved file back into genotypes, in any of these text formats.

Command line¶

The c4 command converts and compares HLA files without writing Python. Every command reads .csv, .tsv or space-separated files. Each file can have genotype columns (A) or allele columns (A_1, A_2).

`c4 split-hla`¶

Genotype columns to allele columns. With no output file, the result is printed:

c4 split-hla calls.csv --genes A,B --columns ID --digits 4

  ID     A_1     A_2     B_1     B_2
S001 A*01:01       . B*07:02 B*08:01
S002 A*02:01 A*03:01       .       .
...

Give an output file to save it instead: c4 split-hla calls.csv calls_split.tsv ....

`c4 merge-hla`¶

Allele columns to genotype columns:

c4 merge-hla calls_split.csv calls.tsv --genes A,B --columns ID

Options for both commands:

Option	Default	Description
`--genes`, `-g`	all canonical genes	Comma-separated genes, e.g. `A,B,DRB1`
`--columns`, `-c`	every column	Other columns to keep, comma-separated
`--digits`, `-d`	all fields	Resolution, e.g. `4` for 2-field
`--format`, `-f`	`standard`	`standard` (`A02:01/24:02`), `standard2` (`A02:01/A*24:02`) or `standard_no_gene` (`02:01/24:02`)

`c4 hla_concordance`¶

Compare two typing files, sample by sample:

c4 hla_concordance wes.csv imputed.tsv \
    --file1_id seqID --file2_id IID \
    --genes A,B --digits 4 --output concordance.tsv

3 IDs in both files (1 only in file 1, 1 only in file 2)
gene  compared  full  partial  none  allele_concordance
   A         3     2        1     0               0.833
   B         2     1        1     0               0.750
✓ Concordance results saved: concordance.tsv

concordance.tsv has one row per ID found in both files, with A_file1, A_file2 and A_concordance for each gene. allele_concordance is the number of matching alleles divided by the number compared.

Option	Default	Description
`--file1_id`, `--file2_id`	(required)	ID column in each file
`--genes`, `-g`	every canonical gene in both files	Comma-separated genes
`--digits`, `-d`	all fields	Compare at no more than this resolution
`--output`, `-o`	print	Per-sample results file

The two files can use different ID columns, but the IDs themselves must match. To link exome and array IDs, load both through the catalog instead (How IDs link).

See c4.hla for every class and function.