wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip
unzip snpEff_latest_core.zip

java -jar /media/sf_dna/apps/snpEff/snpEff.jar download -v hg38

java -Xmx8g -jar /media/sf_dna/apps/snpEff/snpEff.jar hg38 /media/sf_dna/nebula/vcf/nebula.vcf > myVariants.ann.vcf

sudo apt install default-jre
wget https://github.com/broadinstitute/gatk/releases/download/4.2.5.0/gatk-4.2.5.0.zip
unzip gatk-4.2.5.0.zip
# set alias in bash profile 
# alias gatk='/media/sf_dna/apps/gatk-4.2.5.0/gatk'

gatk VariantsToTable -V myVariants.ann.vcf -F CHROM -F POS -F TYPE -F ID -F ANN -F LOF -F NMD -GF AD -GF DP -GF GQ -GF GT -O myVariants.ann.txt

nebulaID="##########" 
sed -i "1s/$nebulaID/myVariants/g" myVariants.ann.txt

annotated_vcf_table <- "C:/Users/jmcgirr/dna/nebula/vcf/myVariants.ann.txt"

ann <- vroom(annotated_vcf_table, guess_max = 1000000)

# table based on SnpEFF documentation
ann_types <- data.frame(impact = c("HIGH","HIGH","HIGH","HIGH","HIGH","HIGH","HIGH","HIGH","HIGH","HIGH","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","MODERATE","LOW","LOW","LOW","LOW","LOW","LOW","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER","MODIFIER"),
                     ontology_term = c("chromosome_number_variation","exon_loss_variant","frameshift_variant","rare_amino_acid_variant","splice_acceptor_variant","splice_donor_variant","start_lost","stop_gained","stop_lost","transcript_ablation","3_prime_UTR_truncation&exon_loss","5_prime_UTR_truncation&exon_loss_variant","coding_sequence_variant-moderate","conservative_inframe_deletion","conservative_inframe_insertion","disruptive_inframe_deletion","disruptive_inframe_insertion","missense_variant","regulatory_region_ablation","splice_region_variant-moderate","TFBS_ablation","5_prime_UTR_premature_start_codon_gain_variant","initiator_codon_variant","splice_region_variant-low","start_retained","stop_retained_variant","synonymous_variant","3_prime_UTR_variant","5_prime_UTR_variant","coding_sequence_variant-modifier","conserved_intergenic_variant","conserved_intron_variant","downstream_gene_variant","exon_variant","feature_elongation","feature_truncation","gene_variant","intergenic_region","intragenic_variant","intron_variant","mature_miRNA_variant","miRNA","NMD_transcript_variant","non_coding_transcript_exon_variant","non_coding_transcript_variant","regulatory_region_amplification","regulatory_region_variant","TF_binding_site_variant","TFBS_amplification","transcript_amplification","transcript_variant","upstream_gene_variant"))

ann_types |> gt() |> tab_header(title = "Putative impact for sequence ontology terms output by SnpEFF") 

# count the number of variants annotated for each ontology
n_variants <- c()
for(ot in ann_types$ontology_term){
n_variants <- c(n_variants,(filter(ann, grepl(ot, ANN)) |> nrow()))
#print(ot)
}

ann_types$n_variants <- n_variants
ann_types <- arrange(ann_types, n_variants)

ann_types$ontology_term <- factor(ann_types$ontology_term, levels = ann_types$ontology_term)
ann_types$impact <- factor(ann_types$impact, levels = c("HIGH","MODERATE","LOW","MODIFIER"))
ann_types$log_n_variants <- log10(ann_types$n_variants)

p1 <- ann_types|> 
   ggplot(aes(ontology_term, log_n_variants, fill = impact)) +
   geom_col() +
   coord_flip() +
   theme_minimal() +
   theme(axis.title.x=element_text(size=14),
    axis.title.y=element_text(size=12),
    axis.title=element_text(size=14),
    axis.text=element_text(size=12),
    plot.title=element_text(size=18))+
   ylab("\nlog(number of annotated variants)")+ xlab("")+ ylim(0,8)+
   ggtitle("Number of variants showing each snpEFF annotation (log scale)")+
   geom_text(aes(label=n_variants), position=position_stack(), hjust=-0.5)+
   scale_fill_manual(values = c(red,yel,blu,gre))
print(p1)

p1 <- ann_types|> filter(impact %in% c("HIGH"), ontology_term != "missense_variant") |>
   ggplot(aes(ontology_term, n_variants, fill = impact)) +
   geom_col() +
   coord_flip() +
   theme_minimal() +
   theme(axis.title.x=element_text(size=14),
    axis.title.y=element_text(size=14),
    axis.title=element_text(size=18),
    axis.text=element_text(size=12),
    plot.title=element_text(size=18))+
   ylab("\nnumber of annotated variants")+ xlab("")+ylim(0,400)+
   ggtitle("Number of variants showing high impact snpEFF annotations")+
   geom_text(aes(label=n_variants), position=position_stack(), hjust=-0.5)+
   scale_fill_manual(values = c(red,yel,blu,gre))
print(p1)

# get table with rsids for an ontology term
annotations.for.ontology.term <- function(ontology_term){
 
  ot_table <- filter(ann, grepl(ontology_term, ANN)) 
  anns <- ot_table$ANN
  
  ot_anns <- c()
  ot_alleles <- c()
  for(ann_string in anns){
    ot_anns <- c(ot_anns,paste0(paste(strsplit(grep(ontology_term, strsplit(ann_string, ",")[[1]], value = TRUE)[1], "\\|")[[1]][c(1,2,3,4,7)], collapse = "|"),"|"))
  }
  
  ot_table$ANN <- ot_anns
  ot_table <- separate(ot_table,ANN, c("allele","annotation","impact","gene_name","feature_id"), sep = "\\|", extra = "drop")
  return(ot_table)
}

# search for SNPedia annotations for rsids output by annotations.for.ontology.term()
get.snpedia.urls <- function(ot_table){
  
  snpedia_urls <- c()
  myGTs  <- myADs <- myGQs <- c() 
  ot_table <- filter(ot_table, grepl("rs",ID))
  
  for(rsid in ot_table$ID){
  snpedia_table <- annotations(snp = rsid, output = 'snpedia')
    if(nrow(snpedia_table)>0){
      snpedia_urls <- c(snpedia_urls,gsub("\\s*\\([^\\)]+\\)","",snpedia_table$url[1]))
      myGTs <- c(myGTs, filter(ot_table, ID == rsid) |> pull(myVariants.GT))
      myADs <- c(myADs, filter(ot_table, ID == rsid) |> pull(myVariants.AD))
      myGQs <- c(myGQs, filter(ot_table, ID == rsid) |> pull(myVariants.GQ))

    }
  }
  
  return(data.frame(SNPedia_urls = snpedia_urls, Genotype_PlusOrientation = myGTs,
                    Allele_Depth = myADs, Genotype_Quality = myGQs))
  
}

make.hyperlink <-  function(myurl,mytext=myurl) {
  paste('<a href="',myurl,'">',mytext,'</a>')
}

for(i_ontology_term in (filter(ann_types, impact == "HIGH", n_variants > 0) |> pull(ontology_term))){
print(get.snpedia.urls(annotations.for.ontology.term(i_ontology_term)) |> 
        gt() |> 
        tab_header(title = paste0("SNPedia entries for variants annotated as ", i_ontology_term)) |> 
        fmt (columns = 'SNPedia_urls',fns = make.hyperlink))
}

annotated_vcf_table <- "C:/Users/jmcgirr/dna/nebula/vcf/myVariants.clinvar.txt"

ann <- vroom(annotated_vcf_table, guess_max = 1000000) |>
      filter(!is.na(ALLELEID))

patho <- filter(ann, grepl("Pathogenic",CLNSIG) | grepl("Likely_pathogenic",CLNSIG)) |> separate("ID", c("ID","unknown"), sep = ";")
#unique(patho$CLNSIG)
#nrow(patho)
as.data.frame(patho) |> gt()

print(get.snpedia.urls(patho) |> 
      gt() |> 
      tab_header(title = "SNPedia entries for variants annotated as pathoginic or likely pathogenic") |> 
      fmt (columns = 'SNPedia_urls',fns = make.hyperlink))