XCI_Escape_Removal (20231018).Rmd

For the first step for removing the XCI escape genes, we need to annotate our variants. To make things easier, Jonathan’s code (as seen in the “Plotting VCF” file) filters out the SNPs that match our criteria. These SNPs are biallelic - which serves our purpose. Now we start with loading required libraries.

# Load libraries
library(magrittr)
library(readr)
library(ggplot2)
library(tibble)
library(tidyr)
library(dplyr)
library(readr)
library(rtracklayer)
library(stringr)
library(plyranges)
library(GenomicFeatures)
library(plyranges)
library(cowplot)
library(Seurat)
library(SeuratDisk)
library(pheatmap)
library(reshape2)
library(RColorBrewer)
library(rhdf5)
library(matrixStats)
library(EGAD)

Setting options('download.file.method.GEOquery'='auto')
Setting options('GEOquery.inmemory.gpl'=FALSE)

library(GO.db)

library(org.Hs.eg.db)

# installation
# BiocManager::install("plyranges")

Let’s regenerate the combined dataframe that includes the XCI skew information.

# Let's do this analysis for both Pilot and Second data

Warning message:
In fun(libname, pkgname) : couldn't connect to display ":0"

merged_df_1 <- read.csv("/home/youkim/Bulk_RNA_seq/XCI_analysis/Pilot_Data/Merged_analysis/20230920_XCI_genes(Pilot_Data).csv")

merged_df_2 <- read.csv("/home/youkim/Bulk_RNA_seq/XCI_analysis/Second_Data/Merged_analysis_After_QC/20230921_XCI_genes(Second_Data).csv")

merged_df <- read.csv("/home/youkim/Bulk_RNA_seq/XCI_analysis/Third_Data/XCI_genes(Third_Data).csv")
merged_df <- merged_df %>% 
  dplyr::select(-X.1)

combined <- rbind(merged_df_1, merged_df_2) %>% dplyr::select(-X.1)

# merge it with the third data
combined <- rbind(combined, merged_df)

Now let’s see if we could generate a heatmap that could encompass all 22 samples based on their genes

Since there were too many NAs in the first heatmap, let’s generate a new one that has a specific NA count filter.

heatmap_data <- dcast(combined, gene_name ~ sample_id, value.var = "mean_XCI_ratio")
rownames(heatmap_data) <- heatmap_data$gene_name  # Set row names to gene names
heatmap_data <- heatmap_data[, -1]  # Remove the gene_name column

# Let's remove B5_12764 and see what happens
heatmap_data <- heatmap_data %>% 
  dplyr::select(-B5_12764)

# Count the number of NAs in each row
na_count <- rowSums(is.na(heatmap_data))

# Create a list to store the heatmaps
heatmap_list <- list()

# Create a function to generate a heatmap for a given threshold
generate_heatmap <- function(threshold) {
    
    # Filter rows with a number of NAs less than or equal to the threshold
    filtered_heatmap_data <- heatmap_data[na_count <= threshold, ]
    
    drows_3 <- dist(filtered_heatmap_data, method = "euclidean")
    # Check that we have no NAs in our distance matrix now!
    which(is.na(drows_3), arr.ind = T)
    
    colors <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100)
    
    annotation_helper <- combined %>% 
      filter(gene_name %in% rownames(filtered_heatmap_data)) %>% 
      arrange(gene_name) %>% 
      dplyr::select(gene_name, mean_exprs_ratio) 
    
    rownames(annotation_helper) <- NULL
    interim <- annotation_helper$gene_name %>% duplicated()
    annotation_helper <- annotation_helper[!interim,]
    
    # double check the order
    summary(rownames(filtered_heatmap_data) == annotation_helper$gene_name)
    
    # Color for annotations
    expression_ratio_palette <- colorRampPalette(c("blue", "white", "red"))(100)
    
    # Annotation dataframe
    annotation_df <- data.frame(Expression_Ratio = log10(annotation_helper$mean_exprs_ratio + 1))
    rownames(annotation_df) <- rownames(filtered_heatmap_data)
    
    p <- pheatmap(
      filtered_heatmap_data,
      clustering_distance_rows = drows_3,  # You can choose a different clustering method if needed
      clustering_distance_cols = "euclidean",  # You can choose a different clustering method if needed
      annotation_row = annotation_df,  # Add the expression ratio annotation
      annotation_colors = list(Expression_Ratio = expression_ratio_palette),
      col = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100),  # Color palette
      na_col = "grey",  # Set the color for NA values
      main = "Gene Expression Heatmap",
      clustering_method = "complete",
      fontsize_row = 6,  # Adjust row label size
      fontsize_col = 9,  # Adjust column label size
      cutree_rows = 2
    )
  return(p)
}

# Use lapply to generate the heatmaps for different thresholds
heatmaps <- lapply(thresholds, generate_heatmap)


pdf(file = "Heatmap_22_samples_NA_thresholds(20231019).pdf",
    height = 14,
    width = 11)
generate_heatmap(11)
generate_heatmap(10)
generate_heatmap(9)
generate_heatmap(8)
generate_heatmap(7)
generate_heatmap(6)
generate_heatmap(5)
generate_heatmap(4)
generate_heatmap(3)
generate_heatmap(2)
generate_heatmap(1)
generate_heatmap(0)
  dev.off()

# Create the heatmap using pheatmap
pdf(file = "Heatmap_22_samples_fitered_NA(20231019).pdf",
    height = 15,
    width = 12)
p <- pheatmap(
  filtered_heatmap_data,
  clustering_distance_rows = drows_3,  # You can choose a different clustering method if needed
  clustering_distance_cols = "euclidean",  # You can choose a different clustering method if needed
  annotation_row = annotation_df,  # Add the expression ratio annotation
  annotation_colors = list(Expression_Ratio = expression_ratio_palette),
  col = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100),  # Color palette
  na_col = "grey",  # Set the color for NA values
  main = "Gene Expression Heatmap",
  clustering_method = "complete",
  fontsize_row = 4,  # Adjust row label size
  fontsize_col = 9,  # Adjust column label size
  cutree_rows = 2
)
p
dev.off()

object <- generate_heatmap(11)

Some major notes from the Supplementary Material of the EGAD publication: - The neighbor voting method takes a gene network and a set of annotations (gene label vectors) as input, and outputs prediction scores for each of the annotation sets, as well as the performance metric of the network.

Let’s try to bring out the EGAD RMarkdown file and try to run it in this file!

testing[[1]]

           auc avg_node_degree degree_null_auc
[1,] 0.6306979        15228.83       0.6945955

Let’s also test it for the gene list derived based on the NA count filter!

read.csv("NA_filter11_genes.csv",header = 0) %>% pull(V1)

  [1] "CHRDL1"    "MIR223"    "RAP2C−AS1" "FHL1"      "MAP7D3"    "GPC4"      "PHKA2"     "RTL8B"     "BCAP31"   
 [10] "TAB3"      "RLIM"      "SMS"       "MTMR1"     "CD99L2"    "TAF1"      "RPS6KA3"   "TMEM164"   "SOWAHD"   
 [19] "MOSPD1"    "SLC35A2"   "WWC3"      "RBM41"     "RTL8A"     "RTL5"      "ZNF185"    "DLG3"      "MORF4L2"  
 [28] "IL13RA1"   "OTUD5"     "MAGT1"     "GAB3"      "XIAP"      "GK"        "PLXNA3"    "ACOT9"     "TLR8"     
 [37] "GK−AS1"    "UBE2A"     "TFE3"      "SSR4"      "IDS"       "TMEM185A"  "IRAK1"     "ELF4"      "SLC9A7"   
 [46] "TBC1D25"   "BEX4"      "NBDY"      "EBP"       "MAGED2"    "TSPYL2"    "SEPTIN6"   "WDR13"     "SLC25A5"  
 [55] "AIFM1"     "FAM3A"     "PDHA1"     "MAMLD1"    "PRPS2"     "CCDC120"   "GNL3L"     "CRLF2"     "TRMT2B"   
 [64] "STAG2"     "ANOS1"     "UXT"       "CCDC22"    "CYSLTR1"   "PGK1"      "TIMM17B"   "TSIX"      "SLC6A8"   
 [73] "FMR1"      "ZBTB33"    "ITM2A"     "TSPAN7"    "DIPK2B"    "ZNF275"    "GRIPAP1"   "PIGA"      "WDR44"    
 [82] "BCOR"      "WASH6P"    "FAM199X"   "LONRF3"    "TLR7"      "SLC25A43"  "SAT1"      "ARAF"      "HUWE1"    
 [91] "SASH3"     "TMEM187"   "LAMP2"     "FTX"       "SYN1"      "ABCD1"     "MECP2"     "BTK"       "MPP1"     
[100] "ZXDB"      "FLNA"      "INE1"      "TAF9B"     "MBTPS2"    "TSC22D3"   "BRWD3"     "PUDP"      "TMSB4X"   
[109] "ZNF75D"    "ARSD−AS1"  "RPL10"     "NXT2"      "VMA21"     "DOCK11"    "SNX12"    

testing_3[[1]]

           auc avg_node_degree degree_null_auc
[1,] 0.6938566        16374.84       0.7693189

Let’s check the co-expression for genes showing high XCI skew from the heatmap when NA count filter = 11

testing_4[[1]]

           auc avg_node_degree degree_null_auc
[1,] 0.6301435        15577.32       0.7222064

Let’s also compare when we set the read count filter as 5 (NA counts <= 6)

testing_5[[1]]

           auc avg_node_degree degree_null_auc
[1,] 0.7194659        14861.01       0.6549716

As a last point of reference, let’s perform EGAD on all genes in the heatmap for read count filter = 5.

Compare the expression values of these genes from the count matrix!

# Let's check for the expression from mapping_XCI_5
merged_exprs[mapping_XCI_5$ENSEMBL,] %>% 
  rowSums() %>% 
  mean()

[1] 567950.8