In yasmina-mekki/gwhap: Genome wide haplotype analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(gwhap)

RESOURCE_ROOT = "<YOUR_LOCAL_RESOURCE_FS>"

Installation of `gwhap`

There are two ways to obtain the gwhap source package:

Using the CRAN

install.packages('gwhap')

Using github

devtools::install_github('yasmina-mekki/gwhap')

Blocs distribution plot

Assuming that df_blocs is your bloc data frame already loaded

Blocs distribution for all chromosomes and optionnaly for different values of delta

haplotype_bloc_distribution_per_delta = haplotype_bloc_distribution(df_blocs)
save_plot(haplotype_bloc_distribution_per_delta, paste0('output_path/', 'haplotype_bloc_distribution_per_delta.png'))

Blocs distribution per chromosome and optionnaly for different values of delta

haplotype_bloc_distribution_per_chr = haplotype_bloc_distribution(df_blocs, per_chromosome=TRUE)
save_plot(haplotype_bloc_distribution_per_chr, paste0('output_path/', 'haplotype_bloc_distribution_per_chr.png'))

Genome coverage of haplotype blocs

For practical reasons, we had chosen to not include the karyplote as a function of the package. However, the code that implement this plot is available below.

#' Genome coverage of haplotype blocs
#'
#' @description Genome coverage of haplotype blocs
#'
#' @param df_blocs data frame structure containing all the blocs created
#' @param colors the colors desired represented by an integer value
#' @param verbose silent warning messages. FALSE by default.
#'
#' @return karyotype plot object
#' 
#' @import karyoploteR
#' @importFrom grDevices rainbow
#' @importFrom dplyr group_by summarise
#' @export
karyotype_plot <- function(df_blocs, colors=200, verbose=FALSE){

  # silent warning messages
  if(verbose == TRUE){options(warn=0)} else{options(warn=-1)}

  # params blocs setting
  df_blocs$lbp = df_blocs$to_bp - df_blocs$from_bp
  df_blocs$mb  = round(df_blocs$from_bp/1e6)
  #sum_l_bp_chr = df_blocs %>% group_by(chr=df_blocs$chr, mb=df_blocs$mb) %>% summarise(coverage=sum(df_blocs$lbp)/1e6)

  sum_l_bp_chr = df_blocs %>% group_by(chr=chr, mb=mb) %>% summarise(coverage=sum(lbp)/1e6)

  # params plot setting
  pp <- getDefaultPlotParams(plot.type = 1)
  pp$data1height = 200

  # init karyotype object with params plot chosen above
  karyotype_plot_obj <- plotKaryotype(chromosomes=c('autosomal'), plot.params = pp)
  kpHeatmap(karyoplot = karyotype_plot_obj,
            chr       = sum_l_bp_chr$chr,
            x0        = sum_l_bp_chr$mb*1e6,
            x1        = (sum_l_bp_chr$mb+1)*1e6,
            y         = round(sum_l_bp_chr$coverage*100),
            colors    = rainbow(colors)[round(sum_l_bp_chr$coverage*100)]) #rainbow(200, start=0.2, end=0.7)
  # add base numbers
  kpAddBaseNumbers(karyotype_plot_obj)

  return(karyotype_plot_obj)
}

Now that you loaded the function, you can produce the karyoplot as following:

karyotype_plot_obj = karyotype_plot(df_blocs)
save_plot(karyotype_plot_obj, paste0('output_path/', 'karyotype_plot.png'))

Phenotypes distribution plot

phenotype_distribution_plot = phenotype_distribution(phenotypes_path = c('phenotypes_path/phenotype.csv'))
save_plot(phenotype_distribution_plot, paste0(output, 'phenotype_distribution_plot.png'))