<jupyter_start><jupyter_text>Download VCFs<jupyter_code>!mkdir -pv data/vcf
%%bash
cd data/vcf
for c in {1..22}
do
    wget --no-clobber ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/ALL.chr${c}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
done<jupyter_output>File ‘ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr2.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr3.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr4.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr5.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr7.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr8.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retrieving.
File ‘ALL.chr9.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz’ already there; not retr[...]<jupyter_text>Find max number of ALT alleles<jupyter_code>import sys
import numpy as np
import dask
from dask import delayed
from dask import distributed
import allel
allel.__version__
cluster = distributed.LocalCluster(n_workers=4, threads_per_worker=1)
client = distributed.Client(cluster)
client
def find_max_alt(chrom):
    import allel
    vcf_path = 'data/vcf/ALL.chr{}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'.format(chrom)
    callset = allel.read_vcf(vcf_path, fields=['variants/numalt'])
    return callset['variants/numalt'].max()
# results = [delayed(find_max_alt)(chrom) for chrom in range(1, 23)]
# max_alt = dask.compute(results)
# max_alt
# [12, 7, 8, 6, 7, 9, 8, 8, 7, 7, 9, 6, 7, 7, 8, 7, 6, 7, 7, 6, 6, 8]<jupyter_output><empty_output><jupyter_text>Parse VCF to Zarr<jupyter_code>def parse_vcf_to_zarr(chrom):
    import allel
    import zarr
    import numcodecs
    compressor = numcodecs.Blosc(cname='zstd', clevel=1, shuffle=-1)
    vcf_path = 'data/vcf/ALL.chr{}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'.format(chrom)
    zarr_path = 'data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes'
    root = zarr.open(zarr_path, mode='a')
    if chrom in root:
        if root[chrom].attrs.get('complete', False):
            print('already done, skipping')
            return
        else:
            del root[chrom]
    allel.vcf_to_zarr(vcf_path, zarr_path, group=chrom, fields='*', alt_number=12,
                      compressor=compressor, log=sys.stdout, chunk_length=2**18,
                      chunk_width=128)
    root[chrom].attrs['complete'] = True

tasks = [delayed(parse_vcf_to_zarr)(chrom) for chrom in range(1, 23)]
dask.compute(tasks)<jupyter_output><empty_output><jupyter_text>Add sample data<jupyter_code>import pandas as pd
df_samples = pd.read_csv('data/vcf/integrated_call_samples_v3.20130502.ALL.panel', 
                         sep='\t', 
                         usecols=['sample', 'pop', 'super_pop', 'gender'])
df_samples
import zarr
zarr_path = 'data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes'
callset = zarr.open(zarr_path, mode='a')
print(callset.tree(level=1))
print(callset['1'].tree())
# Check samples from VCF same as panel
samples_from_vcf = callset['1/samples'][:].tolist()
samples_from_vcf[:5]
samples_from_panel = df_samples['sample'].tolist()
samples_from_panel[:5]
samples_from_vcf == samples_from_panel
callset.require_group('samples')
import numcodecs
callset.create_dataset('samples/ID', data=df_samples['sample'], compressor=None, chunks=None, object_codec=numcodecs.VLenUTF8())
callset.create_dataset('samples/pop', data=df_samples['pop'], compressor=None, chunks=None, object_codec=numcodecs.VLenUTF8())
callset.create_dataset('samples/super_pop', data=df_samples['super_pop'], compressor=None, chunks=None, object_codec=numcodecs.VLenUTF8())
callset.create_dataset('samples/gender', data=df_samples['gender'], compressor=None, chunks=None, object_codec=numcodecs.VLenUTF8())
print(callset.tree(level=1))
print(callset['samples'].tree())
callset['samples/ID'][:]<jupyter_output><empty_output><jupyter_text>Consolidate metadata<jupyter_code>zarr.consolidate_metadata(zarr_path)
print(zarr.open_consolidated(zarr_path).tree(level=1))<jupyter_output>/
 ├── 1
 ├── 10
 ├── 11
 ├── 12
 ├── 13
 ├── 14
 ├── 15
 ├── 16
 ├── 17
 ├── 18
 ├── 19
 ├── 2
 ├── 20
 ├── 21
 ├── 22
 ├── 3
 ├── 4
 ├── 5
 ├── 6
 ├── 7
 ├── 8
 ├── 9
 └── samples<jupyter_text>Upload to GCS<jupyter_code>!gsutil -m rsync -r data/zarr/ gs://1000genomes-zarr/<jupyter_output>Building synchronization state...
At source listing 10000...
Starting synchronization...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.11.0 [Content-Type=application/octet-stream]...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.12.0 [Content-Type=application/octet-stream]...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.13.0 [Content-Type=application/octet-stream]...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.14.0 [Content-Type=application/octet-stream]...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.15.0 [Content-Type=application/octet-stream]...
Copying file://data/zarr/ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes/1/calldata/GT/0.18.0 [Content-Type=application/octet-stream]...
Copyi[...]
{'fraction_non_alphanumeric': 0.1390818228019651, 'fraction_numerical': 0.0642046417076063, 'mean_word_length': 5.581939799331104, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '391840', 'license': 'MIT', 'n_tokens_mistral': 2547, 'n_tokens_neox': 2266, 'n_words': 340, 'path': 'data/1000genomes-phase3/vcf-to-zarr.ipynb', 'repo_name': 'alimanfoo/sgkit-requirements', 'size': '5903'}