<jupyter_start><jupyter_text>Example Notebook This notebook is an example that shows how to use [VariantSpark](https://github.com/aehrc/VariantSpark) with the [Hail v0.2 library](https://hail.is/) and compares the results with PCA and logistic regression. For demonstration purposes this notebook uses the sample dataset available in [ViGWAS](https://github.com/aehrc/VIGWAS). [We're always looking for suggestions and feedback. Please click here for a 1 minute survey](https://docs.google.com/forms/d/e/1FAIpQLScWoazw3-rgNFrZ5vcHL9JUmO0AX6Ji2P54Z2jNJZ-RAObuPg/viewform?usp=sf_link)  User Block<jupyter_code>S3_Path='s3://REPLACE_WITH_NOTEBOOK_PATH/'

## Some configs
numCPU = 32
memory = '60G'
numPartitions = numCPU*4<jupyter_output><empty_output><jupyter_text>Environment initialization<jupyter_code>## Environment init

import os
from pyspark import SparkContext
sc = SparkContext()

import hail as hl
import varspark.hail as vshl
vshl.init(sc=sc)
from bokeh.io import show, output_notebook
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource, FactorRange, LabelSet, Label
from bokeh.transform import factor_cmap
from bokeh.palettes import d3
from bokeh.core.properties import value
from bokeh.embed import file_html
from bokeh.resources import CDN
from bokeh.layouts import gridplot
from bokeh.models.mappers import CategoricalColorMapper

from pprint import pprint
output_notebook()

import re
import numpy as np
import math as math
import sys
import operator
from collections import OrderedDict
import subprocess
from itertools import cycle
import shutil<jupyter_output><empty_output><jupyter_text>Load VCF files<jupyter_code>mt = hl.import_vcf(path=S3_Path+'sample_input/V1.vcf.bgz',
                   skip_invalid_loci=True,
                   min_partitions=int(numPartitions))<jupyter_output><empty_output><jupyter_text>Sample Annotation Data Analysis<jupyter_code>Annot = hl.import_table(S3_Path+'sample_input/hipster.csv',
                        impute=True, delimiter=',').key_by('Sample')<jupyter_output><empty_output><jupyter_text>Annotate dataset with sample annotation<jupyter_code>mt = mt.annotate_cols(pheno = Annot[mt.s])<jupyter_output><empty_output><jupyter_text>PCA analysis<jupyter_code>eigenvalues, pcs, loadings = hl.hwe_normalized_pca(mt.GT, k=2)
mt = mt.annotate_cols(pcs = pcs[mt.s].scores)
p = hl.plot.scatter(pcs.scores[0], pcs.scores[1],
                    label=mt.cols()[pcs.s].pheno.Hipster,
                    title='PCA Case/Control', xlabel='PC1', ylabel='PC2', collect_all=True)
show(p)<jupyter_output><empty_output><jupyter_text>Logistic Regression<jupyter_code>covariates = [mt.pheno.isFemale, mt.pcs[0], mt.pcs[1]]

result = hl.logistic_regression_rows(test ='wald', 
                                          y=mt.pheno.isCase,
                                          x=mt.GT.n_alt_alleles(),
                                          covariates=covariates)

mt = mt.annotate_rows( logreg = result[mt.locus, mt.alleles])
p = hl.plot.manhattan(result.p_value)
show(p)<jupyter_output><empty_output><jupyter_text>Variant-Spark RandomForest<jupyter_code>rf_model = vshl.random_forest_model(y=mt.pheno.isCase, x=mt.GT.n_alt_alleles(),
                                    seed = 13, mtry_fraction = 0.1,
                                    min_node_size = 10, max_depth = 15)

rf_model.fit_trees(n_trees=100, batch_size=25)

impTable = rf_model.variable_importance()

mt = mt.annotate_rows(vs_score = impTable[mt.locus, mt.alleles].importance)
mt = mt.annotate_rows(vs_stats = mt.aggregate_rows(hl.agg.stats(mt['vs_score'])))
mt = mt.annotate_rows(z_score = (mt['vs_score'] - mt.vs_stats.mean)/mt.vs_stats.stdev)
mt = mt.annotate_rows(vs_score_converted = 10** -mt.z_score)
title = 'Variant-Spark Manhattan plot'
hover_fields = {'rsid': mt.rsid, 'vs_score': mt.vs_score}
p = hl.plot.manhattan(pvals=mt.vs_score_converted, hover_fields=hover_fields, title=title)
p.yaxis.axis_label = 'Z score of importantce score by VS'
show(p)<jupyter_output><empty_output><jupyter_text>Describe matrix Table<jupyter_code>mt.describe()<jupyter_output><empty_output><jupyter_text>Write MT to S3<jupyter_code>mt.write(S3_Path+'sample_output/dataset.mt',overwrite=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1089924670433145, 'fraction_numerical': 0.011299435028248588, 'mean_word_length': 4.726415094339623, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '426266', 'license': 'MIT', 'n_tokens_mistral': 1565, 'n_tokens_neox': 1505, 'n_words': 269, 'path': 'databricks/notebooks/VariantSpark_example_old.ipynb', 'repo_name': 'krisbock/variant-databricks', 'size': '4248'}