Prioritizing non-coding variants

README.md

@@ -1,5 +1,117 @@
 # Variant Interpretation Pipeline
 
+This is VIP with added non-coding functionality based on the [GREEN-DB](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8934622/) method. 
+It has extra annotations resources which includes:
+- TFBS regions (The TF that can bind the region)
+- UCNE regions
+- DNase regions
+- FATHMM-MKL scores
+- ReMM scores
+- ncER scores
+- GREEN-DB constraint values 
+- HPO phenotypes
+
+This branch has the annotations that are expected input for the [vip-decision-tree](https://github.com/molgenis/vip-decision-tree/tree/feat/annotation) score annotation branch. This score annotation tool calculates and annotates a variants score based on the added annotations. 
+
+These added annotations require resources that are not downloaded by the install.sh and need to be downloaded manually. 
+They can be downloaden from the [GREEN-DB download website](https://green-varan.readthedocs.io/en/latest/Download.html). The downloaded resources require some preprocessing follow the following steps for each annotation resource:
+
+[BCFtools](https://samtools.github.io/bcftools/) is required te perform the pre processing steps
+
+## TFBS/ncER
+
+
+File: {Build}_ncER_perc.bed.gz
+File: {Build}_TFBS_merged.bed.gz
+
+Unpack the downloaded gz files and follow these steps:
+
+```
+grep -v "#" file.bed  | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > file.bed.gz
+tabix -p bed file.bed.gz
+```
+
+## UCNE/DNase
+
+File: {Build}_UCNE.bed.gz
+File: {Build}_DNase_merged.bed.gz
+
+unpack the downloaded gz files and follow these steps:
+
+Perform this extra step for the **DNase** and **UCNE** file to add a fourth column.
+Unpack the file and perform this step:
+
+```
+awk 'BEGIN{ FS = OFS = "\t" } { print $1,$2,$3, (NR==1? "id" : "UCNE") }' {Build}_UCNE.bed > tmp && mv tmp GRCh38_UCNE.bed 
+
+awk 'BEGIN{ FS = OFS = "\t" } { print $1,$2,$3, (NR==1? "id" : "DNase") }' {Build}_DNase.merged.bed > tmp && mv tmp GRCh37_DNase.bed
+```
+
+And then for both files:
+```
+grep -v "#" file.bed  | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > file.bed.gz
+tabix -p bed file.bed.gz
+```
+
+## GREEN-DB constraint value
+
+File: {build}_GREEN-DB.bed.gz
+
+Unpack the downloaded gz files and follow these steps:
+Or use this [script](https://github.com/JonathanKlimp/Graduation-scripts/blob/main/data_conversion_scripts/process_bed_constraint_VEP.sh)
+
+First step is to remove all unused columns so that only the constraint value remains.
+Second step removes all values smaller than 0.7 and NA values.
+```
+cut -f 1,2,3,9 {Build}_GREEN-DB.bed > tmp.txt                          
+
+awk '$4 > 0.7 && $4 != "NA"' tmp.txt > {Build}_GREEN-DB-constraint.bed 
+
+grep -v "#" file.bed  | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > file.bed.gz
+tabix -p bed file.bed.gz
+```
+
+## ReMM
+
+File: {Build}_ReMM.tsv.gz
+
+Perform these steps:
+
+```
+gunzip {Build}_ReMM.tsv.gz
+
+awk 'BEGIN{ FS = OFS = "\t" } { print $1,$2,$2,$3 }' {Build}_ReMM.tsv > {Build}_ReMM_new_column.bed
+#grep -v "#" {Build}_ReMM_new_column.bed | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > {Build}_ReMM.bed.gz
+tabix -p bed {Build}_ReMM.bed.gz
+
+```
+
+## FATHMM-MKL
+
+File: {Build}_FATHMM_MKL.tsv.gz
+
+Use this script to convert the tsv file to vcf format: [conversion script](https://github.com/JonathanKlimp/Graduation-scripts/blob/main/data_conversion_scripts/convert_score_tsv_to_vcf.sh)
+
+## Phenotyepes
+
+File: GREEN-DB_v2.5.db.gz
+
+Extract the HPO phenotypes table and the table with the regions from the DB file using [SQLite](https://sqlite.org/index.html) or another preferred DB browser.
+
+Once extracted perform the following:
+This command checks if an ID is matched in both files. If this is the case the chrom, pos, stop, phenotype and phenotypeID are printed to a new file.
+
+```
+awk -v FS='\t' -v OFS='\t' 'FNR==NR{a[$4]=$1 FS $2 FS $3;next} ($1 in a) {print a[$1],$3,$4}' GREEN-DB_{Build}_regions.csv GREEN-DB_phenotypes.csv > GREEN_DB_{Build}_region_phenos.bed
+```
+
+Then perform the following command:
+
+```
+grep -v "#" GREEN_DB_{Build}_region_phenos.bed | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > GREEN_DB_{Build}_region_phenos.bed.gz
+tabix -p bed GREEN_DB_{Build}_region_phenos.bed
+```
+
 ## Requirements
 - POSIX compatible system (Linux, OS X, etc) / Windows through [WSL](https://en.wikipedia.org/wiki/Windows_Subsystem_for_Linux)
 - Bash 3.2 (or later)

modules/annotate.nf

@@ -18,6 +18,14 @@ process annotate {
     capiceInputPath = "${id}_chunk${order}_capice_input.tsv"
     capiceOutputPath = "${id}_chunk${order}_capice_output.tsv.gz"
     capiceModelPath = params[params.assembly + "_annotate_capice_model"]
+    greenDbConstraintPath = params[params.assembly + "_constraint_GREEN_DB"]
+    dnaseRegionsPath = params[params.assembly + "_DNase_regions"]
+    tfbsRegionsPath = params[params.assembly + "_TFBS_regions"]
+    ucneRegionsPath = params[params.assembly + "_UCNE_regions"]
+    fathmmMKLScoresPath = params[params.assembly + "_fathmm_MKL_scores"]
+    ncErScoresPath = params[params.assembly + "_ncER_scores"]
+    reMMScoresPath = params[params.assembly + "_ReMM_scores"]
+    regionPhenosPath = params[params.assembly + "_region_phenotypes"]
 
     template 'annotate.sh'
 }

nextflow.config

@@ -42,6 +42,14 @@ params {
   GRCh37_annotate_vep_plugin_vkgl = "${projectDir}/resources/GRCh37/vkgl_consensus_20211201.tsv"
   GRCh37_annotate_capice_model = "${projectDir}/resources/GRCh37/capice_model_v4.0.0-v2.pickle.dat"
   GRCh37_report_genes = "${projectDir}/resources/GRCh37/GCF_000001405.25_GRCh37.p13_genomic_g1k.gff.gz"
+  GRCh37_constraint_GREEN_DB = "${projectDir}/resources/GRCh37/GRCh37_GREEN_DB_constraint.bed.gz"
+  GRCh37_DNase_regions = "${projectDir}/resources/GRCh37/GRCh37_DNase.bed.gz"
+  GRCh37_TFBS_regions = "${projectDir}/resources/GRCh37/GRCh37_TFBS.bed.gz"
+  GRCh37_UCNE_regions = "${projectDir}/resources/GRCh37/GRCh37_UCNE.bed.gz"
+  GRCh37_fathmm_MKL_scores = "${projectDir}/resources/GRCh37/GRCh37_FATHMM-MKL_NC.vcf.gz"
+  GRCh37_ncER_scores = "${projectDir}/resources/GRCh37/GRCh37_ncER_perc.bed.gz"
+  GRCh37_ReMM_scores = "${projectDir}/resources/GRCh37/GRCh37_ReMM.bed.gz"
+  GRCh37_region_phenotypes = "${projectDir}/resources/GRCh37/GREEN_DB_GRCh37_region_phenos.bed.gz"
 
   GRCh38_reference = "${projectDir}/resources/GRCh38/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz"
   GRCh38_annotate_vep_custom_gnomad = "${projectDir}/resources/GRCh38/gnomad.genomes.v3.1.2.sites.stripped.vcf.gz"
@@ -54,6 +62,14 @@ params {
   GRCh38_annotate_vep_plugin_vkgl = "${projectDir}/resources/GRCh38/vkgl_consensus_20211201.tsv"
   GRCh38_annotate_capice_model = "${projectDir}/resources/GRCh38/capice_model_v4.0.0-v2.pickle.dat"
   GRCh38_report_genes = "${projectDir}/resources/GRCh38/GCF_000001405.39_GRCh38.p13_genomic_mapped.gff.gz"
+  GRCh38_constraint_GREEN_DB = "${projectDir}/resources/GRCh38/GRCh38_GREEN_DB_constraint.bed.gz"
+  GRCh38_DNase_regions = "${projectDir}/resources/GRCh38/GRCh38_DNase.bed.gz"
+  GRCh38_TFBS_regions = "${projectDir}/resources/GRCh38/GRCh38_TFBS.bed.gz"
+  GRCh38_UCNE_regions = "${projectDir}/resources/GRCh38/GRCh38_UCNE.bed.gz"
+  GRCh38_fathmm_MKL_scores = "${projectDir}/resources/GRCh38/GRCh38_FATHMM-MKL_NC.vcf.gz"
+  GRCh38_ncER_scores = "${projectDir}/resources/GRCh38/GRCh38_ncER_perc.bed.gz"
+  GRCh38_ReMM_scores = "${projectDir}/resources/GRCh38/GRCh38_ReMM.bed.gz"
+  GRCh38_region_phenotypes = "${projectDir}/resources/GRCh38/GREEN_DB_GRCh38_region_phenos.bed.gz"
 
   singularity_image_dir = "${projectDir}/images"
   chunk_size = 10000

resources/vep/plugins/VIPVaranLevel.pm

@@ -0,0 +1,196 @@
+package VIPVaranLevel;
+
+use strict;
+use warnings;
+
+use Bio::EnsEMBL::Variation::Utils::BaseVepPlugin;
+use base qw(Bio::EnsEMBL::Variation::Utils::BaseVepPlugin);
+
+=head1 NAME
+ VIPVaranLevel
+=head1 SYNOPSIS
+ mv VIPVaranLevel.pm ~/.vep/Plugins
+=head1 DESCRIPTION
+ Plugin to annotate VIPVaran level per variant based on constraint level, fathmm, ncER, ReMM scores and region.
+=cut
+
+sub version {
+    return '0.1';
+}
+
+# Feature includes regulatory features 
+sub feature_types {
+    return [ 'Transcript', 'Intergenic', 'Feature'];
+}
+
+sub get_header_info {
+    return {
+        VIPVaranLevel => "Score between 1 and 3 based on constraint level, fathmm, ncER, ReMM scores and region"
+    };
+}
+
+my %min_scores = (	
+  fathmm => 0.5,
+  ncER => 0.5,
+  ReMM => 0.499,
+  constraint => 0.7
+);
+
+sub is_in_region {
+    my $region_values = $_[0];
+
+    ## if region contains value
+    unless($region_values eq "") {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+sub tool_min_score {
+    my $fathmm_score = $_[0];
+    my $ncER_score = $_[1];
+    my $ReMM_score = $_[2];
+
+    # add logic for when there are multiple scores for the same variant from the same tool.
+    # example: 99.7852&99.7217  (can be more than 2) ReMM: 0.1710&0.9490&0.9560 (low and high score what to do?)
+    if($fathmm_score >= $min_scores{"fathmm"} || $ncER_score >= $min_scores{"ncER"} || $ReMM_score >= $min_scores{"ReMM"}) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+sub constraint_min_score {
+    my $constraint_score = $_[0];
+
+    if ($constraint_score >= $min_scores{"constraint"}) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+sub run {
+    my ($self, $transcript_variation_allele, $line_hash) = @_;
+
+    my $base_variation_feature = $transcript_variation_allele->base_variation_feature;
+    my @vcf_line = @{$base_variation_feature->{_line}};
+
+    #my $data = @{$self->get_data}; get_data is not a thing
+
+    # my $test3 = $self->ncER; werkt ook niet met $line_hash
+    #my @test_data = @{$self->green_db_tool_scores};
+    # code to write to file
+    my $filename = '/groups/solve-rd/tmp10/projects/vip/feat/non-coding/test/VIPVaranLevel.log';
+    open(my $file, '>>', $filename) or die $!;
+    #print($file "HIER onder is vcf line");
+    # foreach (@vcf_line) {
+    #     print($file "$_\n"); 
+    # }
+    print($file "params\n");
+    my @params_array = $self->{"params"};
+    foreach (@params_array) {
+        print($file "$_\n");
+        foreach (@$_) {
+            print($file "$_\n");
+        }
+    }
+    print($file "variant feature types\n");
+    my @vft_array = $self->{"variant_feature_types"};
+    foreach (@vft_array) {
+        print($file "$_\n");
+        foreach (@$_) {
+            print($file "$_\n");
+        }        
+    }
+    print($file "feature types\n");
+    my @f_types_array = $self->{"feature_types"};
+    foreach (@f_types_array) {
+        print($file "$_\n");
+        foreach (@$_) {
+            print($file "$_\n");
+        }        
+    }
+    print($file "variant feature types wanted\n");
+    foreach my $variant_feature_wanted ($self->{"variant_feature_types_wanted"}) {
+        foreach my $key (keys %$variant_feature_wanted) {
+                print($file "in keys loop variant feature types wanted\n");
+                print($file $key);
+                print($file "=>");
+                print($file $variant_feature_wanted->{$key});
+                print($file "\n")
+        }
+    }
+    print($file "feature types wanted\n");
+    foreach my $feature_wanted ($self->{"feature_types_wanted"}) {
+        foreach my $key (keys %$feature_wanted) {
+                print($file "in keys loop feature types wanted\n");
+                print($file $key);
+                print($file "=>");
+                print($file $feature_wanted->{$key});
+                print($file "\n")
+        }
+    }
+    print($file "config\n");
+    foreach my $config ($self->{"config"}) {
+        my @config_custom_array = $self->{"plugin"};
+        foreach (@config_custom_array) {
+            print($file "$_\n");
+        }
+        my @pick_order = $self->{"pick_order"};
+        foreach (@pick_order) {
+            print($file "$_\n");
+        }
+        foreach my $key (keys %$config) {
+                print($file "in keys loop config\n");
+                print($file $key);
+                print($file "=>");
+                print($file $config->{$key});
+                print($file "\n")
+        }
+    }
+    #print($file @vcf_line); # bevat chrom pos ref alt, 0 en 3x "." 
+    # print($file $transcript_variation_allele); # is een hash
+    close($file);
+    # score is 0 by default
+    my $score = 0;
+
+    # $region = something...
+    # $fathmm = something..
+    # $ReMM = something...
+    # $ncER = something...
+    # $constraint = something...
+    my $region = 0;
+    my $fathmm_score = 0;
+    my $ncER_score = 0; 
+    my $ReMM_score = 0;
+    my $constraint_score = 0;
+
+    # Higher level = higher chance of pathogenicity 
+    # If variant lays in in TFBS,DNase or UCNE = level 1
+    # If variant has score of FATHMM above 0.5, ncER above 0.5 or ReMM above 0.499 and above = level 2
+    # If variant has constraint value above 0.7 and above = level 3
+    # Something with phenotype for level 4?
+    # What if scores are sufficient but it is not in one of the regions?
+    # Something with UTRAnnotator?
+
+    if(is_in_region($region) && tool_min_score($fathmm_score, $ncER_score, $ReMM_score) && constraint_min_score($constraint_score)) {
+        $score = 3;
+    } elsif(is_in_region($region) && tool_min_score($fathmm_score, $ncER_score, $ReMM_score) && !constraint_min_score($constraint_score)) {
+        $score = 2;
+    } elsif(is_in_region($region) && !tool_min_score($fathmm_score, $ncER_score, $ReMM_score) && !constraint_min_score($constraint_score)) {
+        $score = 1;
+    } else {
+        $score = 0;
+    }
+
+    # return {
+    #     #VIPVaranLevel => $results
+    #     VIPVaranLevel => $score
+    # };
+    #return $score;
+    return {};
+}
+
+1;

templates/annotate.sh

@@ -186,6 +186,15 @@ vep() {
   args+=("--plugin" "UTRannotator,!{vepPluginUtrAnnotatorPath}")
   args+=("--custom" "!{vepCustomPhyloPPath},phyloP,bigwig,exact,0")
 
+  args+=("--custom" "!{greenDbConstraintPath},constraint,bed,overlap,0")
+  args+=("--custom" "!{dnaseRegionsPath},region,bed,overlap,0")
+  args+=("--custom" "!{tfbsRegionsPath},region,bed,overlap,0")
+  args+=("--custom" "!{ucneRegionsPath},region,bed,overlap,0")
+  args+=("--custom" "!{fathmmMKLScoresPath},fathmm,vcf,exact,0,fathmm")
+  args+=("--custom" "!{ncErScoresPath},ncER,bed,overlap,0")
+  args+=("--custom" "!{reMMScoresPath},ReMM,bed,overlap,0")
+  args+=("--custom" "!{regionPhenosPath},phenotype,bed,overlap,0")
+
   if [ -n "!{vepPluginArtefact}" ]; then
     args+=("--plugin" "Artefact,!{vepPluginArtefact}")
   fi