<jupyter_start><jupyter_text>This notebook weaves together the haplotype scaffold at biallelic sites phased via shapeit2 with the multiallelic and other extra (e.g., non-PASS) variants phased via mvncall.<jupyter_code>%run setup.ipynb
# haplotype scaffold
callset_phased = phase1_ar31.callset_phased
gt_phased = allel.GenotypeDaskArray(callset_phased['2L/calldata/genotype'])
pos_phased = allel.SortedIndex(callset_phased['2L/variants/POS'])
pos_phased.shape, gt_phased.shape
# define region we're going to analyse
loc_region = pos_phased.locate_range(0, 6000000)
loc_region
# extract data for region, remove colony parents (8 samples) 
pos_phased_region = pos_phased[loc_region]
gt_phased_region = gt_phased[loc_region][:, :-8].compute()
pos_phased_region.shape, gt_phased_region.shape 
# load mvn haplotypes
callset_extras = np.load('../data/phasing_extra_phase1.mvncall.200.npz')
pos_extras = callset_extras['variants']['POS']
gt_extras = callset_extras['calldata']['genotype']
pos_extras.shape, gt_extras.shape
# concatenate
gt_combined = np.concatenate([gt_phased_region, gt_extras], axis=0)
pos_combined = np.concatenate([pos_phased_region, pos_extras], axis=0)

# sort by position
idx_sorted = np.argsort(pos_combined)
gt_combined = gt_combined[idx_sorted]
pos_combined = pos_combined[idx_sorted]
# obtain data from unphased callset - only needed for variant annotations
callset = phase1_ar31.callset
pos_all = allel.SortedIndex(callset['2L/variants/POS'])
ann_all = callset['2L/variants/ANN'][:][['Annotation', 'HGVS_p', 'HGVS_c']]
ann_all
# locate the intersection with unphased callset - needed to tie in annotations
loc1, _ = pos_all.locate_intersection(pos_combined)
np.count_nonzero(loc1)
# extract annotations for the phased variants
ann_combined = ann_all[loc1]
ann_combined
# save
haps_combined = allel.GenotypeArray(gt_combined).to_haplotypes()
np.savez_compressed('../data/haps_phase1.npz', haplotypes=haps_combined, POS=pos_combined, ANN=ann_combined)
haps_combined.nbytes
!ls -lh ../data/haps_phase1.npz<jupyter_output>-rw-rw-r-- 1 chris chris 12M Dec  1 14:26 ../data/haps_phase1.npz
{'fraction_non_alphanumeric': 0.10931558935361217, 'fraction_numerical': 0.019486692015209126, 'mean_word_length': 7.22265625, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '336256', 'license': 'MIT', 'n_tokens_mistral': 830, 'n_tokens_neox': 778, 'n_words': 152, 'path': 'notebooks/data_combined_haplotypes.ipynb', 'repo_name': 'shaghayegh-flower/ag1000g-phase1-vgsc-report', 'size': '2104'}