Functions_for_Microbiome/250730_function for blast result at main · KitHubb/Functions_for_Microbiome · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
replace_tax_with_blast <- function(phyloseq_input, blast_lineage_file, blast_result_file, output_path) {
  # 필요한 패키지 로드
  if (!requireNamespace("crayon", quietly = TRUE)) install.packages("crayon")
  library(crayon)
  library(dplyr)
  library(readr)
  library(tibble)
  library(phyloseq)

  # Input validation
  if (!inherits(phyloseq_input, "phyloseq")) stop("phyloseq_input must be a phyloseq object")
  if (!file.exists(blast_lineage_file)) stop("BLAST lineage file does not exist")
  if (!file.exists(blast_result_file)) stop("BLAST result file does not exist")

  cat(green$bold("[1] Loading BLAST and lineage data...\n"))

  # 1. Import BLAST results and lineage
  lineage_out <- readr::read_table(blast_lineage_file) %>% as.data.frame()
  lineage_out$species <- paste0(lineage_out$species, "_", lineage_out$X9)
  lineage_out <- lineage_out[-1, ]

  blast <- readr::read_delim(blast_result_file, col_names = FALSE)
  colnames(blast) <- c("qacc", "taxid", "sacc", "evalue", "bitscore", "qcovus", "pident", "sscinames", "length")

  # 2. Merge BLAST and lineage data
  blast_out <- merge(blast, unique(lineage_out), by = "taxid")

  cat(green$bold("[2] Filtering BLAST results by identity >99%, coverage >99%, evalue < 1e-10...\n"))

  # 3. Filter BLAST results
  blast_out2 <- blast_out %>%
    unique() %>%
    dplyr::select(-sacc, -taxid) %>%
    dplyr::group_by(qacc) %>%
    dplyr::arrange(desc(pident), desc(qcovus), desc(length), evalue) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    as.data.frame()

  blast_out3 <- blast_out2 %>%
    filter(pident > 99 & qcovus > 99 & evalue < 1e-10)

  cat(green$bold(paste0("    ✔ BLAST hits passing cutoff: ", nrow(blast_out3), " sequences\n")))

  # 4. Format BLAST results for tax_table
  blast_out4 <- blast_out3[, c("qacc", "kingdom", "phylum", "class", "order", "family", "genus", "species")] %>%
    column_to_rownames("qacc")
  colnames(blast_out4) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")

  # 5. Extract unidentified taxa from phyloseq
  cat(green$bold("[3] Extracting unidentified sequences from phyloseq...\n"))
  ps <- phyloseq_input
  tax <- tax_table(ps) %>% as.data.frame()
  tax_un <- tax[grepl("_k$|_p$|_c$|_o$|_f$|_g$", tax$Species), ]
  seq_un <- rownames(tax_un)
  cat(green$bold(paste0("    ✔ Number of unidentified taxa in original tax_table: ", nrow(tax_un), "\n")))

  # 6. Find intersection of unidentified sequences and BLAST results
  cat(green$bold("[4] Matching unidentified sequences with BLAST hits...\n"))
  seq_both <- intersect(rownames(blast_out4), seq_un)
  cat(green$bold(paste0("    ✔ Matched sequences for update: ", length(seq_both), "\n")))

  blast_trans <- blast_out4[seq_both, ]
  cat(green$bold("[5] Preview of updated taxonomy from BLAST:\n"))
  print(blast_trans)

  # 7. Update taxonomic assignments (excluding Kingdom level)
  tax2 <- tax
  tax2[seq_both, 2:7] <- blast_trans[seq_both, 2:7]

  # 8. Combine updated and non-updated taxa
  tax3 <- rbind(tax2[rownames(tax2) %in% seq_both, ], tax[!rownames(tax) %in% seq_both, ])

  # 9. Update phyloseq object
  cat(green$bold("[6] Updating phyloseq object with new taxonomy...\n"))
  tax_table(ps) <- as.matrix(tax3)

  # 10. Save modified phyloseq object
  cat(green$bold(paste0("[7] Saving updated phyloseq object to: ", output_path, "\n")))
  saveRDS(ps, output_path)

  cat(green$bold("[✔] All steps completed successfully!\n"))

  return(ps)
}