<jupyter_start><jupyter_text>Kipoi usage tutorial: Per-Variant PredictionsThis tutorial shows how to apply the model to a vcf file, if you would like to predict the *individual* effect of all variants in a vcf file.For this you need a vcf with all the variants of interest.<jupyter_code># Imports
import kipoi
from kipoi.pipeline import Pipeline
from kipoiseq.dataloaders import SingleVariantUTRDataLoader
import os
import numpy as np
import pandas as pd<jupyter_output><empty_output><jupyter_text>Load the Model<jupyter_code># Source model directly from directory
model = kipoi.get_model("../Framepool", source="dir")<jupyter_output>Using downloaded and verified file: /data/nasif12/home_if12/karollus/projects/models/Framepool/downloaded/model_files/weights/d1e9656725e730d509a09d5371e51bd2
WARNING:tensorflow:From /opt/modules/i12g/anaconda/envs/kipoi-Framepool2/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.
Instructions for updating:
Colocations handled automatically by placer.
WARNING:tensorflow:From /opt/modules/i12g/anaconda/envs/kipoi-Framepool2/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py:3445: calling dropout (from tensorflow.python.ops.nn_ops) with keep_prob is deprecated and will be removed in a future version.
Instructions for updating:
Please use `rate` instead of `keep_prob`. Rate should be set to `rate = 1 - keep_prob`.
WARNING:tensorflow:From /opt/modules/i12g/anaconda/envs/kipoi-Framepool2/lib/python3.6/site-packages/tensorflow/python/ops/math_o[...]<jupyter_text>Optional: download example files and hg19 fasta<jupyter_code>import urllib.request
import gzip
import shutil
# make ExampleFile directory if it does not exist
if not os.path.exists("ExampleFiles"):
    os.makedirs("ExampleFiles")
# Download vcf
urllib.request.urlretrieve("https://zenodo.org/record/3584238/files/patho.vcf.gz?download=1", 'ExampleFiles/patho.vcf.gz')
# Download vcf tabix
urllib.request.urlretrieve("https://zenodo.org/record/3584238/files/patho.vcf.gz.tbi?download=1", 'ExampleFiles/patho.vcf.gz.tbi')
# Download GTF
urllib.request.urlretrieve("https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/genes/hg19.ensGene.gtf.gz", 'ExampleFiles/hg19.gtf.gz')
# Download id mapping file
urllib.request.urlretrieve("https://zenodo.org/record/3584238/files/hg19_idmap.tsv?download=1", 'ExampleFiles/hg19_idmap.tsv')
# Download gzipped hg19 fasta (warning: 900mb)
urllib.request.urlretrieve("https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz", 'ExampleFiles/hg19.fa.gz')
# unzip
with gzip.open('ExampleFiles/hg19.fa.gz', 'rb') as f_in:
    with open('ExampleFiles/hg19.fa', 'wb') as f_out:
        shutil.copyfileobj(f_in, f_out)<jupyter_output><empty_output><jupyter_text>Provide the Parameters This Dataloader requires the following input files:1. GTF which provides 5'UTR as a feature type (e.g. 5UTR in gencode, five_prime_utr in Ensembl)2. fasta file that provides the reference genome3. bgzip compressed (single sample) vcf that provides the variants   The vcf must be sorted (by position) and a tabix index must be present (must lie in the same directory and have the same name + .tbi)All files must agree on whether there is a chr prefix (e.g. is it chr1 or just 1). **NB: Some common issues**1. Files do not agree with respect to the ordering or the chr prefix. This creates all kinds of problems.2. The vcf includes non-DNA characters, e.g. * or /. The model cannot currently account for this<jupyter_code># Path of the vcf file
vcf_path = "ExampleFiles/patho.vcf.gz"

# Path of the fasta file
fasta_path = "ExampleFiles/hg19.fa"

# Path of the gtf file
gtf_path = "ExampleFiles/hg19.gtf.gz"
# How are the 5'UTR called in your gtf
feature_type = "5UTR"

# output file path
output_file_path = "patho.tsv"<jupyter_output><empty_output><jupyter_text>Filter your gtf (if necessary) If your vcf is small, or only covers small parts of the human genome, it is advisable to first subset your gtf to the corresponding regions to speed up the dataloader.Note: The dataloader was designed to handle very large vcf files (eg Gnomad VCF) that do not fit into memory easily. Thus, for very large VCF, filtering the gtf will not be beneficial.<jupyter_code>import pyranges as pr
from cyvcf2 import VCF
# Import gtf with pyranges
gr = pr.read_gtf(gtf_path)
id_set = set()
for var in VCF(vcf_path):
    chrom = var.CHROM
    pos = var.POS
    id_set = set(gr[chrom, pos-1:pos].df.gene_id) | id_set # vcf is 1-based, pyranges is not
gr_subset = gr[gr.gene_id.isin(id_set)]
gtf_path = "ExampleFiles/reduced.gtf"
gr_subset.to_gtf(gtf_path)<jupyter_output><empty_output><jupyter_text>Run Prediction<jupyter_code>pipeline = Pipeline(model, SingleVariantUTRDataLoader)
pipeline.predict_to_file(output_file_path, {"gtf_file":gtf_path, 
                               "fasta_file":fasta_path,
                               "vcf_file":vcf_path,
                                "feature_type":feature_type},
                              batch_size=64);<jupyter_output>1it [00:01,  1.10s/it]<jupyter_text>Load results<jupyter_code># Load data as dataframe
df = pd.read_csv(output_file_path, sep="\t")
df = df.rename(index=str, columns={"metadata/variant/chr":"chr",
          "metadata/exon_positions":"exon_positions",
          "metadata/transcript_id":"transcript_id",
          "metadata/variant/str":"variants",
          "metadata/variant/ref":"ref",
          "metadata/variant/alt":"alt",
          "metadata/variant/chrom":"chrom",
          "metadata/variant/id":"id",  
          "metadata/variant/pos":"pos", 
          "preds/mrl_fold_change":"mrl_fold_change",
          "preds/shift_1":"shift_1",
          "preds/shift_2":"shift_2"}
)<jupyter_output><empty_output><jupyter_text>Optional: id map to get a clearer output<jupyter_code># id map path
# Provide an id mapping file to get a richer output:
id_map_path = "ExampleFiles/hg19_idmap.tsv"

# Id map
df_map = pd.read_csv(id_map_path, sep="\t")
df = df.merge(df_map, on="transcript_id")
df<jupyter_output><empty_output><jupyter_text>Filter results1. Choose strongest result per gene (requires id mapping)2. Reduce to absolute effect > 0.5<jupyter_code>df["abs_mrl_fc"] = np.abs(df["mrl_fold_change"])
idx = df.groupby(['gene_name'])['abs_mrl_fc'].transform(max) == df['abs_mrl_fc']
df_max = df[idx]
df_max[np.abs(df_max.mrl_fold_change) > 0.5]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10940381069452981, 'fraction_numerical': 0.02166564228641672, 'mean_word_length': 4.806422836752899, 'pattern_counts': {'":': 16, '<': 35, '<?xml version=': 0, '>': 37, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '427667', 'license': 'MIT', 'n_tokens_mistral': 2346, 'n_tokens_neox': 2219, 'n_words': 551, 'path': 'Framepool/usage_example_single_effect.ipynb', 'repo_name': 'bernardo-de-almeida/models', 'size': '6508'}