comparing_to_other_XIST_systems.R

## ---------------------------
##
## Purpose of script:
##
## Author: Prakhar Bansal
##
## Date Created: 2021-09-09
##
## Copyright (c) Prakhar Bansal, 2021
## Email: pbansal@uchc.edu
##
## ---------------------------
##
## Notes:
##
##
## ---------------------------

## set working directory
setwd(dirname(unlist(rstudioapi::getSourceEditorContext()['path'])))

## ---------------------------

## load up the packages we will need:  (uncomment as required)
require(tidyverse); theme_set(theme_classic() + theme(plot.title=element_text(hjust=0.5)))
require(ggsci)

## ---------------------------

# lawrence lab neurons
lawrence_neural_counts_df <- read_csv("external_tables/GSE125837_monolayer_rawcounts.txt.gz",
                                      col_names = c(
                                        "index", "gene_name","diff10.day21.c5a.d14dox",
                                        "diff10.day21.c5a.d0dox","diff10.day21.c5a.nd",
                                        "diff10.day28.c5a.d14dox","diff10.day28.c5a.d21dox",
                                        "diff10.day28.c5a.d0dox","diff10.day28.c5a_nd",
                                        "diff10.day28.dis.nd","diff10.day28.par.d0dox",
                                        "diff10.day28.par.nd","diff6.day28.par.d0dox",
                                        "diff6.day28.par.nd","diff7.day21.c5a.d14dox",
                                        "diff7.day21.c5a.d0dox","diff7.day21.c5a.nd",
                                        "diff7.day28.c5a.d14dox","diff7.day28.c5a.d21dox",
                                        "diff7.day28.c5a.d0dox","diff7.day28.c5a.nd",
                                        "diff7.day28.dis.nd","diff8.day21.c5a.d14dox",
                                        "diff8.day21.c5a.d0dox","diff8.day21.c5a.nd",
                                        "diff8.day28.c5a.d14dox","diff8.day28.c5a.d21dox",
                                        "diff8.day28.c5a.d0dox","diff8.day28.c5a.nd",
                                        "diff8.day28.dis.nd","diff8.day28.par.d0dox",
                                        "diff8.day28.par.nd"),
                                      skip = 1) |> 
  dplyr::select(-index)
lawrence_neural_samples_df <- read_csv("external_tables/lawrence_neural_sample_info.csv")
library(DESeq2)
dds_all <- DESeqDataSetFromMatrix(
  as.data.frame(lawrence_neural_counts_df),
  colData = lawrence_neural_samples_df |> dplyr::relocate(sample_name, .before = 1),
  design = ~ dox_status + genotype,
  tidy = T
)
dds_all <- estimateSizeFactors(dds_all)
nc <- counts(dds_all, normalized=TRUE)
# only keep genes with more than 5 counts in at least 10% of samples
filter <- rowSums(nc >= 5) >= floor(ncol(nc)/10)
dds_all <- dds_all[filter,]
dds_all <- DESeq(dds_all)

dds_all <- dds_all[which(mcols(dds_all)$betaConv),]

lawrence_neural_vst_df <- assay(vst(dds_all)) |> 
  as.data.frame() |> 
  rownames_to_column(var="gene_name")
write_csv(lawrence_neural_vst_df, "tables/lawrence_neural_vst_df.csv.gz")

library(biomaRt)
mart <- useEnsembl(biomart = "ensembl", 
                   GRCh = 37,
                   dataset = "hsapiens_gene_ensembl")
G_list <- getBM(filters= "hgnc_symbol", 
                attributes= c("ensembl_gene_id_version",
                              "hgnc_symbol", 
                              "chromosome_name", 
                              "start_position", 
                              "end_position", 
                              "strand", 
                              "description",
                              "wikigene_name",
                              "gene_biotype"),
                values=lawrence_neural_counts_df$gene_name,
                mart= mart)


gene_info_df <- as.data.frame(G_list) |> 
  group_by(hgnc_symbol) |> 
  filter(row_number()==1) |> 
  ungroup()


chr21_rownorm_lawrence_neural_foldchange <- lawrence_neural_vst_df |> 
  left_join(gene_info_df |> 
              dplyr::select(hgnc_symbol, chromosome_name, start_position), by=c("gene_name"="hgnc_symbol")) |> 
  pivot_longer(cols = -c("gene_name", "chromosome_name", "start_position"), 
               names_to = "sample_name", 
               values_to = "vst_count") |>
  mutate(sample_name=case_when(sample_name == "diff10.day28.c5a_nd" ~ "diff10.day28.c5a.nd",
                               T ~ sample_name)) |> 
  left_join(lawrence_neural_samples_df |> 
              dplyr::select(sample_name, 
                            dox_status, 
                            genotype, 
                            effective_genotype), 
            by=c("sample_name")) |>
  filter(chromosome_name == "21") |> 
  group_by(gene_name) |> 
  summarise(d21_mean = mean(vst_count[`genotype` == "D21"]),
            t21_mean = mean(vst_count[`genotype` == "T21"]),
            ct21_mean = mean(vst_count[`genotype` == "T21" | (`genotype` == "T21+XIST" & dox_status == "nd")]),
         dox_count = mean(vst_count[`genotype` == "T21+XIST" & dox_status != "nd"]),
         nodox_count = mean(vst_count[`genotype` == "T21+XIST" & dox_status == "nd"]),
         fold_change_xist = dox_count - ct21_mean,
         fold_change_d21=d21_mean - ct21_mean) |> 
  ungroup() #|> 
  #top_frac(.5, abs(fold_change_d21))

output_dir <- "outputs/2024-02-17 21-39-53.641348 timepoint+dox_status+effective_genotype"
neural_output_dir <- (read_csv(paste0(output_dir, "/params.csv.gz")) |>
                        filter(parameter == "neural_output_dir") |>
                        pull(value))[1]
sample_sheet_df <- read_csv(paste0(output_dir, "/tables/sample_sheet_df.csv.gz"))
batch_corrected_vsd <- read_csv(paste0(output_dir, "/tables/vsd_salmon_df.csv.gz"))
neural_sample_sheet_df <- read_csv(paste0(neural_output_dir, "/tables/sample_sheet_df.csv.gz"))
neural_vsd_df <- read_csv(paste0(neural_output_dir, "/tables/vsd_salmon_df.csv.gz"))
gene_info_df <- read_csv("tables/gene_info_df.csv.gz")

our_chr21_rownorm_vsd <- batch_corrected_vsd |> 
  pivot_longer(cols = -ensembl_id, names_to = "sample_id", values_to = "vst_count") |> 
  left_join(gene_info_df, by=c("ensembl_id"="ensembl_gene_id_version")) |> 
  left_join(sample_sheet_df, by=c("sample_id")) |> 
  filter(chromosome_name == "21" &
           hgnc_symbol %in% chr21_rownorm_lawrence_neural_foldchange$gene_name) |> 
  group_by(hgnc_symbol) |> 
  summarise(dox_mean= mean(vst_count[(genotype == "T21XIST") & (dox_status == "yes")]),
            nodox_mean=mean(vst_count[(genotype == "T21XIST") &( dox_status == "no")]),
            withdrawal_mean=mean(vst_count[(genotype == "T21XIST") & (dox_status == "withdrawal")]),
            dwd_mean=mean(vst_count[(genotype == "T21XIST") & (dox_status != "no")]),
            d21_mean=mean(vst_count[(genotype == "D21")]),
            t21_mean=mean(vst_count[(genotype == "T21")]),
            ct21_mean=mean(vst_count[(genotype == "T21") | (genotype == "T21XIST" & dox_status == "no")]),
            fold_change_xist=dwd_mean-ct21_mean,
            fold_change_wd=withdrawal_mean-nodox_mean,
            fold_change_d21=d21_mean - ct21_mean) #|> 
  # top_frac(.5, abs(fold_change_d21))

our_neural_chr21_rownorm_vsd <- neural_vsd_df |> 
  pivot_longer(cols = -ensembl_id, names_to = "sample_id", values_to = "vst_count") |> 
  left_join(gene_info_df, by=c("ensembl_id"="ensembl_gene_id_version")) |> 
  left_join(neural_sample_sheet_df, by=c("sample_id")) |> 
  filter(chromosome_name == "21" &
           hgnc_symbol %in% chr21_rownorm_lawrence_neural_foldchange$gene_name) |> 
  group_by(hgnc_symbol) |> 
  summarise(dwd_mean= mean(vst_count[(genotype == "T21XIST") & (dox_status != "no")]),
            # withdrawal_mean=mean(vst_count[(genotype == "T21XIST") & (dox_status == "wd")]),
            d21_mean=mean(vst_count[(genotype == "D21")]),
            t21_mean=mean(vst_count[(genotype == "T21")]),
            ct21_mean=mean(vst_count[(genotype == "T21") | (genotype == "T21XIST" & dox_status == "no")]),
            fold_change_xist=dwd_mean-ct21_mean,
            # fold_change_wd=withdrawal_mean-t21_mean,
            fold_change_d21=d21_mean - ct21_mean) #|> 
  # top_frac(.5, abs(fold_change_d21))

comparison_df <- chr21_rownorm_lawrence_neural_foldchange |> 
  dplyr::select(hgnc_symbol=gene_name, fold_change_xist, fold_change_d21) |> 
  pivot_longer(cols = -hgnc_symbol,
               names_to = "comparison",
               values_to = "vst_diff") |> 
  mutate(study="Czermiński et al. Neural") |> 
    rbind(
      our_chr21_rownorm_vsd |> 
        dplyr::select(hgnc_symbol, fold_change_xist, fold_change_d21) |> 
        pivot_longer(cols = -hgnc_symbol,
                     names_to = "comparison",
                     values_to = "vst_diff") |> 
        mutate(study="This Study iPSCs")
    ) |> 
  rbind(
    our_neural_chr21_rownorm_vsd |> 
      dplyr::select(hgnc_symbol, fold_change_xist, fold_change_d21) |> 
      pivot_longer(cols = -hgnc_symbol,
                   names_to = "comparison",
                   values_to = "vst_diff") |> 
      mutate(study="This Study Neural")
  )
  
ggplot(comparison_df, aes(x=vst_diff, color=comparison)) +
  geom_density() +
  scale_color_aaas(limits=c("fold_change_d21", "fold_change_xist"), labels=c("D21 vs. T21", "T21XIST vs T21")) +
  facet_grid(rows = vars(study), scales = "free") +
  geom_vline(xintercept = log2(2/3), linetype=2, size=.3) +
  geom_vline(xintercept = 0, size=.5) +
  coord_cartesian(xlim = c(-1,1)) +
  labs(title="Dosage Restoration Systems",
       x="VST Difference to T21",
       color="Comparison") +
  theme(legend.position = c(.75, .85),
        legend.background = element_rect(color="black"))
ggsave("output_figs/2024-02-17 21-39-53.641348 timepoint+dox_status+effective_genotype/system_comparison_same_genes.pdf", width=3.5, height=3.1)