softwareDevelopment/qtl_seq
0
1library(shiny)2library(shinyjs)3library(shinythemes)4library(shinyBS)5library(plotly)6library(DT)7library(QTLseqr)8library(data.table)9library(dplyr)10library(tidyr)11library(vcfR)12library(ggplot2)13library(dbscan)14library(factoextra)15library(cluster)16library(viridis)17 18options(shiny.maxRequestSize = 100 * 1024^2)19 20# Modified qtlseq function21qtlseq <- function(input_csv = NULL,22 refAlleleFreq = 0.20,23 minTotalDepth = 100,24 maxTotalDepth = 400,25 minSampleDepth = 40,26 minGQ = 99,27 windowSize = NULL,28 popStruc = NULL,29 bulkSize = NULL,30 replications = 10000,31 intervals = c(95, 99),32 alpha = 0.01,33 save_plot = FALSE,34 save_csv = FALSE) {35 data <- read.csv(input_csv)36 output_table <- tempfile(fileext = ".table")37 write.table(data, file = output_table, sep = "\t", row.names = FALSE, quote = FALSE)38 39 col_data <- colnames(data[5:ncol(data)])40 unique_samples <- unique(sub("\\..*", "", col_data))41 42 HighBulk <- unique_samples[1]43 LowBulk <- unique_samples[2]44 45 df <- importFromGATK(file = output_table, highBulk = HighBulk, lowBulk = LowBulk)46 47 df_filt <- filterSNPs(48 SNPset = df,49 refAlleleFreq = refAlleleFreq,50 minTotalDepth = minTotalDepth,51 maxTotalDepth = maxTotalDepth,52 minSampleDepth = minSampleDepth,53 minGQ = minGQ54 )55 56 df_filt <- runGprimeAnalysis(57 SNPset = df_filt,58 windowSize = windowSize,59 outlierFilter = "deltaSNP"60 )61 62 df_filt1 <- runQTLseqAnalysis(63 SNPset = df_filt,64 windowSize = windowSize,65 popStruc = popStruc,66 bulkSize = bulkSize,67 replications = replications,68 intervals = intervals69 )70 71 p <- ggplot(df_filt1, aes(x = POS / 1e6)) +72 geom_point(aes(y = deltaSNP), color = "blue", size = 0.5, alpha = 0.6) +73 geom_line(aes(y = tricubeDeltaSNP), color = "red", size = 0.8) +74 geom_line(aes(y = CI_95), color = "orange", linetype = "solid", size = 0.6) +75 geom_line(aes(y = -CI_95), color = "orange", linetype = "solid", size = 0.6) +76 geom_line(aes(y = CI_99), color = "green", linetype = "solid", size = 0.6) +77 geom_line(aes(y = -CI_99), color = "green", linetype = "solid", size = 0.6) +78 geom_hline(yintercept = c(-1, 0, 1), linetype = "dashed", color = "skyblue", alpha = 0.7) +79 facet_wrap(~CHROM, ncol = 4, scales = "free_x") +80 labs(x = "Position (Mb)", y = expression(Delta *"SNP-index"),81 title = "QTLseq Analysis: Delta SNP-index across Chromosomes") +82 theme_minimal(base_size = 12) +83 theme(84 strip.text = element_text(face = "bold", size = 10),85 axis.title = element_text(face = "bold"),86 plot.title = element_text(face = "bold", hjust = 0.5)87 )88 89 if (save_plot) ggsave("deltaSNPIndex.png", plot = p, dpi = 600, width = 12, height = 10, units = "in")90 91 qtl_table <- getQTLTable(SNPset = df_filt1, alpha = alpha, export = save_csv, fileName = "BSA_QTLseq.csv")92 93 return(list(QTLseqResults = df_filt1, GprimeResults = df_filt, plot = p, qtl_table = qtl_table))94}95 96# UI97ui <- fluidPage(98 theme = shinytheme("cerulean"),99 useShinyjs(),100 101 tags$head(102 tags$style(HTML("103 body { background-color: #f4f4f9; font-family: Arial, sans-serif; }104 .well { background-color: #ffffff; border: 1px solid #ddd; border-radius: 8px; padding: 20px; box-shadow: 0 2px 4px rgba(0,0,0,0.1); }105 .btn-primary { background-color: #4CAF50; border-color: #4CAF50; }106 .btn-primary:hover { background-color: #45a049; }107 h3 { color: #333; font-weight: bold; }108 .infographic { 109 border: 2px solid #2196F3; 110 border-radius: 15px; 111 padding: 20px; 112 margin: 15px 0; 113 background: linear-gradient(135deg, #e3f2fd 0%, #bbdefb 100%);114 box-shadow: 0 4px 8px rgba(0,0,0,0.1);115 }116 .info-header { 117 color: #1976d2; 118 border-bottom: 2px solid #64b5f6; 119 padding-bottom: 10px;120 margin-bottom: 15px;121 }122 .ml-card {123 background: white;124 border-radius: 10px;125 padding: 15px;126 margin: 10px 0;127 box-shadow: 0 2px 4px rgba(0,0,0,0.1);128 border-left: 4px solid #4CAF50;129 }130 "))131 ),132 133 titlePanel("QTLseq Shiny App: Advanced BSA with Machine Learning Insights"),134 135 sidebarLayout(136 sidebarPanel(137 fileInput("input_csv", "Upload Input CSV", accept = ".csv"),138 numericInput("refAlleleFreq", "Reference Allele Frequency", value = 0.20, min = 0, max = 1, step = 0.01),139 bsTooltip("refAlleleFreq", "Fraction of reads matching reference allele", placement = "right"),140 numericInput("minTotalDepth", "Min Total Depth", value = 100),141 bsTooltip("minTotalDepth", "Minimum total read depth for SNPs", placement = "right"),142 numericInput("maxTotalDepth", "Max Total Depth", value = 400),143 numericInput("minSampleDepth", "Min Sample Depth", value = 40),144 numericInput("minGQ", "Min GQ", value = 99),145 numericInput("windowSize", "Window Size", value = 5e5),146 selectInput("popStruc", "Population Structure", choices = c("F2", "RIL", "BC1")),147 textInput("bulkSize", "Bulk Size (e.g., 20,20)", value = "20,20"),148 numericInput("replications", "Replications", value = 10000),149 textInput("intervals", "Intervals (e.g., 95,99)", value = "95,99"),150 numericInput("alpha", "Alpha", value = 0.01, min = 0, max = 1, step = 0.01),151 actionButton("runAnalysis", "Run Analysis", class = "btn-primary")152 ),153 154 mainPanel(155 tabsetPanel(156 tabPanel("QTLseq Results",157 plotlyOutput("deltaPlot"),158 downloadButton("downloadPlot", "Download Plot"),159 DTOutput("qtlTable"),160 downloadButton("downloadCSV", "Download QTL CSV")161 ),162 163 tabPanel("ML Insights",164 h3("Machine Learning Integrations for Enhanced BSA Analysis"),165 166 tabsetPanel(167 tabPanel("K-Means Clustering",168 div(class = "infographic",169 h4("🎯 Automated Region Clustering with K-Means", class = "info-header"),170 p("K-Means clustering automatically groups genomic regions based on their deltaSNP patterns, revealing natural groupings without manual threshold setting."),171 plotlyOutput("clusterPlot"),172 DTOutput("clusterSummary")173 )174 ),175 176 tabPanel("Outlier Detection",177 div(class = "infographic",178 h4("🔍 DBSCAN Outlier SNP Detection", class = "info-header"),179 p("DBSCAN identifies statistical outliers that may represent key causative variants, sequencing artifacts, or novel biological phenomena worth investigating."),180 plotlyOutput("outlierPlot"),181 DTOutput("outlierTable")182 )183 ),184 185 tabPanel("PCA Overview",186 div(class = "infographic",187 h4("📊 Genome-wide Pattern Analysis with PCA", class = "info-header"),188 p("Principal Component Analysis provides a big-picture view of genomic patterns, helping identify batch effects, population structure, or technical artifacts."),189 plotlyOutput("pcaPlot"),190 plotOutput("pcaVariancePlot")191 )192 )193 )194 ),195 196 tabPanel("Advanced Visualizations",197 h3("Enhanced Interactive Infographics"),198 199 div(class = "infographic",200 h4("📈 Interactive Chromosome Explorer", class = "info-header"),201 p("Explore each chromosome individually with enhanced zoom and hover capabilities for detailed QTL region investigation."),202 selectInput("chromSelect", "Select Chromosome:", choices = NULL),203 plotlyOutput("chromosomePlot")204 ),205 206 div(class = "infographic",207 h4("🌊 3D DeltaSNP Landscape", class = "info-header"),208 p("Visualize the deltaSNP landscape across chromosomes in 3D, revealing spatial patterns and relationships between genomic regions."),209 plotlyOutput("threeDPlot")210 )211 )212 )213 )214 )215)216 217# Server218server <- function(input, output, session) {219 220 results <- reactiveVal()221 222 observeEvent(input$runAnalysis, {223 req(input$input_csv)224 showNotification("Running QTLseq analysis...", type = "message")225 226 bulkSize_vec <- as.numeric(unlist(strsplit(input$bulkSize, ",")))227 intervals_vec <- as.numeric(unlist(strsplit(input$intervals, ",")))228 229 tryCatch({230 res <- qtlseq(231 input_csv = input$input_csv$datapath,232 refAlleleFreq = input$refAlleleFreq,233 minTotalDepth = input$minTotalDepth,234 maxTotalDepth = input$maxTotalDepth,235 minSampleDepth = input$minSampleDepth,236 minGQ = input$minGQ,237 windowSize = input$windowSize,238 popStruc = input$popStruc,239 bulkSize = bulkSize_vec,240 replications = input$replications,241 intervals = intervals_vec,242 alpha = input$alpha243 )244 245 results(res)246 showNotification("Analysis completed successfully!", type = "message")247 248 # Update chromosome selector249 chrom_choices <- unique(res$QTLseqResults$CHROM)250 updateSelectInput(session, "chromSelect", choices = chrom_choices)251 252 }, error = function(e) {253 showNotification(paste("Error:", e$message), type = "error")254 })255 })256 257 # QTLseq Plot258 output$deltaPlot <- renderPlotly({259 req(results())260 ggplotly(results()$plot, tooltip = c("x", "y", "CHROM")) %>%261 layout(height = 800)262 })263 264 # Download Plot265 output$downloadPlot <- downloadHandler(266 filename = "deltaSNPIndex.png",267 content = function(file) {268 ggsave(file, plot = results()$plot, dpi = 600, width = 12, height = 10, units = "in")269 }270 )271 272 # QTL Table273 output$qtlTable <- renderDT({274 req(results())275 datatable(results()$qtl_table, 276 options = list(pageLength = 10, scrollX = TRUE),277 caption = "Significant QTL Regions")278 })279 280 # Download CSV281 output$downloadCSV <- downloadHandler(282 filename = "BSA_QTLseq.csv",283 content = function(file) {284 write.csv(results()$qtl_table, file, row.names = FALSE)285 }286 )287 288 # K-Means Clustering289 output$clusterPlot <- renderPlotly({290 req(results())291 df <- results()$QTLseqResults292 293 # Prepare data for clustering294 cluster_data <- df %>% 295 select(POS, tricubeDeltaSNP) %>% 296 drop_na() %>%297 sample_n(min(10000, nrow(.))) # Sample for performance298 299 # Perform K-Means Clustering300 set.seed(123)301 kmeans_result <- kmeans(cluster_data, centers = 3)302 303 # Add cluster assignment304 df_sample <- df %>%305 filter(row_number() %in% as.numeric(rownames(cluster_data))) %>%306 mutate(Cluster = as.factor(kmeans_result$cluster))307 308 p <- ggplot(df_sample, aes(x = POS/1e6, y = tricubeDeltaSNP, color = Cluster,309 text = paste("CHROM:", CHROM, "<br>POS:", POS, 310 "<br>Cluster:", Cluster))) +311 geom_point(alpha = 0.7, size = 2) +312 scale_color_viridis_d(option = "plasma") +313 labs(x = "Position (Mb)", y = "Tricube Delta SNP-index",314 title = "K-Means Clustering of Genomic Regions",315 subtitle = "Colors represent automated grouping based on deltaSNP patterns") +316 theme_minimal() +317 theme(legend.position = "bottom")318 319 ggplotly(p, tooltip = "text") %>%320 layout(height = 500)321 })322 323 output$clusterSummary <- renderDT({324 req(results())325 df <- results()$QTLseqResults326 327 cluster_data <- df %>% select(POS, tricubeDeltaSNP) %>% drop_na()328 set.seed(123)329 kmeans_result <- kmeans(cluster_data, centers = 3)330 331 summary_df <- data.frame(332 Cluster = 1:3,333 Size = as.numeric(table(kmeans_result$cluster)),334 Center_Position = kmeans_result$centers[,1],335 Center_DeltaSNP = kmeans_result$centers[,2]336 )337 338 datatable(summary_df, options = list(dom = 't'),339 caption = "Cluster Summary Statistics")340 })341 342 # DBSCAN Outlier Detection343 output$outlierPlot <- renderPlotly({344 req(results())345 df <- results()$QTLseqResults346 347 # Prepare data348 dbscan_data <- df %>% 349 select(POS, deltaSNP) %>% 350 drop_na() %>%351 sample_n(min(5000, nrow(.))) # Sample for performance352 353 # Scale data354 scaled_data <- scale(dbscan_data)355 356 # Perform DBSCAN357 dbscan_result <- dbscan(scaled_data, eps = 0.5, minPts = 10)358 359 # Add outlier status360 df_sample <- df %>%361 filter(row_number() %in% as.numeric(rownames(dbscan_data))) %>%362 mutate(Outlier = ifelse(dbscan_result$cluster == 0, "Outlier", "Normal"))363 364 p <- ggplot(df_sample, aes(x = POS/1e6, y = deltaSNP, color = Outlier,365 text = paste("CHROM:", CHROM, "<br>POS:", POS,366 "<br>Status:", Outlier))) +367 geom_point(alpha = 0.7, size = 2) +368 scale_color_manual(values = c("Outlier" = "red", "Normal" = "blue")) +369 labs(x = "Position (Mb)", y = "Delta SNP-index",370 title = "DBSCAN Outlier Detection",371 subtitle = "Red points represent statistical outliers") +372 theme_minimal() +373 theme(legend.position = "bottom")374 375 ggplotly(p, tooltip = "text") %>%376 layout(height = 500)377 })378 379 output$outlierTable <- renderDT({380 req(results())381 df <- results()$QTLseqResults382 383 dbscan_data <- df %>% select(POS, deltaSNP) %>% drop_na()384 scaled_data <- scale(dbscan_data)385 dbscan_result <- dbscan(scaled_data, eps = 0.5, minPts = 10)386 387 outlier_df <- df %>%388 filter(row_number() %in% as.numeric(rownames(dbscan_data))) %>%389 mutate(Outlier = ifelse(dbscan_result$cluster == 0, "Yes", "No")) %>%390 filter(Outlier == "Yes") %>%391 select(CHROM, POS, deltaSNP, tricubeDeltaSNP, REF_FRQ)392 393 datatable(outlier_df, options = list(pageLength = 5, scrollX = TRUE),394 caption = "Detected Outlier SNPs")395 })396 397 # PCA Analysis398 output$pcaPlot <- renderPlotly({399 req(results())400 df <- results()$QTLseqResults401 402 # Sample data for PCA (for performance)403 pca_data <- df %>%404 group_by(CHROM) %>%405 sample_n(min(50, n())) %>%406 ungroup() %>%407 select(POS, REF_FRQ, deltaSNP) %>%408 drop_na()409 410 # Perform PCA411 pca_result <- prcomp(pca_data, scale. = TRUE)412 413 pca_df <- as.data.frame(pca_result$x)414 pca_df$CHROM <- df$CHROM[1:nrow(pca_df)]415 416 p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = CHROM,417 text = paste("PC1: ", round(PC1, 2), 418 "<br>PC2: ", round(PC2, 2)))) +419 geom_point(alpha = 0.7, size = 3) +420 scale_color_viridis_d(option = "magma") +421 labs(title = "PCA of Genomic Features",422 subtitle = "Each point represents a genomic region") +423 theme_minimal()424 425 ggplotly(p, tooltip = "text") %>%426 layout(height = 500)427 })428 429 output$pcaVariancePlot <- renderPlot({430 req(results())431 df <- results()$QTLseqResults432 433 pca_data <- df %>%434 group_by(CHROM) %>%435 sample_n(min(50, n())) %>%436 ungroup() %>%437 select(POS, REF_FRQ, deltaSNP) %>%438 drop_na()439 440 pca_result <- prcomp(pca_data, scale. = TRUE)441 442 fviz_eig(pca_result, addlabels = TRUE, 443 main = "PCA - Variance Explained by Principal Components",444 barfill = "#4CAF50", barcolor = "#4CAF50") +445 theme_minimal()446 })447 448 # Chromosome-specific plot449 output$chromosomePlot <- renderPlotly({450 req(results(), input$chromSelect)451 df <- results()$QTLseqResults452 453 chrom_df <- df %>% filter(CHROM == input$chromSelect)454 455 p <- ggplot(chrom_df, aes(x = POS/1e6, y = tricubeDeltaSNP,456 text = paste("Position:", POS, 457 "<br>DeltaSNP:", round(tricubeDeltaSNP, 3)))) +458 geom_line(color = "red", size = 1) +459 geom_ribbon(aes(ymin = -CI_99, ymax = CI_99), fill = "green", alpha = 0.1) +460 geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +461 labs(x = "Position (Mb)", y = "Tricube Delta SNP-index",462 title = paste("Chromosome", input$chromSelect, "Detailed View")) +463 theme_minimal()464 465 ggplotly(p, tooltip = "text") %>%466 layout(height = 400)467 })468 469 # 3D Plot470 output$threeDPlot <- renderPlotly({471 req(results())472 df <- results()$QTLseqResults473 474 # Sample data for 3D plot475 plot_data <- df %>%476 group_by(CHROM) %>%477 sample_n(min(100, n())) %>%478 ungroup()479 480 plot_ly(plot_data, x = ~POS/1e6, y = ~as.numeric(factor(CHROM)), z = ~tricubeDeltaSNP,481 type = 'scatter3d', mode = 'markers',482 color = ~tricubeDeltaSNP,483 colors = viridis::viridis(100),484 marker = list(size = 3),485 text = ~paste("CHROM:", CHROM, "<br>POS:", POS, 486 "<br>DeltaSNP:", round(tricubeDeltaSNP, 3))) %>%487 layout(scene = list(488 xaxis = list(title = 'Position (Mb)'),489 yaxis = list(title = 'Chromosome'),490 zaxis = list(title = 'Delta SNP-index')491 ),492 title = "3D Genomic Landscape of Delta SNP-index"493 )494 })495}496 497# Run the app498shinyApp(ui = ui, server = server)