<jupyter_start><jupyter_text>Annotate genes-outliers with DAVID Authors:<NAME>  <NAME><jupyter_code># # To install GO.db package for R 3.5
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("GO.db", version = "3.8")
# -------------------------------------------------------
library('metap')
suppressMessages(library('GO.db'))
path_class_daniil = '../go_daniil/'
path_ensembl = '../ensembl_annotation/'
path_david_clusters = '../clusters_david/'


file_david_terms = 'david_all_go.txt'

target_tags = c('H3K4me3', 'H3K9ac', 'H3K27ac', 'H3K27me3', 'H3K9me3') 
organisms = c('human', 'mouse')

bad_list = c('general_terms.txt', 'other.txt', 'other_metabolism.txt')<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code># -------------------------------------------------------
# The function returns ancestrors for a GO-term
# -------------------------------------------------------
get_anc_by_term <-function(term)
{
    ontology = Ontology(term)
    if (is.na(ontology)){ return (c())}

    if (ontology == 'CC') {GOxxANCESTOR = GOCCANCESTOR}
    if (ontology == 'BP') {GOxxANCESTOR = GOBPANCESTOR}
    if (ontology == 'MF') {GOxxANCESTOR = GOMFANCESTOR}

    ancestors = get(term, GOxxANCESTOR)
    ancestors = ancestors[ancestors!='all']
    return (ancestors)
}

# -------------------------------------------------------
# This function returns a probable category for the term 
# if this term is in c('general_terms.txt', 'other.txt', 'other_metabolism.txt').
# -------------------------------------------------------
get_class_by_term <- function(term, ancestors_of_terms, classes_of_terms)
{
    bad_list = c('general_terms.txt', 'other.txt', 'other_metabolism.txt')
    if(!(term %in% bad_list)){return (NULL)}
    
    ancestors = get_anc_by_term(term)
    if(is.null(ancestors)){return (NULL)}
    x = intersect(ancestors, names(ancestors_of_terms))
    if (length(x) == 0) {return (NULL)}
    
    y = sapply(x, function(y) classes_of_terms[[y]])
    y = setdiff(y, c('general_terms.txt', 'other.txt', 'other_metabolism.txt'))
    if (length(y) == 0) {return (NULL)}
               
    y = table(y)
    res_class = names(y)
    tmp = res_class[y == max(y)]
    return(tmp)
}  
# -------------------------------------------------------
# This function tries to find the most similat term with by common ancestors
# -------------------------------------------------------
get_class_by_similar_anc <- function(term, ancestors_of_terms, classes_of_terms)
{
    anc = get_anc_by_term(term)
    if(is.null(anc)){return('other.txt')}
    bad_list = c('general_terms.txt', 'other.txt', 'other_metabolism.txt')
    similar_term = 'other.txt'
    n_similar = 0
    for(term2 in names(ancestors_of_terms))
    {
        if(classes_of_terms[term2] %in% bad_list){next}
        tmp_intersect = intersect(anc, ancestors_of_terms[term2])
        n_common = length(tmp_intersect)
        if(n_common == 0) {next}
        
        n_common = length(tmp_intersect[[1]])
        if(n_common > n_similar)
        {
            n_similar = n_common
            similar_term = term2
        }
    }
    if(n_similar == 0) {return('other.txt')}
    return(classes_of_terms[similar_term])
}<jupyter_output><empty_output><jupyter_text>Read all Daniil's annotation<jupyter_code>class_daniil = list.files(path_class_daniil)
ancestors_of_terms = list()
classes_of_terms = list()
for(file in class_daniil)
{
#     print(file)
    terms = read.table(paste0(path_class_daniil, file), header = FALSE, sep = '\t')[,1]
#     print(length(terms))
    for(t in terms)
    {
        term = sprintf("GO:%07d", as.numeric(substr(t, 4, nchar(t))))

        ancestors_of_terms[[term]] = c(get_anc_by_term(term), term)
        classes_of_terms[[term]] = file
    } 
}
length(ancestors_of_terms)
length(classes_of_terms)<jupyter_output>Warning message in scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
“EOF внутри закавыченной строки”Warning message in read.table(paste0(path_class_daniil, file), header = FALSE, sep = "\t"):
“incomplete final line found by readTableHeader on '../go_daniil/electron_transfer_chain_reactions.txt'”<jupyter_text>Read all clustering from DAVID and atribute new terms to existing caregories<jupyter_code>ancestors_of_terms_init = ancestors_of_terms
classes_of_terms_init = classes_of_terms
# Annotation which Anna performed on the David clustering
additional_annot = read.table('additional_annot.txt', header = FALSE, sep=c('\t'), quote = "\"", row.names=1)
# head(additional_annot)

length(ancestors_of_terms_init)
length(classes_of_terms_init)
ancestors_of_terms = ancestors_of_terms_init
classes_of_terms = classes_of_terms_init
length(classes_of_terms)
length(ancestors_of_terms)

david_cluster_id = 1
files_clusters = list.files(path_david_clusters)
for(i in 1:length(files_clusters))
{   
#     if(i == 12){break}
    if(length(classes_of_terms) != length(ancestors_of_terms)) {break}
    if(i %in% c(0)) {next}
    file = files_clusters[i]
    f = paste0(path_david_clusters, file)
    d_tmp = read.table(text = gsub("~", "\t", readLines(f)), skip=2, header = FALSE, sep=c('\t'), quote = "\"")
    d_go = d_tmp[,2]
    d_descriptor = d_tmp[,3]
    d_descriptor
    
    classes_tmp = c()
    for(term in d_go)
    {
        if(term %in% names(classes_of_terms))
            {classes_tmp = c(classes_tmp, classes_of_terms[[term]])}
        else
            {classes_tmp = c(classes_tmp, 'None')}
    }
    classes_res = unique(setdiff(classes_tmp, c('general_terms.txt', 'other.txt', 'None')))
#     print(length(classes_res))
    if(length(classes_res) < 1)
    {
        for(term in d_go)
        {
            ancestors_of_terms[[term]] = c(get_anc_by_term(term), term)
            classes_of_terms[[term]] = paste0('david_cluster_', as.character(david_cluster_id))
        } 
        david_cluster_id = david_cluster_id + 1
        next
    }
    
    if(length(intersect(classes_tmp, c('general_terms.txt', 'other.txt', 'None'))) == 0){ next}
    
    if(length(classes_res) == 1)
    {
        for(term in d_go)
        {
            ancestors_of_terms[[term]] = c(get_anc_by_term(term), term)
            classes_of_terms[[term]] = classes_res
        } 
        next
    }
    # For Annotation that Anna made
    for(term in d_go)
    {
        if(term %in% names(classes_of_terms)){next}
        ancestors_of_terms[[term]] = c(get_anc_by_term(term), term)
        classes_of_terms[[term]] = as.character(additional_annot[term, 2])
    } 
}
length(classes_of_terms)
length(ancestors_of_terms)
# ------------------------------------
# # previously helpful
# ------------------------------------
# classes_tmp = c()
# for(term in d_go)
# {
#     if(term %in% names(classes_of_terms))
#         {classes_tmp = c(classes_tmp, classes_of_terms[[term]])}
#     else
#         {classes_tmp = c(classes_tmp, 'None')}
# }
# cbind(as.character(d_go),as.character( d_descriptor), classes_tmp)
# # d_descriptor
# write.table(file = 'tmp.txt', x=cbind(as.character(d_go),as.character( d_descriptor), classes_tmp), 
#             quote = FALSE, row.names=FALSE, col.names=FALSE, sep='\t')<jupyter_output><empty_output><jupyter_text>Try to attribute general_terms.txt and others.txt to an informative category<jupyter_code># terms_daniil = names(ancestors_of_terms)
# for(term in terms_daniil)
# {
#     if(!(classes_of_terms[[term]] %in% c('general_terms.txt', 'other.txt', 'other_metabolism.txt')))
#         {next}
#     class_tmp = get_class_by_term(term, ancestors_of_terms, classes_of_terms)
#     if( is.null(class_tmp)){ next }
#     print('ok')
#     print(classes_of_terms[[term]], class_tmp)
#     classes_of_terms[[term]] = class_tmp
# }<jupyter_output><empty_output><jupyter_text>Attribute DAVID terms to existing ones by the intersection of ancestrals<jupyter_code>david_data = read.table(paste0(path_ensembl, file_david_terms), header = FALSE, sep='\t', quote = "\"")
david_terms = as.character(david_data[,2])
# head(david_terms)

new_class_flag = c()
res = c()
for(term in david_terms)
{
    if(term %in% names(ancestors_of_terms))
    {
        res = c(res, classes_of_terms[[term]])
        new_class_flag = c(new_class_flag, 0)
    }else
    {
        class_tmp = get_class_by_term(term, ancestors_of_terms, classes_of_terms)
        
        if(is.null(class_tmp))
        {
#             class_tmp = get_class_by_similar_anc(term, ancestors_of_terms, classes_of_terms)
#             res = c(res, class_tmp)
            res = c(res, 'other.txt')
#             print(class_tmp)
        }else
        {
#             print(class_tmp)
            res = c(res, paste0(class_tmp[1], collapse = '_'))
        }
        new_class_flag = c(new_class_flag, 1)
    }
}

z = cbind(david_data, res, new_class_flag)
write.table(file = 'david_anna_annot.txt', x=z, quote = FALSE, 
            col.names = FALSE, row.names = FALSE, sep='\t')

result_david = as.matrix(res)
row.names(result_david) <- david_terms
head(result_david)<jupyter_output><empty_output><jupyter_text>Table with results text<jupyter_code>name_classes = unique(classes_of_terms)
results_pval = matrix(, nrow = 0, ncol = length(name_classes))
results_count = matrix(, nrow = 0, ncol = length(name_classes))
colnames(results_pval) <- name_classes
colnames(results_count) <- name_classes

for(tag in target_tags)
{
    for(outlier in c('up', 'down'))
    {
        for(organism in organisms)
        {
            results_pval = rbind(results_pval, rep(0, length(name_classes)))
            n = nrow(results_pval)
            pattern = paste0(c(organism, '_', tag, '_', outlier), collapse = '')
            row.names(results_pval)[n] <- pattern
            
            results_count = rbind(results_count, rep(0, length(name_classes)))
            n = nrow(results_count)
            pattern = paste0(c(organism, '_', tag, '_', outlier), collapse = '')
            row.names(results_count)[n] <- pattern
            
            file = paste0(c(path_ensembl, 'david_', pattern, '.txt'), collapse = '')
            d = read.table(file, header = TRUE, sep=c('\t'), quote = "\"")
            d_category = d[,'Term']
            d_pvalue = d[,'PValue']
            
            id = sapply(d_category, function(x) (grepl('~', x)) & (substring(x, 1, 1) == 'G'))
            d_category = as.character(d_category[id])
            d_pvalue = d_pvalue[id]
            d_category = sapply(d_category, function(x) gsub( "~.*$", "", x ))
                                
            d_class = result_david[d_category,]
            unique_class = unique(d_class)
            for(class_tmp in unique_class)
            {
                if(class_tmp %in% bad_list) {next}
                if(is.na(class_tmp)) {next}
                p_vals = d_pvalue[d_class %in% class_tmp]
                if(length(p_vals) > 1)
                {
#                     results[n, class_tmp] = length(p_vals)
                    results_pval[n, class_tmp] = allmetap(p_vals, method = 'sumlog')[,'p'][[1]]
                } else if(length(p_vals) == 1){
                    results_pval[n, class_tmp] = p_vals
                }
                results_count[n, class_tmp] = length(p_vals)

            }
        }
    }
}
write.table(file = 'results_pval.txt', x = cbind(sample = rownames(results_pval), results_pval), 
            quote = FALSE, sep='\t', col.names = TRUE, row.names = FALSE)
write.table(file = 'results_count.txt', x = cbind(sample = rownames(results_count), results_count), 
            quote = FALSE, sep='\t', col.names = TRUE, row.names = FALSE)
p_vals = d_pvalue[d_class %in% class_tmp]
                if(length(p_vals) > 1)
                {
#                     results[n, class_tmp] = length(p_vals)
                    results_pval[n, class_tmp] = allmetap(p_vals, method = 'sumlog')[,'p'][[1]]
                } else if(length(p_vals) == 1){
                    results_pval[n, class_tmp] = p_vals
                }
results_pval[n, ]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1724137931034483, 'fraction_numerical': 0.005982550893228085, 'mean_word_length': 2.7483649953285583, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '638291', 'license': 'MIT', 'n_tokens_mistral': 4282, 'n_tokens_neox': 4126, 'n_words': 730, 'path': 'scripts/annotate_genes.ipynb', 'repo_name': 'arsen-z/evo_epigen', 'size': '12058'}