<jupyter_start><jupyter_text>NGRE and GRE data Authors:<NAME>  <NAME>  Libraries & Paths<jupyter_code>library(readxl)
library(ggplot2)
library(gridExtra)
library(grid)
library(biomaRt)
options(repr.plot.width=10, repr.plot.height=7)

path_ngre_human = '../human/ngre/'
path_gre_human = '../human/gre/'
path_ngre_mouse = '../mouse/ngre/'
path_gre_mouse = '../mouse/gre/'

path_results = '../results/'



path_human = '../human/'
path_mouse = '../mouse/'

path_gene_annot = '../gene_annotation/'
path_ensembl = '../ensembl_annotation/'

target_stats = c('gre', 'ngre')

# Our tags only
target_tags = c('H3K4me3', 'H3K27me3')

# One additional tag for heterochromatin
target_tags = c('H3K4me3', 'H3K27me3', 'H3K9me3')

# With additional tags from Dasha: H3K9ac и H3K27ac 
target_tags = c('H3K4me3', 'H3K9ac', 'H3K27ac', 'H3K27me3', 'H3K9me3')<jupyter_output><empty_output><jupyter_text>Read NGRE and GRE data for Human<jupyter_code>stat_all_h = c()  # Table to accumulate all statistics
for(stat in target_stats)  # GRE or NGRE
{
    path_stat = paste0(c(path_human, stat, '/'), collapse = '')
    files = list.files(path_stat)
    
    for(tag in target_tags) 
    {   
        acc_stat = 0  # Accumulate GRE or NGRE, the mean across all cell lines
        for(file in files)
        {
            filename = paste0(path_stat, file)

            data = as.data.frame(read_excel(filename, sheet = tag))
            gene_names = data[,'Gene_name']
            acc_stat = acc_stat + data[,'total_signal']
        }
        acc_stat = acc_stat / length(files)
        length(gene_names)
        length(acc_stat)
        names(acc_stat) <- gene_names
        
        file_acc_stat = paste0(c(path_results, 'human_', stat, '_', tag, '.txt'), collapse = '')

        write.table(file = file_acc_stat, x = acc_stat, quote = FALSE, col.names = FALSE)
        
        # Accumulate all statistics
        stat_all_h = cbind(stat_all_h, acc_stat)
        colnames(stat_all_h)[ncol(stat_all_h)] <- paste0(c(stat, tag), collapse = '_')
    }  
}
head(stat_all_h)
# Save backgrounf info for DAVID annotation
write.table(paste0(path_gene_annot, 'background_genes_human.txt'), 
            x=gene_names, quote = FALSE, row.names = FALSE, col.names = 'human_background')<jupyter_output><empty_output><jupyter_text>Read NGRE and GRE data for Mouse<jupyter_code>stat_all_m = c() 
for(tag in target_tags)
{
    for(stat in target_stats)  # GRE or NGRE
    {
        path_stat = paste0(c(path_mouse, stat, '/'), collapse = '')
        files = list.files(path_stat, pattern = paste0(c('*', tag, '.txt'), collapse=""))
        
        acc_stat = 0  # Accumulate GRE or NGRE, the mean across all cell lines
        for(file in files)
        {
            filename = paste0(path_stat, file)

            data = as.data.frame(read.table(filename,header  = TRUE))
            gene_names = data[,'gene']
            acc_stat = acc_stat + data[,'total_signal']
        }
        acc_stat = acc_stat/length(files)
        length(gene_names)
        length(acc_stat)
        names(acc_stat) <- gene_names

        file_acc_stat = paste0(c(path_results, 'mouse_', stat, '_', tag, '.txt'), collapse = '')

        write.table(file = file_acc_stat, x = acc_stat, quote = FALSE, col.names = FALSE)
        
        # Accumulate all statistics
        stat_all_m = cbind(stat_all_m, acc_stat)
        colnames(stat_all_m)[ncol(stat_all_m)] <- paste0(c(stat, tag), collapse = '_')
    }  
}
head(stat_all_m)
# Save backgrounf info for DAVID annotation
write.table(paste0(path_gene_annot, 'background_genes_mouse.txt'), 
            x=gene_names, quote = FALSE, row.names = FALSE, col.names = 'mouse_background')<jupyter_output><empty_output><jupyter_text>Create Plots for "GRE vs NGRE" for Human & Mouse<jupyter_code>stats = list(Human = stat_all_h, Mouse = stat_all_m)
names(stats)

plot_list = list()
for(istat in 1:length(stats))
{
    stat_all = stats[[istat]]
    for(tag in target_tags) 
    {
        # Get profiles
        gre = stat_all[,paste0('gre_', tag)]
        ngre = stat_all[,paste0('ngre_', tag)]

        # Fit Linear Regressio with zero intercept
        slope = coef(lm(ngre ~ gre - 1))

        # Find outliers
        cutoff = 0.95
        point_dist = abs(-slope * gre + ngre) / slope
        idx_dist_up = ngre > slope * gre
        idx_dist_down = ngre < slope * gre

        point_dist_up = point_dist[idx_dist_up]
        point_dist_down = point_dist[idx_dist_down]

        threshold_up = quantile(point_dist_up, cutoff)
        threshold_down = quantile(point_dist_down, cutoff)


        outliers_up = idx_dist_up & (point_dist > threshold_up)
        outliers_down = idx_dist_down & (point_dist > threshold_down)
        normal_points = !outliers_up & !outliers_down

        # Calculate correlation and p-value and create annotation
        res = cor.test(gre, ngre)[c(3,4)]
        res_pval = as.character(round(res[[1]],2))
        res_cor = as.character(round(res[[2]],2))
        text_cor = paste('cor: ', res_cor)
        text_pval = paste('p-value: ', res_pval)
        annotations <- data.frame(xpos = Inf, ypos =  Inf,
                                  annotateText = paste0(c(text_cor, text_pval), collapse = '\n'),
                                  hjustvar = c(1), vjustvar = c(1)) #<- adjust


        # Create plot
        p <- ggplot(as.data.frame(cbind(gre, ngre)), aes(gre, ngre)) + 
             geom_point(alpha=0.1, color = 2*outliers_up + 3*outliers_down + normal_points)+ 
             geom_abline(intercept = 0, slope = slope, colour = '#858585') +
             geom_text(data=annotations,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText))+
             ggtitle(paste0(c(names(stats)[istat], tag), collapse = ' ')) + labs(x='GRE', y='NGRE')
        plot_list[[length(plot_list) + 1]] = p
        
        # Save genes for the further DAVID annotation
        file_up = paste0(c(path_gene_annot, tolower(names(stats)[istat]), '_', tag, '_up.txt'), collapse = '')
        file_down = paste0(c(path_gene_annot, tolower(names(stats)[istat]), '_', tag, '_down.txt'), collapse = '')
        
        gene_names = rownames(stat_all)
        
        write.table(file = file_up, x = gene_names[outliers_up], 
                    quote = FALSE, col.names = FALSE, row.names = FALSE)
        write.table(file = file_down, x = gene_names[outliers_down], 
                    quote = FALSE, col.names = FALSE, row.names = FALSE)
    }
}<jupyter_output><empty_output><jupyter_text>Plotting<jupyter_code>options(repr.plot.width=10, repr.plot.height=5)

n <- length(plot_list)
nCol <- length(target_tags)
do.call("grid.arrange", c(plot_list, ncol=nCol))<jupyter_output><empty_output><jupyter_text>Convert Gene_symbol to ENSG<jupyter_code>files = list.files(path_gene_annot)

ensembl <- useMart('ENSEMBL_MART_ENSEMBL', 
                   dataset = 'hsapiens_gene_ensembl', 
                   host = 'may2009.archive.ensembl.org')

for(i in 1:length(files))
{
    file_in = paste0(path_gene_annot, files[i])
    file_out = paste0(path_ensembl, files[i])
    gene_symbols = read.table(file_in, header = FALSE)

    gene_annot = getBM(attributes=c('hgnc_symbol', 'ensembl_gene_id'), 
                       filters = 'hgnc_symbol', 
                       values= gene_symbols, 
                       mart= ensembl)
    if(grepl('background', files[i]))
        write.table(file = file_out, gene_annot[,2], quote = FALSE, row.names = FALSE, col.names = 'background')
    else
        write.table(file = file_out, gene_annot[,2], quote = FALSE, row.names = FALSE, col.names = FALSE)
}<jupyter_output>Batch submitting query [==>----------------------------]   9% eta:  2s<jupyter_text>Prepare files for multi-list upload to DAVID<jupyter_code>organisms = c('human', 'mouse')
for(org in organisms)
{
    files_genes = list.files(path_ensembl, pattern = paste0(org, '_H*'))
    genes = c()
    for(file in files_genes)
    {
        tmp = read.table(paste0(path_ensembl, file))
        if (is.null(nrow(genes)))
        {
            genes = tmp
            next
        }
        
        n = max(nrow(genes), nrow(tmp))
        n_add_genes = n - nrow(genes)
        genes_add = matrix(rep('', n_add_genes * ncol(genes)), 
                                nrow = n_add_genes, ncol = ncol(genes))
        colnames(genes_add) = colnames(genes)
        genes = rbind(genes, genes_add)

        n_add_tmp = n - nrow(tmp)
        tmp = rbind(tmp, matrix(rep('', n_add_tmp)))

        genes = cbind(genes, tmp)
    }
    colnames(genes) <- files_genes
    
    file_multi_list = paste0(c(path_ensembl,'multi_list_',  org, '.txt'), collapse = '')
    write.table(file = file_multi_list, genes, row.names = FALSE, quote = FALSE, sep = '\t')
}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13529210796036897, 'fraction_numerical': 0.011615989067304407, 'mean_word_length': 2.833260584897425, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '618983', 'license': 'MIT', 'n_tokens_mistral': 3191, 'n_tokens_neox': 3013, 'n_words': 591, 'path': 'scripts/analyse_genes.ipynb', 'repo_name': 'arsen-z/evo_epigen', 'size': '8804'}