=head1 LICENSE Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute Copyright [2016-2024] EMBL-European Bioinformatics Institute Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. =head1 CONTACT Xiaolei Zhang Ensembl =cut =head1 NAME UTRAnnotator =head1 SYNOPSIS mv UTRAnnotator.pm ~/.vep/Plugins vep -i variations.vcf --plugin UTRAnnotator,file=/path/to/uORF_starts_ends_GRCh38_PUBLIC.txt # skip annotation for variants with a 80% or higher overlap of the UTR vep -i variations.vcf --plugin UTRAnnotator,file=/path/to/uORF_starts_ends_GRCh38_PUBLIC.txt,max_overlap=80 =head1 DESCRIPTION A VEP plugin that annotates the effect of 5' UTR variant especially for variant creating/disrupting upstream ORFs. Available for both GRCh37 and GRCh38. Options are passed to the plugin as key=value pairs: file : (Required) Path to UTRAnnotator data file: - Download 'uORF_5UTR_GRCh37_PUBLIC.txt' or 'uORF_5UTR_GRCh38_PUBLIC.txt' from https://github.com/Ensembl/UTRannotator - Download from http://sorfs.org max_overlap : (Optional) Maximum percentage of overlap between variant and UTR for UTR annotation (default: 100) Citation About the role of 5'UTR variants in human genetic disease: Whiffin, N., Karczewski, K.J., Zhang, X. et al. Characterising the loss-of-function impact of 5’ untranslated region variants in 15,708 individuals. Nat Commun 11, 2523 (2020). https://doi.org/10.1038/s41467-019-10717-9 About UTRAnnotator: The original UTRAnnotator plugin is written by Xiaolei Zhang et al. Later adopted by Ensembl VEP plugins with some changes. You can find the original plugin here - https://github.com/ImperialCardioGenetics/UTRannotator Please cite the UTRannotator publication alongside the Ensembl VEP if you use this resource - Annotating high-impact 5'untranslated region variants with the UTRannotator Zhang, X., Wakeling, M.N., Ware, J.S, Whiffin, N. Bioinformatics; doi: https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btaa783/5905476 =cut package UTRAnnotator; use base qw(Bio::EnsEMBL::Variation::Utils::BaseVepPlugin); use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp); use List::Util qw(min max); use Scalar::Util qw(looks_like_number); use strict; sub feature_types { return ['Transcript']; } sub variant_feature_types { return ['VariationFeature', 'StructuralVariationFeature']; } sub new { my $class = shift; my $self = $class->SUPER::new(@_); printf "\nWarning: no FASTA file specified (--fasta). VEP running will take a long time." if ($self->config->{fasta} eq ""); my $param_hash = $self->params_to_hash(); if (-e $param_hash->{file}){ open my $fh, "<", $param_hash->{file} or die $!; my %uORF_evidence; while (<$fh>) { chomp; my ($chr, $pos, $gene, $strand, $type, $stop_pos) = split /\t/; next if $chr eq "chr" && $pos eq "pos"; die "$chr does not exist; Check your reference file." unless ($chr =~ m/^(chr*)/); die "Position $pos is not a number; Check your reference file." unless looks_like_number($pos); my $key = $chr . ":" . $pos; # chr has 'chr' prefix $uORF_evidence{$key} = 1; } close $fh; $self->{uORF_evidence} = \%uORF_evidence; } else { printf "Warning: small ORF file not found. For human, you could use the curated list of uORFs found in the repository (https://github.com/Ensembl/UTRannotator):\n" . "'uORF_starts_ends_GRCh37_PUBLIC.txt' for GRCh37 or 'uORF_starts_ends_GRCh38_PUBLIC.txt' for GRCh38.\n"; } # maximum percentage UTR overlap $self->{max_overlap} = defined $param_hash->{max_overlap} ? $param_hash->{max_overlap} : 100; # Kozak-Strength Class my %kozak_strength; $kozak_strength{1}='Weak'; $kozak_strength{2}='Moderate'; $kozak_strength{3}='Strong'; $self->{kozak_strength} = \%kozak_strength; return $self; } sub get_header_info { my $self->{_header_info} = { '5UTR_consequence' => 'Variant consequence from UTRAnnotator', '5UTR_annotation' => 'Variant annotation from UTRAnnotator', Existing_uORFs => 'The number of existing uORFs with a stop codon within the 5 prime UTR', Existing_OutOfFrame_oORFs => 'The number of existing out-of-frame overlapping ORFs (OutOfFrame oORF) at the 5 prime UTR', Existing_InFrame_oORFs => 'The number of existing inFrame overlapping ORFs (inFrame oORF) at the 5 prime UTR', }; return $self->{_header_info}; } sub run { my ($self, $tva) = @_; #only annotate the effect if the variant is 5_prime_UTR_variant return {} unless grep {$_->SO_term eq '5_prime_UTR_variant'} @{$tva->get_all_OverlapConsequences}; #retrieve the variant info my $vf = $tva->can('variation_feature') ? $tva->variation_feature : $tva->structural_variation_feature; my $chr = ($vf->{chr} || $vf->seq_region_name); my $pos = ($vf->{start} || $vf->seq_region_start); my $pos_end = ($vf->{end} || $vf->seq_region_end); my @alleles = split /\//, $vf->allele_string; my $ref = shift @alleles; my $alt = $tva->can('variation_feature_seq') ? $tva->variation_feature_seq : ''; my %variant = ( "chr" => $chr, "pos" => $pos, "ref" => $ref, "alt" => $alt, ); #retrieve the UTR info: transcript id, strand, five prime UTR sequence, start and end genomic coordinates. my $t = $tva->transcript; my $transcript_id = (defined $t? $t->stable_id: undef); #retrieve the gene symbol of the transcript my $symbol = $t->{_gene_symbol} || $t->{_gene_hgnc}; #retrieve the strand of the transcript my $tr_strand = $t->strand + 0; my $cds = $t->translateable_seq(); #retrieve the five prime utr sequence my $five_prime_seq = (defined $t->five_prime_utr? $t->five_prime_utr->seq(): undef); #Type: five_prime_feature - Bio::EnsEMBL::Feature my $UTRs = $t->get_all_five_prime_UTRs(); my @five_utr_starts; my @five_utr_ends; foreach my $utr (@$UTRs){ my $UTR_start = $utr->start(); # this will return the absolute starting positions in chromosome of the UTR exons my $UTR_end = $utr->end(); # this will return the absolute ending positions in chromosome of the UTR exons #skip if UTR and variant do not overlap next if $pos > $UTR_end || $pos_end < $UTR_start; #avoid storing UTRs if variant goes over the maximum allowed UTR overlap percentage my $UTR_length = $UTR_end - $UTR_start + 1; my $UTR_overlap_start = max($UTR_start, $pos); my $UTR_overlap_end = min($UTR_end, $pos_end); my $UTR_overlap_perc = 100 * ($UTR_overlap_end - $UTR_overlap_start) / $UTR_length; next if $UTR_overlap_perc >= $self->{max_overlap}; push(@five_utr_starts, $UTR_start); push(@five_utr_ends, $UTR_end); } return {} unless (@five_utr_starts || @five_utr_ends); my @sorted_starts = sort {$a <=> $b} @five_utr_starts; my @sorted_ends = sort {$a <=> $b} @five_utr_ends; my %UTR_info = ( "gene" => $symbol, "start" => \@sorted_starts, "end" => \@sorted_ends, "seq" => $five_prime_seq, #the exon sequences in 5'UTR "strand" => $tr_strand, "cds_seq" => $cds, ); $self->{ref_coding} = $self->get_ref_coding($ref); $self->{alt_coding} = $self->get_alt_coding($alt,$tr_strand); my %uAUG_gained = {}; my %uSTOP_lost = {}; my %uAUG_lost = {}; my %uSTOP_gained = {}; my %uFrameshift = {}; my ($mut_pos, $end_pos) = $self->get_allele_exon_pos($tr_strand, $pos, $self->{ref_coding}, \%UTR_info); if(defined($mut_pos) && defined($end_pos)) { my @sequence = split //, $five_prime_seq; my $mut_utr_seq = $self->mut_utr_sequence( \@sequence,$mut_pos, $self->{ref_coding}, $self->{alt_coding}, $tr_strand ); my @utr_sequence = split //, $UTR_info{seq}; my %existing_uORF = %{$self->existing_uORF(\@utr_sequence)}; my @mut_utr_seq = split //,$mut_utr_seq; my %existing_utr_uORF = %{$self->existing_uORF(\@mut_utr_seq)}; my @start = @{$self->get_ATG_pos(\@utr_sequence)}; %uAUG_gained = %{$self->uAUG_gained(\%variant,\%UTR_info, $mut_pos, $mut_utr_seq, \%existing_utr_uORF)}; %uSTOP_lost = %{$self->uSTOP_lost(\%variant,\%UTR_info, $mut_pos, $mut_utr_seq, \%existing_uORF, \%existing_utr_uORF)}; %uAUG_lost = %{$self->uAUG_lost(\%variant,\%UTR_info, $mut_pos, $mut_utr_seq, \%existing_uORF, \@start)}; %uSTOP_gained = %{$self->uSTOP_gained(\%variant,\%UTR_info, $mut_pos, $end_pos, $mut_utr_seq, \%existing_uORF, \%existing_utr_uORF, \@start)}; %uFrameshift = %{$self->uFrameshift(\%variant,\%UTR_info, $mut_pos, $end_pos, $mut_utr_seq, \%existing_uORF, \%existing_utr_uORF, \@start)}; } my %five_prime_flag = ( "uAUG_gained" => $uAUG_gained{'uAUG_gained_flag'}, "uSTOP_lost" => $uSTOP_lost{'uSTOP_lost_flag'}, "uAUG_lost" => $uAUG_lost{'uAUG_lost_flag'}, "uSTOP_gained" => $uSTOP_gained{'uSTOP_gained_flag'}, "uFrameshift" =>$uFrameshift{'uFrameShift_flag'}, ); my %five_prime_annotation = ( "uAUG_gained" => $uAUG_gained{"uAUG_gained_effect"}, "uSTOP_lost" => $uSTOP_lost{"uSTOP_lost_effect"}, "uAUG_lost" => $uAUG_lost{"uAUG_lost_effect"}, "uSTOP_gained" => $uSTOP_gained{'uSTOP_gained_effect'}, "uFrameshift" =>$uFrameshift{'uFrameShift_effect'}, ); my $output_five_prime_flag; my $output_five_prime_annotation; foreach my $flag (keys %five_prime_flag){ if($five_prime_flag{$flag}){ $output_five_prime_flag = $five_prime_flag{$flag}; $output_five_prime_annotation = $five_prime_annotation{$flag}; } }; my %utr_effect; if ($self->{config}->{output_format} eq "json" || $self->{config}->{rest}){ %utr_effect = ( "5UTR_consequence" => $output_five_prime_flag, "5UTR_annotation" => $output_five_prime_annotation, ); } else { my @consequences; foreach my $value (keys %{$output_five_prime_annotation}){ push @consequences, join(":", map { "$_=$output_five_prime_annotation->{$value}{$_}" } keys %{$output_five_prime_annotation->{$value}}); }; my $consequence; $consequence = join(",", @consequences); %utr_effect = ( "5UTR_consequence" => defined($output_five_prime_flag) ? $output_five_prime_flag : "-", "5UTR_annotation" => defined($consequence) && $consequence ne '' ? $consequence : "-", ); } my $existing_uORF_num = $self->count_number_ATG($five_prime_seq); my $output ={%utr_effect, %$existing_uORF_num}; return $output? $output: {}; } ################## # Main functions # ################## sub uAUG_gained { # Description: annotate if a five_prime_UTR_variant creates ATG my ($self, $variant_info,$UTR_info, $mut_pos, $mut_utr_seq, $existing_utr_uorf) = @_; my $pos = $variant_info->{pos}; my $ref = $variant_info->{ref}; my $alt = $variant_info->{alt}; my @sequence = split //, $UTR_info->{seq}; my $strand = $UTR_info->{strand}; #return annotators my $uAUG_gained_DistanceToCDS = ""; # the distance between the gained uAUG to CDS my $uAUG_gained_KozakContext = ""; # the Kozak context sequence of the gained uAUG my $uAUG_gained_KozakStrength = ""; # the Kozak strength of the gained uAUG my $uAUG_gained_type = ""; # the type of uORF created - any of the following: uORF, inframe_oORF,OutOfFrame_oORF my $uAUG_gained_DistanceFromCap = ""; # the distance between the gained uAUG to the start of the five prime UTR my $uAUG_gained_DistanceToStop = ""; #the distance between the gained uAUG to stop codon (could be in CDS) #indicate whether the variant creates a ATG my $flag = 0; my $current_kozak = ""; my $current_kozak_strength =""; #the relative position of input variant in the UTR sequence my %result = (); #result is a hash table with two elements: $flag and $output_effects my $output_flag = ""; my %output_effects; my $ref_coding = $self->{ref_coding}; my $alt_coding = $self->{ref_coding}; my @mut_utr_seq = split //,$mut_utr_seq; my $mut_utr_length = @mut_utr_seq; #get the nt sequence that might have the pos_A for a new ATG codon my @mut_seq = @mut_utr_seq[$mut_pos-2..$mut_pos+length($alt_coding)+1]; my @mut_atg_pos = @{$self->get_ATG_pos(\@mut_seq)}; # check whether there is a ATG codon is to check the length of the array #if there is a ATG codon, flag it and output the relative positive of A my $pos_A; if (@mut_atg_pos){ #get the pos of A (in ATG) in the UTR sequence $pos_A=$mut_pos-2+$mut_atg_pos[0]; $flag=1; } if ($flag){ ################################################################################ #annotator 1: get the distance to the start codon of the main ORF ################################################################################ $uAUG_gained_DistanceToCDS = $mut_utr_length-$pos_A; $uAUG_gained_DistanceFromCap = $pos_A; ################################################################################ #annotator 2: determine kozak context; ################################################################################ if ((($pos_A-3)>=0)&&($mut_utr_seq[($pos_A+3)])){ $current_kozak = $mut_utr_seq[($pos_A-3)].$mut_utr_seq[($pos_A-2)].$mut_utr_seq[$pos_A-1]."ATG".$mut_utr_seq[$pos_A+3]; } else{ $current_kozak = '-'; } #get the strength of kozak context if ($current_kozak !~ /-/){ my @split_kozak = split //, $current_kozak; $current_kozak_strength = 1; if ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))&&($split_kozak[6] eq 'G')){ $current_kozak_strength = 3; } elsif ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))||($split_kozak[6] eq 'G')){ $current_kozak_strength = 2; } } $uAUG_gained_KozakContext=$current_kozak; $uAUG_gained_KozakStrength=$self->{kozak_strength}{$current_kozak_strength}? $self->{kozak_strength}{$current_kozak_strength}:$current_kozak_strength; ################################################################################ #annotator 3: Type of new ORF with respect the main ORF: uORF/Overlapping_inFrame/Overlapping_outofframe ################################################################################ if(exists($existing_utr_uorf->{$pos_A})){ #if there is stop codon within 5'UTR $uAUG_gained_type = "uORF"; } else{ if(($mut_utr_length - $pos_A) % 3){ $uAUG_gained_type = "OutOfFrame_oORF"; } else{ $uAUG_gained_type = "inFrame_oORF"; } } my @overlapping_seq = split //, $mut_utr_seq.$UTR_info->{cds_seq}; my %existing_uORF = %{$self->existing_uORF(\@overlapping_seq)}; if(exists($existing_uORF{$pos_A})){ my @stop_pos_array = sort{$a<=>$b}@{$existing_uORF{$pos_A}}; my $stop_pos = $stop_pos_array[0]; $uAUG_gained_DistanceToStop = $stop_pos-$pos_A; } else { $uAUG_gained_DistanceToStop = "NA"; } my %uORF_effect = ( "KozakContext" => $uAUG_gained_KozakContext, "KozakStrength" => $uAUG_gained_KozakStrength, "DistanceToCDS" => $uAUG_gained_DistanceToCDS, "type" => $uAUG_gained_type, "DistanceToStop" => $uAUG_gained_DistanceToStop, "CapDistanceToStart" => $uAUG_gained_DistanceFromCap, ); $output_flag = "5_prime_UTR_premature_start_codon_gain_variant"; my $size = (keys %output_effects) + 1; $output_effects{$size} = \%uORF_effect; } $result{'uAUG_gained_flag'} = $output_flag; $result{'uAUG_gained_effect'} = \%output_effects; return \%result; } sub uSTOP_gained { # Description: annotate whether a five_prime_UTR_variant creates new stop codon. It only evaluate SNVs. my ($self, $variant_info,$UTR_info, $mut_pos, $end_pos, $mut_utr_seq, $existing_ref_uORF, $mut_uORF, $start) = @_; my $chr = $variant_info->{chr}; my $pos = $variant_info->{pos}; my $ref = $variant_info->{ref}; my $alt = $variant_info->{alt}; my @sequence = split //, $UTR_info->{seq}; my $strand = $UTR_info->{strand}; my $utr_length = @sequence; #return annotators my $uSTOP_gained_ref_StartDistanceToCDS = ""; # the distance between the uAUG of the disrupting uORF to CDS my $uSTOP_gained_KozakContext = ""; # the Kozak context sequence of the disrupting uORF my $uSTOP_gained_KozakStrength = ""; # the Kozak strength of the the disrupting uORF my $uSTOP_gained_ref_type = ""; # the type of uORF being disrupted - any of the following: uORF, inframe_oORF,OutOfFrame_oORF my $uSTOP_gained_newSTOPDistanceToCDS = ""; # the distance between the gained uSTOP to the start of the CDS my $uSTOP_gained_evidence = ""; # the translation evidence of the disrupted uORF #indicate whether the variant creates a stop codon my $flag = 0; my $current_kozak = ""; my $current_kozak_strength =""; #the relative position of input variant in the UTR sequence my %result = (); #result is a hash table with two elements: $flag and $output_effects my $output_flag = 0; my %output_effects; my $ref_coding = $self->{ref_coding}; my $alt_coding = $self->{alt_coding}; #only evaluate SNVs and MNVs, for deletion and insertion it would be evaluated in the uframeshift return{} unless(length($ref_coding) eq length($alt_coding)); # the length of the 5'UTR won't change so as the coordinates of the nts my @mut_utr_seq = split //,$mut_utr_seq; my $mut_utr_length = @mut_utr_seq; my @start = @{$start}; if(@start){ for (my $i=0;$i<@start;$i++){ $flag=0; my $start_pos = $start[$i]; # to set the range for searching new STOP codon: not including start_codon, (1) for uORFs - start_codon...stop_codon (2) for oORFs - start_codon...end of 5'UTR my $check_point = $start_pos; # the checking end point of eligible area #For uORF: start_pos .. check_point #For overlapping ORF: start_pos .. 3' end of 5'UTR sequence #if the variant is entirely in the uORF (within 5'UTR if(exists($existing_ref_uORF->{$start_pos})) { my @stops = sort {$a <=> $b} @{$existing_ref_uORF->{$start_pos}}; $check_point = $stops[0]-1; } else{ #if the existing uORF is an oORF $check_point = $utr_length-1; } # only check the ones between start and stop codons # ignore the cases at the boundary of start and stop codon since the effect on start/stop would be evaluated anyway next if(($mut_pos<$start_pos+3)|($end_pos>$check_point)); #check whether there are new stop codon induced by this mutation my @mut_stops; next unless(exists($mut_uORF->{$start_pos})); @mut_stops = sort {$a <=> $b} @{$mut_uORF->{$start_pos}}; my $mut_stop = $mut_stops[0]; if($mut_stop<$check_point){ $flag=1; } if($flag){ #getting the Kozak context and Kozak strength of the start codon if ((($start_pos-3)>=0)&&($sequence[($start_pos+3)])){ $current_kozak = $sequence[($start_pos-3)].$sequence[($start_pos-2)].$sequence[$start_pos-1]."ATG".$sequence[$start_pos+3]; } else{ $current_kozak = '-'; } if ($current_kozak !~ /-/){ my @split_kozak = split //, $current_kozak; $current_kozak_strength = 1; if ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))&&($split_kozak[6] eq 'G')){ $current_kozak_strength = 3; } elsif ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))||($split_kozak[6] eq 'G')){ $current_kozak_strength = 2; } } $uSTOP_gained_KozakContext=$current_kozak; $uSTOP_gained_KozakStrength=$self->{kozak_strength}{$current_kozak_strength}? $self->{kozak_strength}{$current_kozak_strength}:$current_kozak_strength; #the annotation of the original uORF if (exists($existing_ref_uORF->{$start_pos})){ #if there is stop codon within 5'UTR $uSTOP_gained_ref_type = "uORF" } elsif (($utr_length-$start_pos) % 3){ $uSTOP_gained_ref_type = "OutOfFrame_oORF"; } else { $uSTOP_gained_ref_type = "InFrame_oORF"; } $uSTOP_gained_ref_StartDistanceToCDS = $utr_length - $start_pos; $uSTOP_gained_newSTOPDistanceToCDS = $mut_utr_length - $mut_stop; #find evidence in added reference file $uSTOP_gained_evidence = (exists $self->{uORF_evidence})? $self->find_uorf_evidence($UTR_info,$chr,$start_pos): "NA"; my %uORF_effect = ( "ref_type" => $uSTOP_gained_ref_type, "ref_StartDistanceToCDS" => $uSTOP_gained_ref_StartDistanceToCDS, "newSTOPDistanceToCDS" => $uSTOP_gained_newSTOPDistanceToCDS, "KozakContext" => $uSTOP_gained_KozakContext, "KozakStrength" => $uSTOP_gained_KozakStrength, "Evidence" => $uSTOP_gained_evidence, ); $output_flag = "5_prime_UTR_uORF_stop_codon_gain_variant"; my $size = (keys %output_effects) + 1; $output_effects{$size} = \%uORF_effect; } } } $result{'uSTOP_gained_flag'} = $output_flag; $result{'uSTOP_gained_effect'} = \%output_effects; return \%result; } sub uSTOP_lost { # Description: annotate if a five_prime_UTR_varint removes a stop codon of an existing uORF (given that uORF does not not change) my ($self, $variant_info, $UTR_info, $mut_pos, $mut_utr_seq, $existing_uORF, $mut_uORF) = @_; my $chr = $variant_info->{chr}; my $pos = $variant_info->{pos}; my $ref = $variant_info->{ref}; my $alt = $variant_info->{alt}; my @sequence = split //, $UTR_info->{seq}; my $strand = $UTR_info->{strand}; #return annotators my $uSTOP_lost_AltStop = ""; #whether there is an alternative stop codon downstream my $uSTOP_lost_AltStopDistanceToCDS = ""; # the distance between alternative stop codon and the CDS my $uSTOP_lost_FrameWithCDS = ""; # the frame of the uORF with respect to CDS my $uSTOP_lost_KozakContext = ""; # the Kozak context sequence of the uORF with the lost stop codon my $uSTOP_lost_KozakStrength = ""; # the Kozak strength of the uORF with the lost stop codon my $uSTOP_lost_evidence = ""; # whether this uORF has any evidence of being translated ### my $current_kozak = ""; my $current_kozak_strength = ""; #indicate whether the variant ever disrupts a stop codon my $output_flag = ""; #indicate whether the variant ever disrupts the uORF that being analyzed my $flag_uORF=0; #the relative position of input variant in the UTR sequence my %result = (); #result is a hash table with two elements: $flag and $output_effects my %output_effects; my @stop_codons = ("TAA","TGA","TAG"); my $ref_coding = $self->{ref_coding}; my $alt_coding = $self->{alt_coding}; my @mut_utr_seq = split //,$mut_utr_seq; my $length = @mut_utr_seq; my @start = sort {$a <=> $b} (keys %{$existing_uORF}); if(%{$existing_uORF}){ for (my $i=0;$i<@start;$i++){ $flag_uORF=0; my $start_pos = $start[$i]; my @stops = sort {$a <=> $b} @{$existing_uORF->{$start_pos}}; my $stop_pos=$stops[0]; next if ($mut_pos-$stop_pos>2); #for snps and deletion next if (defined($ref_coding) & $mut_pos+length($ref_coding)-1<$stop_pos); next if ($mut_pos<$stop_pos & !defined($ref_coding)); #for insertion #for deletion, it definitely disrupting the stop codon. if (length($alt_coding) eq 0) { $flag_uORF=1; } else { my $mut_codon = $mut_utr_seq[$stop_pos].$mut_utr_seq[$stop_pos+1].$mut_utr_seq[$stop_pos+2]; next if (grep( /^$mut_codon$/, @stop_codons)); $flag_uORF=1; } if($flag_uORF){ #getting the Kozak context and Kozak strength of the start codon if ((($start_pos-3)>=0)&&($sequence[($start_pos+3)])){ $current_kozak = $sequence[($start_pos-3)].$sequence[($start_pos-2)].$sequence[$start_pos-1]."ATG".$sequence[$start_pos+3]; } else{ $current_kozak = '-'; } if ($current_kozak !~ /-/){ my @split_kozak = split //, $current_kozak; $current_kozak_strength = 1; if ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))&&($split_kozak[6] eq 'G')){ $current_kozak_strength = 3; } elsif ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))||($split_kozak[6] eq 'G')){ $current_kozak_strength = 2; } } $uSTOP_lost_KozakContext=$current_kozak; $uSTOP_lost_KozakStrength=$self->{kozak_strength}{$current_kozak_strength}? $self->{kozak_strength}{$current_kozak_strength}:$current_kozak_strength; # the sequence before mut_pos should not be changed. Thus start_pos shall still correspond to a start codon; # if the sequence is indeed very short as such ATGTGA my @mut_stops; if(exists($mut_uORF->{$start_pos})){ @mut_stops = sort {$a <=> $b} @{$mut_uORF->{$start_pos}} } if (@mut_stops>0){ $uSTOP_lost_AltStop = "True"; $uSTOP_lost_AltStopDistanceToCDS = $length-$mut_stops[0]; } #if there is no alternative stop codon else{ $uSTOP_lost_AltStop = "False"; $uSTOP_lost_AltStopDistanceToCDS = "NA"; } if (($length-$start_pos) % 3){ $uSTOP_lost_FrameWithCDS = "outOfFrame"; } else { $uSTOP_lost_FrameWithCDS = "inFrame"; } #find evidence from sorf $uSTOP_lost_evidence= (exists $self->{uORF_evidence})? $self->find_uorf_evidence($UTR_info,$chr,$start_pos): "NA"; my %uORF_effect = ( "AltStop" => $uSTOP_lost_AltStop, "AltStopDistanceToCDS" => $uSTOP_lost_AltStopDistanceToCDS, "FrameWithCDS" => $uSTOP_lost_FrameWithCDS, "KozakContext" => $uSTOP_lost_KozakContext, "KozakStrength" => $uSTOP_lost_KozakStrength, "Evidence" => $uSTOP_lost_evidence, ); $output_flag = "5_prime_UTR_uORF_stop_codon_loss_variant"; my $size = (keys %output_effects) + 1; $output_effects{$size} = \%uORF_effect; } } } $result{'uSTOP_lost_flag'} = $output_flag; $result{'uSTOP_lost_effect'} = \%output_effects; return \%result; } sub uAUG_lost { # Description: annotate if a five_prime_UTR_varint removes a start codon of an existing uORF my ($self, $variant_info, $UTR_info, $mut_pos, $mut_utr_seq, $existing_ref_uORF, $start) = @_; my $chr = $variant_info->{chr}; my $pos = $variant_info->{pos}; my $ref = $variant_info->{ref}; my $alt = $variant_info->{alt}; my @sorted_starts = @{$UTR_info->{start}}; my @sequence = split //, $UTR_info->{seq}; my $utr_length = @sequence; my $strand = $UTR_info->{strand}; #return annotators my $uAUG_lost_type = ""; # the uORF type with the reference allele - uORF, inframe_oORF, outOfFrame_oORF my $uAUG_lost_DistanceToCDS = ""; # the distance of the lost uAUG to CDS my $uAUG_lost_DistanceToSTOP = ""; # the distance of the lost uAUG to stop codon (could be in CDS) my $uAUG_lost_KozakContext = ""; # the Kozak context sequence of the lost uAUG my $uAUG_lost_KozakStrength = ""; # the strength of KozakContext of the lost uAUG my $uAUG_lost_evidence = ""; # whether this uAUG has any evidence of being translated my $current_kozak = ""; my $current_kozak_strength = ""; #indicate whether the variant ever disrupts a stop codon my $output_flag = ""; #indicate whether the variant ever disrupts the uORF that being analyzed my $flag_uORF=0; #the relative position of input variant in the UTR sequence my %result = (); #result is a hash table with two elements: $flag and $output_effects my %output_effects; my $ref_coding = $self->{ref_coding}; my $alt_coding = $self->{alt_coding}; my @mut_utr_seq = split //,$mut_utr_seq; my @start = @{$start}; for (my $i=0;$i<@start;$i++){ $flag_uORF=0; my $start_pos = $start[$i]; next if ($mut_pos-$start_pos>2); if (length($ref_coding) ne 0){ next if ($mut_pos+length($ref_coding)-1<$start_pos); } next if ($mut_pos<$start_pos); #for insertion #for deletion, it definitely disrupting the uORF. if (length($alt_coding) eq 0) { $flag_uORF=1; } else { my $mut_codon = $mut_utr_seq[$start_pos].$mut_utr_seq[$start_pos+1].$mut_utr_seq[$start_pos+2]; next if ($mut_codon eq "ATG"); $flag_uORF=1; } if($flag_uORF){ #getting the Kozak context and Kozak strength of the lost start codon if ((($start_pos-3)>=0)&&($sequence[($start_pos+3)])){ $current_kozak = $sequence[($start_pos-3)].$sequence[($start_pos-2)].$sequence[$start_pos-1]."ATG".$sequence[$start_pos+3]; } else { $current_kozak = '-'; } if ($current_kozak !~ /-/){ my @split_kozak = split //, $current_kozak; $current_kozak_strength = 1; if ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))&&($split_kozak[6] eq 'G')){ $current_kozak_strength = 3; } elsif ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))||($split_kozak[6] eq 'G')){ $current_kozak_strength = 2; } } $uAUG_lost_KozakContext=$current_kozak; $uAUG_lost_KozakStrength=$self->{kozak_strength}{$current_kozak_strength}? $self->{kozak_strength}{$current_kozak_strength}:$current_kozak_strength; #check what kind of uORF does that correspond to? #first check whether it's overlapping with CDS if (exists($existing_ref_uORF->{$start_pos})){ #if there is stop codon within 5'UTR $uAUG_lost_type = "uORF" } elsif(($utr_length-$start_pos) % 3){ $uAUG_lost_type = "OutOfFrame_oORF"; } else { $uAUG_lost_type = "InFrame_oORF"; } $uAUG_lost_DistanceToCDS = $utr_length - $start_pos; my @overlapping_seq = split //, $UTR_info->{seq}.$UTR_info->{cds_seq}; my %existing_overlapping_uORF = %{$self->existing_uORF(\@overlapping_seq)}; if(exists($existing_overlapping_uORF{$start_pos})){ my @stop_pos_array = sort{$a<=>$b}@{$existing_overlapping_uORF{$start_pos}}; my $stop_pos = $stop_pos_array[0]; $uAUG_lost_DistanceToSTOP = $stop_pos-$start_pos; } else { $uAUG_lost_DistanceToSTOP = "NA" } $uAUG_lost_evidence= (exists $self->{uORF_evidence})? $self->find_uorf_evidence($UTR_info,$chr,$start_pos): "NA"; my %uORF_effect = ( "type" => $uAUG_lost_type, "CapDistanceToStart" =>$start_pos, "DistanceToCDS" => $uAUG_lost_DistanceToCDS, "DistanceToStop" => $uAUG_lost_DistanceToSTOP, "KozakContext" => $uAUG_lost_KozakContext, "KozakStrength" => $uAUG_lost_KozakStrength, "Evidence" => $uAUG_lost_evidence, ); $output_flag = "5_prime_UTR_premature_start_codon_loss_variant"; my $size = (keys %output_effects) + 1; $output_effects{$size} = \%uORF_effect; } } $result{'uAUG_lost_flag'} = $output_flag; $result{'uAUG_lost_effect'} = \%output_effects; return \%result; } sub uFrameshift { # Description: annotate if a five_prime_UTR_varint create a frameshift in existing uORFs my ($self, $variant_info, $UTR_info, $mut_pos, $end_pos, $mut_utr_seq, $existing_uORF, $mut_uORF, $start) = @_; my $chr = $variant_info->{chr}; my $pos = $variant_info->{pos}; my $ref = $variant_info->{ref}; my $alt = $variant_info->{alt}; my @sequence = split //, $UTR_info->{seq}; my $strand = $UTR_info->{strand}; my $utr_length = @sequence; #return annotators my $uFrameshift_ref_type = ""; # the type of uORF with the reference allele my $uFrameshift_ref_type_length = ""; # the length of the uORF with of the reference allele. my $uFrameshift_StartDistanceToCDS = ""; # the distance between the start codon of the disrupted uORF and CDS my $uFrameshift_alt_type = ""; # the type of uORF with the alternative allele my $uFrameshift_alt_type_length = ""; # the length of the uORF with of the alternative allele my $uFrameshift_KozakContext = ""; # the Kozak context sequence of the disrupted uORF my $uFrameshift_KozakStrength = ""; # the Kozak strength of the disrupted uORF my $uFrameshift_evidence = ""; # whehter there is translation evidence for the disrupted uORF #indicate whether the variant ever introduce a frameshift variant my $flag_uORF; my $output_flag = ""; my $current_kozak=""; my $current_kozak_strength=""; my %result = (); #result is a hash table with two elements: $flag and $output_effects my %output_effects; my $ref_coding = $self->{ref_coding}; my $alt_coding = $self->{alt_coding}; #skip alleles with same length return {} unless(length($ref_coding) ne length($alt_coding)); #if it's a deletion at the boundary of exon and intron, we would skip the annotation my @mut_utr_seq = split //,$mut_utr_seq; my $length = @mut_utr_seq; my @start = @{$start}; #check for each uORF if(@start){ for (my $i=0;$i<@start;$i++){ $flag_uORF=0; my $start_pos = $start[$i]; my $check_point = $start_pos; # the checking end point of eligible area #For uORF: start_pos .. check_point #For overlapping ORF: start_pos .. 3' end of 5'UTR sequence #if the variant is entirely in the uORF (within 5'UTR) if(exists($existing_uORF->{$start_pos})) { my @stops = sort {$a <=> $b} @{$existing_uORF->{$start_pos}}; $check_point = $stops[0]-1; } else { #if the existing uORF is an oORF $check_point = $utr_length-1; } # only check the ones between start and stop codons and # ignore the cases at the boundary of start and stop codon since the effect on start/stop would be evaluated anyway if(($mut_pos>=$start_pos+3)&($end_pos<=$check_point)){ #snp and insertion could only be annotated here $flag_uORF = abs(length($ref_coding)-length($alt_coding))%3; } if($flag_uORF){ #getting the Kozak context and Kozak strength of the start codon if ((($start_pos-3)>=0)&&($sequence[($start_pos+3)])){ $current_kozak = $sequence[($start_pos-3)].$sequence[($start_pos-2)].$sequence[$start_pos-1]."ATG".$sequence[$start_pos+3]; } else{ $current_kozak = '-'; } if ($current_kozak !~ /-/){ my @split_kozak = split //, $current_kozak; $current_kozak_strength = 1; if ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))&&($split_kozak[6] eq 'G')){ $current_kozak_strength = 3; } elsif ((($split_kozak[0] eq 'A')||($split_kozak[0] eq 'G'))||($split_kozak[6] eq 'G')){ $current_kozak_strength = 2; } } $uFrameshift_KozakContext=$current_kozak; $uFrameshift_KozakStrength=$self->{kozak_strength}{$current_kozak_strength}? $self->{kozak_strength}{$current_kozak_strength}:$current_kozak_strength; #the annotation of the original uORF my @ref_overlapping_seq = split //, $UTR_info->{seq}.$UTR_info->{cds_seq}; my %ref_existing_oORF = %{$self->existing_uORF(\@ref_overlapping_seq)}; if (exists($existing_uORF->{$start_pos})){ #if there is stop codon within 5'UTR $uFrameshift_ref_type = "uORF"; } elsif (($utr_length-$start_pos) % 3) { $uFrameshift_ref_type = "OutOfFrame_oORF"; } else{ $uFrameshift_ref_type = "InFrame_oORF"; } if (exists($ref_existing_oORF{$start_pos})){ my @stops = sort {$a <=> $b} @{$ref_existing_oORF{$start_pos}}; $uFrameshift_ref_type_length = $stops[0]-$start_pos+3; } else { $uFrameshift_ref_type_length = "NA"; } $uFrameshift_StartDistanceToCDS = $utr_length - $start_pos; my @alt_overlapping_seq = split //, $mut_utr_seq.$UTR_info->{cds_seq}; my %alt_existing_oORF = %{$self->existing_uORF(\@alt_overlapping_seq)}; if (exists($mut_uORF->{$start_pos})){ $uFrameshift_alt_type = "uORF"; } elsif(($length-$start_pos)%3){ $uFrameshift_alt_type = "OutOfFrame_oORF"; } else { $uFrameshift_alt_type = "InFrame_oORF"; } #get the length if (exists($alt_existing_oORF{$start_pos})){ my @stops = sort {$a <=> $b} @{$alt_existing_oORF{$start_pos}}; $uFrameshift_alt_type_length = $stops[0]-$start_pos+3; } else { $uFrameshift_alt_type_length = "NA"; } #find evidence in added reference file $uFrameshift_evidence= (exists $self->{uORF_evidence})? $self->find_uorf_evidence($UTR_info,$chr,$start_pos): "NA"; my %uORF_effect = ( "ref_type" => $uFrameshift_ref_type, "ref_type_length" => $uFrameshift_ref_type_length, "ref_StartDistanceToCDS" => $uFrameshift_StartDistanceToCDS, "alt_type" => $uFrameshift_alt_type, "alt_type_length" => $uFrameshift_alt_type_length, "KozakContext" => $uFrameshift_KozakContext, "KozakStrength" => $uFrameshift_KozakStrength, "Evidence" => $uFrameshift_evidence, ); $output_flag = "5_prime_UTR_uORF_frameshift_variant"; my $size = (keys %output_effects) + 1; $output_effects{$size} = \%uORF_effect; } } } $result{'uFrameShift_flag'} = $output_flag; $result{'uFrameShift_effect'} = \%output_effects; return \%result; } ######### # Utils # ######### =head2 count_number_ATG Arg [1] : $seq Description: Count the number of existing ATGs in the five prime UTR sequence. Return: returns the number of uORFs and the number of oORFs Returntype : hashref =cut sub count_number_ATG { my ($self,$seq) = @_; my @sequence = split //, $seq; my $length = @sequence; my @atg_pos = @{$self->get_ATG_pos(\@sequence)}; my @mes_pos = @{$self->get_stopcodon_pos(\@sequence)}; my $inframe_stop_num=0; my $outofframe_atg_num=0; my $inframeORF_num=0; foreach my $atg (@atg_pos){ my $flag=0; #indicate whether there is a stop codon with respect to this ATG foreach my $mes (@mes_pos){ unless (($mes-$atg) % 3 || ($mes-$atg) < 0 ){ $flag=1; last; }; } #if there is no stop codon, then look at whether it's Out_of_frame or Inframe if($flag==1){ $inframe_stop_num++; } else { (($length-$atg) % 3)? $outofframe_atg_num++ : $inframeORF_num++; } } my %existing_uORF_num = ( "Existing_uORFs" => $inframe_stop_num, "Existing_OutOfFrame_oORFs" => $outofframe_atg_num, "Existing_InFrame_oORFs" => $inframeORF_num, ); return \%existing_uORF_num; } =head2 existing_uORF Arg [1] : $seq Description: obtaining the relative coordinates of start and end pos of existing uORF in the five prime UTR sequence. Return: returns as the key the position of the first nucleotide of start codon and the value is all the positions of the first nucleotide of the stop codon. Returntype : hashref =cut sub existing_uORF { my ($self,$seq) = @_; my @atg_pos = @{$self->get_ATG_pos($seq)}; my @mes_pos = @{$self->get_stopcodon_pos($seq)}; my %uORF; foreach my $atg (@atg_pos){ my @stop_pos; foreach my $mes (@mes_pos){ push @stop_pos, $mes unless ((($mes-$atg) % 3) || ($mes-$atg)<0); } $uORF{$atg} = [@stop_pos] if (@stop_pos); } return \%uORF; } =head2 get_ATG_pos Arg [1] : $seq Description: get all the relative position of A in ATG of the five prime UTR sequence Return: returns ATG position in sequence. Returntype : listref =cut sub get_ATG_pos { my ($self,$seq) = @_; my @sequence = @{$seq}; my $length = @sequence; my $seq_str=join '', @sequence; my @atg_pos = grep {(substr ($seq_str, $_,3) eq 'ATG')} 0..($length); return \@atg_pos; } =head2 get_stopcodon_pos Arg [1] : $seq Description: get all the relative position of stop codons in the five prime UTR sequence Return: returns stop-codon position in sequence. Returntype : listref =cut sub get_stopcodon_pos { my ($self,$seq) = @_; my @sequence = @{$seq}; my $length = @sequence; my @met_pos; if($length){ for (my $seq_n=0; $seq_n<$length; $seq_n++){ if ((($sequence[$seq_n] eq 'T')&&($sequence[$seq_n+1] eq 'A')&&($sequence[$seq_n+2] eq 'A')) ||(($sequence[$seq_n] eq 'T')&&($sequence[$seq_n+1] eq 'A')&&($sequence[$seq_n+2] eq 'G')) ||(($sequence[$seq_n] eq 'T')&&($sequence[$seq_n+1] eq 'G')&&($sequence[$seq_n+2] eq 'A'))){ push @met_pos,$seq_n; } } } return \@met_pos; } sub get_ref_coding { my ($self,$ref) = @_; my $ref_coding = $ref; if($ref_coding eq "-"){ #for insertion $ref_coding = ""; } return $ref_coding; } sub get_alt_coding { my ($self,$alt,$strand) = @_; my $alt_coding = $alt; $alt_coding = "" if($alt_coding eq "-"); reverse_comp(\$alt_coding) if($strand < 0); return $alt_coding; } =head2 mut_utr_sequence Description: get the mutated five prime UTR sequence Return: Mutated UTR sequence Returntype : String =cut sub mut_utr_sequence { my ($self,$seq,$mut_pos,$ref_coding,$alt_coding,$strand) = @_; my @sequence = @{$seq}; my $length = @sequence; #order from 5' to 3' my $mut_sequence; $mut_sequence = (join "", @sequence[0..$mut_pos-1]).$alt_coding.(join "", @sequence[$mut_pos+length($ref_coding)..$length-1]); return $mut_sequence; } sub transform_hash_to_string { my ($self,$hash) = @_; my %hash = %{$hash}; my $output_str = ""; foreach my $key (sort keys %hash){ $output_str = $output_str.$key.":".$hash{$key}.","; }; chop($output_str) if($output_str); return $output_str; } sub utr_exon_position { my ($self,$UTR_info) = @_; my @sorted_starts = @{$UTR_info->{start}}; my @sorted_ends = @{$UTR_info->{end}}; my $num_exons = @sorted_starts; my $strand = $UTR_info->{strand}; my %chr_position; my $utr_position = 0; if ($strand == 1){ #create a map from chromosome position to UTR position for (my $m=0; $m<$num_exons; $m++){ for (my $p=$sorted_starts[$m]; $p<=$sorted_ends[$m]; $p++){ $chr_position{$p}=$utr_position; $utr_position++; } } } if ($strand == -1){ #create a map from chromosome position to UTR position #the exons were arranged in increasing order above for (my $m=$num_exons-1; $m>=0; $m--){ for (my $p=$sorted_ends[$m]; $p>=$sorted_starts[$m]; $p--){ $chr_position{$p}=$utr_position; $utr_position++; } } } return \%chr_position; } sub chr_position { my ($self,$UTR_info) = @_; my @sorted_starts = @{$UTR_info->{start}}; my @sorted_ends = @{$UTR_info->{end}}; my $num_exons = @sorted_starts; my $strand = $UTR_info->{strand}; my %utr_position; my $utr_position = 0; if ($strand == 1){ #create a map from chromosome position to UTR position for (my $m=0; $m<$num_exons; $m++){ for (my $p=$sorted_starts[$m]; $p<=$sorted_ends[$m]; $p++){ $utr_position{$utr_position}=$p; $utr_position++; } } } if ($strand == -1){ #create a map from chromosome position to UTR position #the exons were arranged in increasing order above for (my $m=$num_exons-1; $m>=0; $m--){ for (my $p=$sorted_ends[$m]; $p>=$sorted_starts[$m]; $p--){ $utr_position{$utr_position}=$p; $utr_position++; } } } return \%utr_position; } sub get_allele_exon_pos { # return the 5' start and end position at the exon of the ref allele my ($self, $strand, $pos, $ref_coding, $UTR_info) = @_; my %utr_exon_pos = %{$self->utr_exon_position($UTR_info)}; my $ref_start; my $ref_end; if ($strand == 1) { $ref_start = $utr_exon_pos{$pos}; if(length($ref_coding)){ $ref_end = $utr_exon_pos{$pos + length($ref_coding) - 1}; } else { $ref_end = $ref_start; } } else { $ref_start = $utr_exon_pos{$pos+length($ref_coding)-1}; if(length($ref_coding)){ $ref_end = $utr_exon_pos{$pos}; } else { $ref_end = $ref_start; } } return ($ref_start, $ref_end); } sub find_uorf_evidence { my ($self,$UTR_info,$chr,$start_pos)=@_; my %utr_pos = %{$self->chr_position($UTR_info)}; my $start_chr_pos = $utr_pos{$start_pos}; my $query = ($chr=~/chr/i)? $chr . ":" . $start_chr_pos : "chr" . $chr . ":" . $start_chr_pos; my $evidence = $self->{uORF_evidence}->{$query}? "True" : "False"; return $evidence; } 1;