How do I show a (section of) an alignment for a phased genome

[hyanwong]

E.g. if I want to show the alignment for certain ranges (e.g. pos 100-120) for some of the sample genomes, like this:

    sample  alignment
    0       A..G......C..G.....AT
    0       T..G......A..C.....AG
    1       A..G......A..C.....AC
    1       A..C......C..G.....CT
            |         |         |
            100       110       120

I assume it could be useful to pick out snippets of sample genomes like this, e.g. for verification purposes.

[jeromekelleher]

[hyanwong]

Sorry, this is in sgkit, not tskit (so there is no alignments method). This is what I have at the moment, but it's probably not the recommended way to do it.

import sgkit
import numpy as np
ds = sgkit.simulate_genotype_call_dataset(n_variant=8, n_sample=3, missing_pct=0, phased=True, seed=123)

alleles = ds['variant_allele'].values.astype(str)
sites = np.arange(ds['call_genotype'].shape[0])
for sample, genome in [(2,0), (2,1)]:  # just pick the first 2 genomes of sample 2
    genotypes = ds['call_genotype'][:,sample, genome]
    print(f"sample {sample} ({genome}): " + "".join(alleles[sites, genotypes.values]))

[jeromekelleher]

Doh! I just deleted my comment.

[jeromekelleher]

Alignments aren't something we're considering at the moment, it's all about the variant matrix.

[hyanwong]

Right, but instead of an alignment, it would be helpful to get haplotypes out, I think? E.g. see my example above. Presumably people might want to know (a snippet of) the haplotype for a particular sample? E.g. for SARS-CoV-2 genomes

[jeromekelleher]

Sure - it's just not a focus currently