Multiallelic in multisample MatrixTable

Hello. What is the proper way to query multiallelic sites in a MatrixTable containing hundreds of samples? I know you can ask for hl.len(mt.alleles) > 3 , but what can you do for multiple samples?

Can you be specific about the query you want to write? You’re correct that:

mt.filter_rows(hl.len(mt.alleles) >= 3)

will return a matrixtable with only multiallelic variants.

Let’s say we have a row, for which we have hl.len(mt.alleleS) == 8. We also have 3 samples at this position each with the following GT 1/2 0/1 2/3. How can we query what alleles each sample has and what they are?

@matin ,

I directly answer your question below, but I think we can give you a better answer if you tell us what you’re ultimately trying to accomplish. Are you trying to run a GWAS on these alleles? Are you doing a burden test? Something else?

Hail MatrixTables are not indexed by number of alleles, so this pipeline necessary scans across the entire dataset. For large datasets, this will take a long time. If you plan to do a lot of work on this subset of the variants, you might want to filter once and write to a smaller file and use the smaller file for your analysis.

mt = hl.read_matrix_table(...)
mt = mt.filter_rows(hl.len(mt.alleles) >= 3)

Will print the top-left corner of the filtered matrix of genotypes to your Notebook or terminal.

If you want to export those genotypes to a sharded VCF you can replace with:

hl.export_vcf(mt, 'gs://path/to/file.vcf.bgz', parallel=True)

If you want to collect those genotypes as an array of arrays of genotypes, replace with:

mt = mt.annotate_rows(genotypes = hl.agg.collect(mt.GT))
array_of_arrays_of_gts = mt.genotypes.collect()

To elaborate on my previous comment: say we have three samples in our MatrixTable. It also has a row at locus 5:5000 for which we have an alleles field that contains 8 values ( e.g, ["CTGTGTGTGTGTG","<*>", "C", "CTGTG", "CTGTGTG", "CTGTGTGTG", "CTGTGTGTGTGTGT... ). We have called variants for each sample at this locus. Assume that the GT values are 1/2, 0/1, and 2/3, respectively for each sample. Now I want to query the ref/alt values for say, the first sample. my ultimate goal is to map alt/ref’s at this row to each sample.

Are you asking how to index into the alleles array? Meaning if you have these GTs with your shared alleles field:

1/2 0/1 2/3

You want these values:

("<*>", "C") ("CTGTGTGTGTGTG", "<*>") ("C", "CTGTG")


You can index both the alleles array and the call like normal Python arrays:

mt = mt.annotate_entries(
    first_allele_sequence = mt.alleles[mt.GT[0]],
    second_allele_sequence = mt.alleles[mt.GT[1]]