SQLite Database - Results

library(RSQLite)
library(DBI)
library(dplyr)
library(ggplot2)

#misc code

#close any open SQLite connections library(DBI) library(RSQLite) lapply(dbListConnections(SQLite()), dbDisconnect)

#check if old database exists if (file.exists(“rnaseq_analysis_results.db”)) { cat(“Database file exists”) cat(“File size:”, file.size(“rnaseq_analysis_results.db”), “bytes”) cat(“Last modified:”, file.mtime(“rnaseq_analysis_results.db”), “”) } else { cat(“Database file does not exist”) }

#delete the old database file.remove(“rnaseq_analysis_results.db”) cat(“Old database deleted”)

#verify cat(“File exists after deletion:”, file.exists(“rnaseq_analysis_results.db”), “”)

#after restarting R, delete file.remove(“rnaseq_analysis_results.db”) file.exists(“rnaseq_analysis_results.db”) # Should be FALSE

#create database schema function 
create_database_schema <- function(db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  
  #table 1: statistical analysis with gene annotations
  dbExecute(con, "
    CREATE TABLE IF NOT EXISTS differential_expression (
      dataset_id TEXT,
      comparison_name TEXT,
      gene_id TEXT,
      gene_symbol TEXT,
      log2FoldChange REAL,
      pvalue REAL,
      padj REAL,
      baseMean REAL,
      lfcSE REAL,
      stat REAL,
      PRIMARY KEY (dataset_id, comparison_name, gene_id)
    )
  ")
  
  #table 2: normalized expression for whole dataset
  dbExecute(con, "
    CREATE TABLE IF NOT EXISTS normalized_expression (
      sample_id TEXT,
      gene_id TEXT,
      expression_value REAL,
      PRIMARY KEY (sample_id, gene_id)
    )
  ")
  
  #table with extra information
  dbExecute(con, "
    CREATE TABLE IF NOT EXISTS datasets (
      dataset_id TEXT PRIMARY KEY,
      dataset_name TEXT,
      description TEXT,
      organism TEXT,
      n_samples INTEGER,
      n_genes INTEGER,
      date_added DATE,
      study_design TEXT,
      geo_id TEXT
    )
  ")
  
  dbExecute(con, "
    CREATE TABLE IF NOT EXISTS samples (
      sample_id TEXT PRIMARY KEY,
      dataset_id TEXT,
      original_sample_name TEXT,
      condition TEXT,
      treatment TEXT,
      time_point TEXT,
      batch TEXT,
      cell_line TEXT,
      tissue_type TEXT,
      FOREIGN KEY (dataset_id) REFERENCES datasets (dataset_id)
    )
  ")
  
  #indices for fast querying
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_de_gene_id ON differential_expression (gene_id)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_de_gene_symbol ON differential_expression (gene_symbol)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_de_dataset ON differential_expression (dataset_id)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_de_padj ON differential_expression (padj)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_de_lfc ON differential_expression (log2FoldChange)")
  
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_norm_gene ON normalized_expression (gene_id)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_norm_sample ON normalized_expression (sample_id)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_norm_expr ON normalized_expression (expression_value)")
  
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_samples_condition ON samples (condition)")
  dbExecute(con, "CREATE INDEX IF NOT EXISTS idx_samples_dataset ON samples (dataset_id)")
  
  #creating this for common queries
  dbExecute(con, "
    CREATE VIEW IF NOT EXISTS significant_genes AS
    SELECT dataset_id, comparison_name, gene_id, gene_symbol, 
           log2FoldChange, padj, baseMean
    FROM differential_expression 
    WHERE padj < 0.05 AND ABS(log2FoldChange) > 1
  ")
  
  dbExecute(con, "
    CREATE VIEW IF NOT EXISTS gene_expression_summary AS
    SELECT 
      ne.gene_id,
      COUNT(DISTINCT ne.sample_id) as n_samples,
      AVG(ne.expression_value) as mean_expression,
      MIN(ne.expression_value) as min_expression,
      MAX(ne.expression_value) as max_expression,
      (MAX(ne.expression_value) - MIN(ne.expression_value)) as expression_range
    FROM normalized_expression ne
    GROUP BY ne.gene_id
  ")
  
  dbDisconnect(con)
  message("Database schema created with advanced features!")
}

#create the fresh database
create_database_schema()

#verification
con <- dbConnect(SQLite(), "rnaseq_analysis_results.db")
message("Tables created:")
print(dbListTables(con))
dbDisconnect(con)

process_dataset <- function(dataset_id, final_results, vst_data, comparison_name, 
                           sample_metadata = NULL, db_path = "rnaseq_analysis_results.db") {
  
  con <- dbConnect(SQLite(), db_path)
  
  # DUPLICATE PREVENTION: Check if this exact combination already exists
  existing_de <- dbGetQuery(con, paste0("
    SELECT COUNT(*) as count 
    FROM differential_expression 
    WHERE dataset_id = '", dataset_id, "' AND comparison_name = '", comparison_name, "'
  "))
  
  if (existing_de$count > 0) {
    message("Dataset ", dataset_id, " with comparison ", comparison_name, " already exists. Skipping...")
    dbDisconnect(con)
    return()
  }
  
  message(paste("Processing dataset:", dataset_id, "-", comparison_name))
  
  #1st: Statistical analysis results - Table 1
  
  # reading the csv
  de_results <- read.csv(final_results, stringsAsFactors = FALSE, row.names = 1)
  
  # need to convert the rownames to a separate column
  de_results$gene_id <- rownames(de_results)
  
  # adding required database columns
  de_results$dataset_id <- dataset_id
  de_results$comparison_name <- comparison_name
  
  # handle missing columns 
  required_cols <- c("dataset_id", "comparison_name", "gene_id", "gene_symbol", 
                    "log2FoldChange", "pvalue", "padj", "baseMean")
  optional_cols <- c("lfcSE", "stat")
  
  for (col in optional_cols) {
    if (!col %in% names(de_results)) {
      de_results[[col]] <- NA
    }
  }
  
  # handle missing gene_symbol column
  if (!"gene_symbol" %in% names(de_results)) {
    de_results$gene_symbol <- NA
  }
  
  # select and clean data
  de_table <- de_results[, c(required_cols, optional_cols)]
  de_table <- de_table[!is.na(de_table$gene_id), ]
  
  # insert statistical results
  dbWriteTable(con, "differential_expression", de_table, append = TRUE, row.names = FALSE)
  message(paste("  Added", nrow(de_table), "differential expression results"))
  
  #2nd: Normalized expression data - Table 2
  
  # read the long-format VST data 
  norm_expr <- read.csv(vst_data, stringsAsFactors = FALSE)
  
  # Handle the X column issue - remove row numbers column
  if ("X" %in% names(norm_expr)) {
    norm_expr <- norm_expr[, !names(norm_expr) %in% "X"]
  }
  
  # Clean expression data
  norm_table <- norm_expr[!is.na(norm_expr$expression_value), ]
  
  # DUPLICATE PREVENTION: Remove any existing expression data from these samples
  sample_ids <- unique(norm_table$sample_id)
  if (length(sample_ids) > 0) {
    sample_list <- paste0("'", sample_ids, "'", collapse = ",")
    deleted_rows <- dbExecute(con, paste0("DELETE FROM normalized_expression WHERE sample_id IN (", sample_list, ")"))
    if (deleted_rows > 0) {
      message("  Removed ", deleted_rows, " existing expression values to prevent duplicates")
    }
  }
  
  # Insert expression data in chunks
  chunk_size <- 10000
  total_rows <- nrow(norm_table)
  for (i in seq(1, total_rows, chunk_size)) {
    end_i <- min(i + chunk_size - 1, total_rows)
    chunk <- norm_table[i:end_i, ]
    dbWriteTable(con, "normalized_expression", chunk, append = TRUE, row.names = FALSE)
    
    # Progress indicator for large datasets
    if (total_rows > 50000) {
      message("  Processed ", end_i, "/", total_rows, " expression values")
    }
  }
  
  message(paste("  Added", total_rows, "normalized expression values"))
  
  #3rd: Dataset metadata 
  
  # Check if dataset metadata already exists
  existing_dataset <- dbGetQuery(con, paste0("SELECT COUNT(*) as count FROM datasets WHERE dataset_id = '", dataset_id, "'"))
  
  if (existing_dataset$count == 0) {
    n_samples <- length(unique(norm_table$sample_id))
    n_genes <- length(unique(de_table$gene_id))
    
    # GEO ID mapping for each dataset
  geo_ids <- c("dataset1" = "GSE243564", 
               "dataset2" = "GSE130160",
               "dataset3" = "GSE129221", 
               "dataset4" = "GSE94405",
               "dataset5" = "GSE79688")
  
  dataset_num <- substr(dataset_id, nchar(dataset_id), nchar(dataset_id))
    
    dataset_record <- data.frame(
      dataset_id = dataset_id,
      description = comparison_name,
      organism = "Homo sapiens",
      n_samples = n_samples,
      n_genes = n_genes,
      date_added = as.character(Sys.Date()),
      study_design = comparison_name,
      geo_id = geo_ids[dataset_id],
      stringsAsFactors = FALSE
    )
    
    dbWriteTable(con, "datasets", dataset_record, append = TRUE, row.names = FALSE)
    message("  Added dataset metadata")
  } else {
    message("  Dataset metadata already exists, skipping")
  }
  
  dbDisconnect(con)
  message(paste("  Dataset", dataset_id, "processing complete!"))
}

#Dataset 1 - single comparison
process_dataset(
  dataset_id = "dataset1",
  final_results = "D1_results.csv",
  vst_data = "D1_vst_data.csv",
  comparison_name = "untreated_vs_osimertinib"
)

#Dataset 2 - main comparison (wild type vs exon19 deletion)
process_dataset(
  dataset_id = "dataset2",
  final_results = "D2_results.csv",
  vst_data = "D2_vst_data.csv",
  comparison_name = "wildtype_vs_exon19del"
)

#Dataset 3 - single comparison
process_dataset(
  dataset_id = "dataset3",
  final_results = "D3_results.csv",
  vst_data = "D3_vst_data.csv",
  comparison_name = "control_vs_apatinib"
)

#Dataset 4 - single comparison
process_dataset(
  dataset_id = "dataset4",
  final_results = "D4_results.csv",
  vst_data = "D4_vst_data.csv",
  comparison_name = "resistant_vs_parental"
)

#Dataset 5 - main comparison (treatment effect)
process_dataset(
  dataset_id = "dataset5",
  final_results = "D5_results.csv",
  vst_data = "D5_vst_data.csv",
  comparison_name = "gefitinib_vs_dmso"
)

#database verification

con <- dbConnect(SQLite(), "rnaseq_analysis_results.db")

message("=== DATABASE STATUS ===")
message("File size: ", round(file.size("rnaseq_analysis_results.db")/1024/1024, 2), " MB")
message("Tables in database:")
print(dbListTables(con))

message("\n=== DATASETS LOADED ===")
datasets_info <- dbGetQuery(con, "SELECT dataset_id, dataset_name, n_samples, n_genes, study_design FROM datasets")
print(datasets_info)

message("\n=== COMPARISONS IN DATABASE ===")
comparisons <- dbGetQuery(con, "
  SELECT dataset_id, comparison_name, COUNT(*) as n_genes 
  FROM differential_expression 
  GROUP BY dataset_id, comparison_name
")
print(comparisons)

message("\n=== DATA SUMMARY ===")
de_count <- dbGetQuery(con, "SELECT COUNT(*) as total_de_results FROM differential_expression")
expr_count <- dbGetQuery(con, "SELECT COUNT(*) as total_expression_values FROM normalized_expression")
print(paste("Total differential expression results:", de_count$total_de_results))
print(paste("Total expression measurements:", expr_count$total_expression_values))

message("\n=== SIGNIFICANT GENES PREVIEW ===")
sig_genes <- dbGetQuery(con, "SELECT * FROM significant_genes LIMIT 5")
print(sig_genes)

message("\n=== SAMPLE DATABASE QUERIES ===")
# Test query - genes across datasets
test_query <- dbGetQuery(con, "
  SELECT dataset_id, comparison_name, COUNT(*) as significant_genes
  FROM significant_genes 
  GROUP BY dataset_id, comparison_name
  ORDER BY significant_genes DESC
")
print(test_query)

dbDisconnect(con)

message("\n=== DATABASE READY FOR WEB VIEWER ===")
message("Your database file: rnaseq_analysis_results.db")
message("You can now open this file in any SQLite web viewer")
message("Database creation complete!")

# Function to query any gene across all datasets
query_gene_across_datasets <- function(gene_symbol, db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, paste0("
    SELECT dataset_id, comparison_name, gene_symbol, gene_id,
           log2FoldChange, padj, baseMean
    FROM differential_expression 
    WHERE gene_symbol = '", gene_symbol, "' OR gene_id = '", gene_symbol, "'
    ORDER BY padj ASC
  "))
  dbDisconnect(con)
  return(result)
}

# Function to find consistently regulated genes
find_consistent_genes <- function(min_datasets = 3, max_padj = 0.05, db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, paste0("
    SELECT gene_symbol, 
           COUNT(*) as n_datasets,
           AVG(log2FoldChange) as avg_lfc,
           MIN(padj) as best_padj,
           GROUP_CONCAT(dataset_id || ':' || comparison_name, '; ') as found_in
    FROM differential_expression 
    WHERE padj <= ", max_padj, " AND gene_symbol IS NOT NULL AND gene_symbol != ''
    GROUP BY gene_symbol
    HAVING COUNT(*) >= ", min_datasets, "
    ORDER BY n_datasets DESC, best_padj ASC
  "))
  dbDisconnect(con)
  return(result)
}

message("\n=== EXAMPLE ADVANCED QUERIES ===")
message("Try these functions:")
message("query_gene_across_datasets('BRCA1')")
message("find_consistent_genes(min_datasets = 3)")

#advanced queries functions - good to have


# Function to query any gene across all datasets
query_gene_across_datasets <- function(gene_symbol, db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, paste0("
    SELECT dataset_id, comparison_name, gene_symbol, gene_id,
           log2FoldChange, padj, baseMean
    FROM differential_expression 
    WHERE gene_symbol = '", gene_symbol, "' OR gene_id = '", gene_symbol, "'
    ORDER BY padj ASC
  "))
  dbDisconnect(con)
  return(result)
}

# Function to find consistently regulated genes
find_consistent_genes <- function(min_datasets = 3, max_padj = 0.05, db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, paste0("
    SELECT gene_symbol, 
           COUNT(*) as n_datasets,
           AVG(log2FoldChange) as avg_lfc,
           MIN(padj) as best_padj,
           GROUP_CONCAT(dataset_id || ':' || comparison_name, '; ') as found_in
    FROM differential_expression 
    WHERE padj <= ", max_padj, " AND gene_symbol IS NOT NULL AND gene_symbol != ''
    GROUP BY gene_symbol
    HAVING COUNT(*) >= ", min_datasets, "
    ORDER BY n_datasets DESC, best_padj ASC
  "))
  dbDisconnect(con)
  return(result)
}

message("\n=== EXAMPLE ADVANCED QUERIES ===")
message("Try these functions:")
message("query_gene_across_datasets('BRCA1')")
message("find_consistent_genes(min_datasets = 3)")

#advanced queries - notes
# Query 1: Dataset overview with GEO information
showcase_dataset_overview <- function(db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, "
    SELECT 
      d.dataset_name,
      d.geo_id,
      d.description,
      d.n_samples,
      d.n_genes,
      COUNT(de.gene_id) as total_tested_genes,
      COUNT(CASE WHEN de.padj < 0.05 THEN 1 END) as significant_genes,
      ROUND(COUNT(CASE WHEN de.padj < 0.05 THEN 1 END) * 100.0 / COUNT(de.gene_id), 2) as percent_significant
    FROM datasets d 
    LEFT JOIN differential_expression de ON d.dataset_id = de.dataset_id
    GROUP BY d.dataset_id
    ORDER BY d.dataset_id
  ")
  dbDisconnect(con)
  return(result)
}

# Query 2: Top upregulated and downregulated genes per dataset
showcase_top_genes_per_dataset <- function(top_n = 5, db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  
  # Top upregulated
  up_genes <- dbGetQuery(con, paste0("
    SELECT 
      dataset_id,
      comparison_name,
      gene_symbol,
      log2FoldChange,
      padj,
      'upregulated' as direction
    FROM differential_expression 
    WHERE padj < 0.05 AND log2FoldChange > 0
    GROUP BY dataset_id
    HAVING log2FoldChange = MAX(log2FoldChange)
    ORDER BY dataset_id
  "))
  
  # Top downregulated  
  down_genes <- dbGetQuery(con, paste0("
    SELECT 
      dataset_id,
      comparison_name,
      gene_symbol,
      log2FoldChange,
      padj,
      'downregulated' as direction
    FROM differential_expression 
    WHERE padj < 0.05 AND log2FoldChange < 0
    GROUP BY dataset_id
    HAVING log2FoldChange = MIN(log2FoldChange)
    ORDER BY dataset_id
  "))
  
  dbDisconnect(con)
  
  # Combine results
  result <- rbind(up_genes, down_genes)
  return(result[order(result$dataset_id, result$direction), ])
}

# Query 3: Cross-dataset gene expression correlation
showcase_expression_patterns <- function(genes = c("EGFR", "TP53", "MYC"), db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  
  gene_list <- paste0("'", genes, "'", collapse = ",")
  result <- dbGetQuery(con, paste0("
    SELECT 
      de.gene_symbol,
      de.dataset_id,
      d.geo_id,
      de.comparison_name,
      de.log2FoldChange,
      de.padj,
      CASE 
        WHEN de.padj < 0.001 THEN 'highly_significant'
        WHEN de.padj < 0.05 THEN 'significant'
        ELSE 'not_significant'
      END as significance_level
    FROM differential_expression de
    JOIN datasets d ON de.dataset_id = d.dataset_id
    WHERE de.gene_symbol IN (", gene_list, ")
    ORDER BY de.gene_symbol, de.dataset_id
  "))
  
  dbDisconnect(con)
  return(result)
}

# Query 4: Sample expression distribution analysis
showcase_expression_distribution <- function(db_path = "rnaseq_analysis_results.db") {
  con <- dbConnect(SQLite(), db_path)
  result <- dbGetQuery(con, "
    SELECT 
      ROUND(expression_value) as expression_level,
      COUNT(*) as frequency
    FROM normalized_expression 
    WHERE expression_value BETWEEN 0 AND 15
    GROUP BY ROUND(expression_value)
    ORDER BY expression_level
  ")
  dbDisconnect(con)
  return(result)
}

# Query 5: Pathway-relevant genes (focusing on cancer/drug resistance)
showcase_pathway_genes <- function(db_path = "rnaseq_analysis_results.db") {
  # Cancer and drug resistance related genes
  cancer_genes <- c("TP53", "BRCA1", "BRCA2", "EGFR", "KRAS", "PIK3CA", "APC", "MYC", 
                   "PTEN", "RB1", "CDKN2A", "MLH1", "BRAF", "ERBB2", "MDM2")
  
  con <- dbConnect(SQLite(), db_path)
  
  gene_list <- paste0("'", cancer_genes, "'", collapse = ",")
  result <- dbGetQuery(con, paste0("
    SELECT 
      de.gene_symbol,
      COUNT(*) as datasets_found_in,
      AVG(de.log2FoldChange) as avg_log2fc,
      MIN(de.padj) as best_pvalue,
      COUNT(CASE WHEN de.padj < 0.05 THEN 1 END) as significant_datasets
    FROM differential_expression de
    WHERE de.gene_symbol IN (", gene_list, ")
    GROUP BY de.gene_symbol
    HAVING COUNT(*) >= 2
    ORDER BY significant_datasets DESC, best_pvalue ASC
  "))
  
  dbDisconnect(con)
  return(result)
}

message("\n=== SHOWCASE QUERIES FOR DISSERTATION ===")
message("Use these functions to demonstrate your database capabilities:")
message("1. showcase_dataset_overview() - Complete dataset summary with GEO IDs")
message("2. showcase_top_genes_per_dataset() - Most regulated genes per study") 
message("3. showcase_expression_patterns(c('EGFR', 'TP53')) - Track specific genes across datasets")
message("4. showcase_expression_distribution() - Expression level distributions")
message("5. showcase_pathway_genes() - Cancer-relevant genes across all studies")
message("\nThese queries showcase cross-dataset analysis, statistical insights, and biological relevance!")

file.exists("rnaseq_analysis_results.db")