R/functions.R

Defines functions fast_map circos_plot depth_approx cpm_approx tpm_approx

Documented in circos_plot fast_map

tpm_approx <- function(data,length){

  RPK = data / (length/1000)
  TPM = RPK/(colSums(RPK)/1e6)
  return(TPM)
}

cpm_approx <- function(data){

  CPM = (data/colSums(data))* 1e6
  return(CPM)
}

depth_approx <- function(data){

  CPM = data / colSums(data)
  return(CPM)
}


#' Draw final circle plot
#'

circos_plot <- function(expression,cnv,variant, arriba, bed,
                        SNP_max_rad = 1.145, SNP_min_rad =1,
                        exp_max_rad = 0.97, exp_min_rad = 0.80,
                        CNV_max_rad = 0.75, CNV_min_rad = 0.4,
                        fusion_max_rad = 0.35,
                        CNV_cutoff = 1.5
                        ){


  ########## ========== Global parameters



  ### Set genes coords
  gene_order <- gsub("chr","",unique(bed$chr))

  genes_coords <- bed %>%
    #mutate(chr = factor(chr, levels = unique(bed$chr))) %>%
    distinct(gene_name, .keep_all = T) %>%
    group_by(chr) %>%
    mutate(order = seq(1:n()),
           chr = factor(gsub("chr","",chr), levels = gene_order))

  ### Set custom genome coords
  custom_genome <- genes_coords %>%
    group_by(chr) %>%
    summarise(length = length(unique(gene_name))) %>%
    column_to_rownames("chr") %>% t() %>% as.data.frame(check.names = F) %>% as.list()

  ########## ========== Prepare data

  ##### ===== Expression

  salmon_res <- expression %>%
    #select(gene_name, expression = all_of(input$sample)) %>% # DEV TO REMOVE
    inner_join(genes_coords, by = "gene_name") %>%
    arrange(chr)

  ##### ===== CNV

  cnv_res <- cnv %>% left_join(genes_coords, by = c("gene"="gene_name"),
                                 suffix = c("_cnv","_coords"))


  cnv_res_notsig <- cnv_res %>% filter(abs(CNV) < CNV_cutoff)
  cnv_res_sig <- cnv_res %>% filter(abs(CNV) >= CNV_cutoff)


  ##### ===== Mutations

  RNAmut_res <- variant %>% filter(Info == "Oncogenic", type == "Variant") %>%
    #filter(sample_id == input$sample) %>% # DEV TO REMOVE
    arrange(desc(VAF2)) %>%
    distinct(gene_mut, .keep_all = T) %>%
    left_join(genes_coords, by = c("Gene"="gene_name"))

  ##### ===== Fusions

  fusion_res <- arriba %>% select(-read_identifiers, -fusion_transcript) %>%
    #filter(sample_id == input$sample) %>% # DEV TO REMOVE
    filter(confidence == "high") %>%
    left_join(genes_coords, by = c("gene1"="gene_name")) %>%
    left_join(genes_coords, by = c("gene2"="gene_name"), suffix = c("_g1","_g2")) %>%
    na.omit()
    # separate(breakpoint1, into = c("chr1","linkpos1")) %>%
    # separate(breakpoint2, into = c("chr2","linkpos2"))



  ########## ========== Draw tracks


  ##### ===== SNP tracks

  if(nrow(RNAmut_res) == 0){

    tracks = BioCircosTextTrack("text1",
                                         'No variants found',
                                         weight = "normal",
                                         x = -0.2, y = -1.07)
    tracks = tracks + BioCircosTextTrack("text2",
                                         'by RNAMut',
                                         weight = "normal",
                                         x = -0.1, y = -1.02)

  } else {

    tracks =  BioCircosArcTrack("RNAmut_variants",
                                        chromosome = as.character(RNAmut_res$chr),
                                        starts = RNAmut_res$order,
                                        ends = RNAmut_res$order+4,
                                        labels = RNAmut_res$gene_mut,

                                        maxRadius = SNP_max_rad, minRadius = SNP_min_rad,
                                        color = "black")

  }


  ##### ===== Expression track
  tracks = tracks + BioCircosHeatmapTrack("heatmap1",
                                          as.character(salmon_res$chr),
                                          salmon_res$order-1,
                                          salmon_res$order,
                                          round(salmon_res$expression,2),
                                          labels = salmon_res$gene_name,
                                          minRadius = exp_min_rad,
                                          maxRadius = exp_max_rad)

  ##### ===== CNV tracks

  tracks = tracks + BioCircosCNVTrack("cnv_backm2",
                                      as.character(cnv_res$chromosome),
                                      cnv_res$order_coords-1,
                                      cnv_res$order_coords,
                                      rep(-2,length(as.character(cnv_res$chromosome))),
                                      color = "black",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_backm1",
                                      as.character(cnv_res$chromosome),
                                      cnv_res$order_coords-1,
                                      cnv_res$order_coords,
                                      rep(-1,length(as.character(cnv_res$chromosome))),
                                      color = "lightgrey",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_back0",
                                      as.character(cnv_res$chromosome),
                                      cnv_res$order_coords-1,
                                      cnv_res$order_coords,
                                      rep(0,length(as.character(cnv_res$chromosome))),
                                      color = "grey",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_back1",
                                      as.character(cnv_res$chromosome),
                                      cnv_res$order_coords-1,
                                      cnv_res$order_coords,
                                      rep(1,length(as.character(cnv_res$chromosome))),
                                      color = "lightgrey",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_back2",
                                      as.character(cnv_res$chromosome),
                                      cnv_res$order_coords-1,
                                      cnv_res$order_coords,
                                      rep(2,length(as.character(cnv_res$chromosome))),
                                      color = "black",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_notsig",width = 2,
                             as.character(cnv_res_notsig$chromosome),
                             cnv_res_notsig$order_coords-1,
                             cnv_res_notsig$order_coords,
                             cnv_res_notsig$CNV,
                             color = "cornflowerblue",
                             range = c(-2,2),
                             minRadius = CNV_min_rad,
                             maxRadius = CNV_max_rad)

  tracks = tracks + BioCircosCNVTrack("cnv_sig", width = 5,
                                      as.character(cnv_res_sig$chromosome),
                                      cnv_res_sig$order_coords-1,
                                      cnv_res_sig$order_coords,
                                      cnv_res_sig$CNV,
                                      color = "firebrick",
                                      range = c(-2,2),
                                      minRadius = CNV_min_rad,
                                      maxRadius = CNV_max_rad)





  ##### ===== Fusions track

  if(nrow(fusion_res) == 0){

    tracks = tracks + BioCircosTextTrack("text3",
                                        'No fusions found',
                                        weight = "normal",
                                        x = -0.2, y = 0)
    tracks = tracks + BioCircosTextTrack("text4",
                                         'by Arriba',
                                         weight = "normal",
                                         x = -0.10, y = 0.1)

  } else {

    tracks = tracks + BioCircosLinkTrack("fusions",
                                         gene1Chromosomes = as.character(fusion_res$chr_g1),
                                         gene1Starts = as.numeric(fusion_res$order_g1),
                                         gene1Ends = as.numeric(fusion_res$order_g1)+4,
                                         gene2Chromosomes = as.character(fusion_res$chr_g2),
                                         gene2Starts = as.numeric(fusion_res$order_g2),
                                         gene2Ends = as.numeric(fusion_res$order_g2)+4,
                                         labels = fusion_res$fusion,
                                         displayAxis = F, axisPadding = 6,
                                         color = "firebrick", width = "0.5em",
                                         displayLabel = F,
                                         maxRadius = fusion_max_rad)
  }

  ##### ===== Backgrounds
#
  # tracks = tracks + BioCircosBackgroundTrack("BG_SNP",
  #                                            borderColors = "black",
  #                                            fillColor = "white", borderSize = 0.6,
  #                                            minRadius = CNV_max_rad/2, maxRadius = CNV_max_rad)
  # tracks = tracks + BioCircosBackgroundTrack("BG_SNP_2",
  #                                            borderColors = "black",
  #                                            fillColor = "white", borderSize = 0.6,
  #                                            minRadius = CNV_min_rad, maxRadius = CNV_max_rad/2)
  #


  ##### ===== Final viz
  p <- BioCircos(tracks,
                 genome = custom_genome,
                 genomeFillColor = "Spectral",
                 yChr = F,
                 chrPad = 0,
                 displayGenomeBorder = F,
            genomeTicksLen = 3,
            genomeTicksTextSize = 0,
            genomeTicksScale = 50000000,
            genomeLabelTextSize = 18,
            genomeLabelDy = 0)

  return(p)



}




#' Input = df with samples as rownames and features as cols
#' Annotation = vector of annotation (in the same order as input's rownames!!)
#' anno_pal = viridis_pal, hue_pal, grey_pal...
#'


fast_map <- function(df, center_scale = T,
                     heat_colors=colorRampPalette(c('#3361A5','#248AF3','#14B3FF','#88CEEF','#C1D5DC',
                                                    '#EAD397','#FDB31A','#E42A2A','#A31D1D'))(100),
                     anno_pal = viridis_pal(),
                     Zlim = 3, method = "ward.D2", metric = "euclidean", title = "title",
                     Data_input = "Z-score expression",
                     anno_1 = NULL, anno_1_name = NULL, anno_1_cols = NULL,
                     cluster_rows = T, cluster_cols = T,
                     reorder_rows = T, reorder_cols = T,
                     row_names = T,col_names=T,
                     row_fontsize = 12, col_fontsize = 12){



#### Parameters

  message("Parameters")
  heat_col <- circlize::colorRamp2(seq(-Zlim, +Zlim, by =(2*Zlim)/99),heat_colors)


  sup_leg = list(legend_1 = Legend(col_fun = heat_col,
                                   title = Data_input,
                                   labels_gp = gpar(fontsize = 12),
                                   direction = c("horizontal"),
                                   title_gp = gpar(fontsize = 12, fontface="bold"))
  )


  message("Anno if ")
  if(is.null(anno_1)){

    sup_ha = NULL

  } else {

    if(is.numeric(anno_1) | is.double(anno_1)){

      ##### ===== 1 unique continuous annotation

      col_list <- list(supha_1 = circlize::colorRamp2(c(min(anno_1, na.rm=T),
                                                        median(anno_1,na.rm=T),
                                                        max(anno_1,na.rm=T)),
                                                      c(anno_pal(100)[1],
                                                        anno_pal(100)[50],
                                                        anno_pal(100)[100])))

      sup_ha =  HeatmapAnnotation(supha_1 = anno_1,
                                  show_annotation_name = F,
                                  col = col_list,
                                  show_legend = F,
                                  na_col = "grey")

      sup_leg$legend_2 = Legend(col_fun = col_list$supha_1,
                                title=anno_1_name,
                                direction = c("horizontal"), ncol=2,
                                labels_gp = gpar(fontsize = 12),
                                title_gp = gpar(fontsize = 12, fontface="bold"))


    } else {

      ##### ===== 1 unique categorical annotation

      ha_1 <- anno_pal(length(unique(na.omit(anno_1))))
      names(ha_1) <- unique(na.omit(anno_1))

      col_list <- list(supha_1 = ha_1)

      sup_ha =  HeatmapAnnotation(supha_1 = anno_1,
                                  show_annotation_name = F,
                                  col = col_list,
                                  show_legend = F)

      sup_leg$legend_2 = Legend(labels=names(col_list[[1]]),
                                title=anno_1_name,
                                direction = c("horizontal"), ncol=2,
                                labels_gp = gpar(fontsize = 12),
                                title_gp = gpar(fontsize = 12, fontface="bold"),
                                legend_gp = gpar(fill = as.vector(col_list[[1]])))

    }



  }

  ##### ===== Prepare df

  message("Prepare df ")

  matrix <- df
  if(center_scale == T){
    matrix <- scale(matrix)
  }


  message("set heatmap ")
  ht <-   ComplexHeatmap::Heatmap(t(matrix),
                  name = "heat",
                  col = heat_col,
                  show_column_names = col_names,
                  show_row_names = row_names,
                  #column_split = k,
                  cluster_columns = cluster_cols,
                  show_column_dend = cluster_cols,
                  column_dend_reorder =   reorder_cols,
                  clustering_distance_columns = metric,
                  clustering_method_column = method,

                  show_heatmap_legend = F,

                  show_row_dend = cluster_rows,
                  cluster_rows = cluster_rows,
                  row_dend_reorder = reorder_rows,
                  clustering_distance_rows = metric,
                  clustering_method_rows = method,
                  use_raster = TRUE,
                  raster_device = c("png"),
                  raster_quality = 2,
                  top_annotation = sup_ha,
                  column_title = title,
                  column_title_gp = gpar(fontsize = 14, fontface = "bold"),
                  row_title_gp = gpar(fontsize = 12),

                  row_names_gp = gpar(fontsize = row_fontsize),
                  column_names_gp = gpar(fontsize = col_fontsize)
                  )

### draw
  message("Plot heatmap")
  draw(ht, annotation_legend_list = sup_leg, annotation_legend_side = "bottom", newpage = T)



}
VincentAlcazer/hemRNA documentation built on Aug. 26, 2022, 4:50 a.m.