Analysis of the RNA binding protein motifs for RNA Seq

gene_x 0 like s 152 view s

Tags: pipeline

Analysis of the RNA binding protein motifs for RNA Seq

TODOs_NEXTWEEK: how to use RBPmap results to enrich the significant RBPs for RNA-seq (3' UTRs?) and miRNAs?

  1.     🧬 1. RBPmap (Web + Command-line)
  3.         Type: Web tool and downloadable Perl scripts
  5.         Purpose: Scans RNA sequences for known RBP motifs
  7.         Database: Uses curated RBP motif sets (e.g., from literature, CLIP data)
  9.         URL: https://rbpmap.technion.ac.il/
  11.     ✅ Good for: Biologists without strong R background
  12.     🚫 Not a package, but scriptable via CLI or HTML results

  1. filtering out RBP from the database

  2. improve the output of the 3utr_down_cleaned.fasta

  3. multiple testing

  4. perform this for miRNA.

http://xgenes.com/article/article-content/146/peak-and-motif-analyses-in-promoters/

http://xgenes.com/article/article-content/154/density-of-motif-plots-and-its-statistical-tests/

There are several alternative R packages and tools to perform motif enrichment analysis for RNA-binding proteins (RBPs), beyond PWMEnrich::motifEnrichment(). Here are the most notable ones:

| Tool / Package | Enrichment | Custom Motifs | CLI or R? | RNA-specific? |

| ------------------------ | ----------------- | --------------- | --------- | -------------- |

| PWMEnrich | ✅ | ✅ | R | ✅ |

| RBPmap | ✅ | ❌ (uses own db) | Web/CLI | ✅ | ----> try RBPmap_results + enrichments!

| Biostrings/TFBSTools | ❌ (only scanning) | ✅ | R | ❌ | #ATtRACT + Biostrings / TFBSTools

| rmap | ✅ (CLIP-based) | ❌ | R | ✅ |

| Homer | ✅ | ✅ | CLI | ⚠ RNA optional |

| MEME (AME, FIMO) | ✅ | ✅ | Web/CLI | ⚠ Generic |

  1. Using R to get 3UTR_sequences.fasta, 5UTR

    1. setwd("~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/motif_analysis")
    3. # === STEP 0: Load libraries ===
    4. if (!requireNamespace("biomaRt", quietly = TRUE)) install.packages("biomaRt")
    5. library(biomaRt)
    7. # === STEP 1: Load DEG files ===
    8. ~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/MKL-1_wt.EV_vs_parental-up.txt    #20086
    9. ~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/MKL-1_wt.EV_vs_parental-down.txt  #634
    10. ~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/WaGa_wt.EV_vs_parental-up.txt     #23832
    11. ~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/WaGa_wt.EV_vs_parental-down.txt   #375
    12. up_file <- "~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/MKL-1_wt.EV_vs_parental-up.txt"
    13. down_file <- "~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/results_2025_1/degenes/MKL-1_wt.EV_vs_parental-down.txt"
    15. up_degs <- read.table(up_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
    16. down_degs <- read.table(down_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
    18. # === STEP 2: Clean & extract Ensembl gene IDs ===
    19. up_ids <- unique(na.omit(up_degs[[1]]))
    20. down_ids <- unique(na.omit(down_degs[[1]]))
    22. # Optional: protein-coding filter
    23. # up_ids <- unique(up_degs$external_gene_name[up_degs$gene_biotype == "protein_coding"])
    24. # down_ids <- unique(down_degs$external_gene_name[down_degs$gene_biotype == "protein_coding"])
    26. # === STEP 3: Connect to Ensembl Asia mirror ===
    27. ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = "https://asia.ensembl.org")
    29. # === STEP 4: Robust, minimal-loss sequence fetch ===
    30. get_sequences <- function(ids, seqType, batch_size = 300, retry_depth = 3) {
    31. all_seqs <- list()
    32. failed_ids <- c()
    34. safe_batch_query <- function(batch_ids, depth = 1) {
    35.     if (length(batch_ids) == 0) return(NULL)
    37.     tryCatch({
    38.     seqs <- getSequence(id = batch_ids,
    39.                         type = "ensembl_gene_id",
    40.                         seqType = seqType,
    41.                         mart = ensembl)
    42.     seqs <- seqs[seqs[[seqType]] != "", ]
    43.     return(seqs)
    44.     }, error = function(e) {
    45.     if (depth < retry_depth && length(batch_ids) > 1) {
    46.         # Split batch and retry each half
    47.         split_batches <- split(batch_ids, ceiling(seq_along(batch_ids)/ceiling(length(batch_ids)/2)))
    48.         result <- lapply(split_batches, safe_batch_query, depth = depth + 1)
    49.         return(do.call(rbind, result))
    50.     } else {
    51.         message("❌ Final failure for ", length(batch_ids), " genes. Logging IDs.")
    52.         failed_ids <<- c(failed_ids, batch_ids)
    53.         return(NULL)
    54.     }
    55.     })
    56. }
    58. id_batches <- split(ids, ceiling(seq_along(ids) / batch_size))
    59. for (i in seq_along(id_batches)) {
    60.     batch <- id_batches[[i]]
    61.     cat(sprintf("    🔹 Batch %d/%d (%d genes)\n", i, length(id_batches), length(batch)))
    62.     result <- safe_batch_query(batch)
    63.     if (!is.null(result) && nrow(result) > 0) {
    64.     all_seqs[[length(all_seqs) + 1]] <- result
    65.     }
    66. }
    68. if (length(failed_ids) > 0) {
    69.     writeLines(failed_ids, paste0("failed_ids_", seqType, ".txt"))
    70.     message("⚠️ Failed IDs saved to failed_ids_", seqType, ".txt")
    71. }
    73. if (length(all_seqs) == 0) return(data.frame())
    74. return(do.call(rbind, all_seqs))
    75. }
    77. # === STEP 5: Write FASTA ===
    78. write_fasta <- function(sequences, names, file) {
    79. con <- file(file, "w")
    80. for (i in seq_along(sequences)) {
    81.     cat(">", names[i], "\n", sequences[i], "\n", file = con, sep = "")
    82. }
    83. close(con)
    84. }
    86. # === STEP 6: Background sampling ===
    87. cat("📦 Fetching background gene list...\n")
    88. all_deg_ids <- unique(c(up_ids, down_ids))
    89. all_genes <- getBM(attributes = c("ensembl_gene_id"), mart = ensembl)
    90. bg_ids <- setdiff(all_genes$ensembl_gene_id, all_deg_ids)
    91. bg_sample <- sample(bg_ids, 1000)
    93. # === STEP 7: Fetch sequences by group and type ===
    94. seq_types <- c("3utr", "5utr", "cds", "transcript")
    95. groups <- list(
    96. up = up_ids,
    97. down = down_ids,
    98. background = bg_sample
    99. )
    101. for (seq_type in seq_types) {
    102. cat(sprintf("\n⏳ Fetching %s sequences...\n", seq_type))
    103. for (group_name in names(groups)) {
    104.     cat(sprintf("  🔸 Group: %s\n", group_name))
    105.     ids <- groups[[group_name]]
    106.     seqs <- get_sequences(ids, seq_type)
    107.     if (nrow(seqs) > 0) {
    108.     out_file <- paste0(seq_type, "_", group_name, ".fasta")
    109.     write_fasta(seqs[[seq_type]], seqs$ensembl_gene_id, out_file)
    110.     cat("    ✅ Written to", out_file, "\n")
    111.     } else {
    112.     message("⚠️ No sequences found for ", group_name, " (", seq_type, ")")
    113.     }
    114. }
    115. }
    117. cat("\n✅ All FASTA files created for: 3′UTR, 5′UTR, CDS, and Transcript.\n")

  2. Using https://rbpmap.technion.ac.il/1747409685/results.html, how to predict p-value for a specific RBP?

    1. #Try using RBPmap
    2. ~/DATA/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/motif_analysis/run2/3utr_down.fasta

  3. Significant RBP enrichments

    1. install.packages("BiocManager")
    2. BiocManager::install("BSgenome")
    3. BiocManager::install("Biostrings")
    4. BiocManager::install("motifmatchr")
    5. #BiocManager::install("MotifDb")
    6. #BiocManager::install("motifStack")
    7. BiocManager::install("TFBSTools")
    8. #install.packages("pheatmap")
    10. library(Matrix)
    11. library(pheatmap)
    12. library(Biostrings)
    13. library(motifmatchr)
    14. library(GenomicRanges)
    15. #library(MotifDb)
    16. #library(motifStack)
    17. library(TFBSTools)
    18. library(SummarizedExperiment)
    19. library(ggplot2)
    20. library(pheatmap)
    22. library(PWMEnrich)
    23. library(PWMEnrich.Hsapiens.background)
    25. # ------------- Methods -----------
    27. clean_fasta_file <- function(infile, outfile) {
    28. lines <- readLines(infile)
    30. # Initialize variables to store cleaned lines
    31. cleaned_lines <- list()
    32. add_line <- FALSE  # To track when we need to add a sequence
    34. for (line in lines) {
    35.     # If the line is a header and contains "Sequence unavailable", skip it
    36.     if (startsWith(line, ">") && grepl("Sequence unavailable", line)) {
    37.     add_line <- FALSE  # Don't add this header or its sequence
    38.     } else if (startsWith(line, ">")) {
    39.     # If it's a new header, start adding the following sequence
    40.     add_line <- TRUE
    41.     cleaned_lines <- c(cleaned_lines, line)
    42.     } else if (add_line) {
    43.     # Add the sequence line if it follows a valid header
    44.     cleaned_lines <- c(cleaned_lines, line)
    45.     }
    46. }
    48. # Write cleaned lines to the output file
    49. writeLines(cleaned_lines, outfile)
    50. }
    52. # Function to clean sequences by replacing invalid characters with 'N'
    53. clean_invalid_sequences <- function(dna_string_set) {
    54. valid_bases <- c("A", "T", "C", "G", "N")  # Define valid DNA bases
    55. clean_sequences <- lapply(dna_string_set, function(seq) {
    56.     # Remove invalid characters, replace them with 'N'
    57.     seq <- gsub("[^ATCGN]", "N", as.character(seq))  # Replace non-ATCGN characters with 'N'
    58.     DNAString(seq)  # Return as a valid DNAString object
    59. })
    60. return(DNAStringSet(clean_sequences))  # Return cleaned sequences as a DNAStringSet
    61. }
    63. # ============================================
    64. # 2. Read and clean 3′ UTR sequences
    65. # ============================================
    66. # Clean your FASTA file
    67. clean_fasta_file("run2/3utr_up.fasta", "run2/3utr_up_cleaned.fasta")
    68. clean_fasta_file("run2/3utr_down.fasta", "run2/3utr_down_cleaned.fasta")
    69. clean_fasta_file("run2/3utr_background.fasta", "run2/3utr_background_cleaned.fasta")
    70. clean_fasta_file("run2/5utr_up.fasta", "run2/5utr_up_cleaned.fasta")
    71. clean_fasta_file("run2/5utr_down.fasta", "run2/5utr_down_cleaned.fasta")
    72. clean_fasta_file("run2/5utr_background.fasta", "run2/5utr_background_cleaned.fasta")
    74. # Read in the cleaned FASTA files
    75. utr_3_up_raw <- readDNAStringSet("run2/3utr_up_cleaned.fasta")
    76. utr_3_down_raw <- readDNAStringSet("run2/3utr_down_cleaned.fasta")
    77. utr_3_background_raw <- readDNAStringSet("run2/3utr_background_cleaned.fasta")
    78. utr_5_up_raw <- readDNAStringSet("run2/5utr_up_cleaned.fasta")
    79. utr_5_down_raw <- readDNAStringSet("run2/5utr_down_cleaned.fasta")
    80. utr_5_background_raw <- readDNAStringSet("run2/5utr_background_cleaned.fasta")
    82. # Clean sequences
    83. utr_3_up <- clean_invalid_sequences(utr_3_up_raw)
    84. utr_3_down <- clean_invalid_sequences(utr_3_down_raw)
    85. utr_3_background <- clean_invalid_sequences(utr_3_background_raw)
    86. utr_5_up <- clean_invalid_sequences(utr_5_up_raw)
    87. utr_5_down <- clean_invalid_sequences(utr_5_down_raw)
    88. utr_5_background <- clean_invalid_sequences(utr_5_background_raw)
    90. # Filter out sequences shorter than the shortest PWM (e.g. 6 bp)
    91. # and any sequences containing ambiguous bases (non‑ACGT)
    92. min_pwm_length <- 6
    93. utr_3_down_clean <- utr_3_down[
    94. width(utr_3_down) >= min_pwm_length &
    95. !grepl("[^ACGT]", as.character(utr_3_down))
    96. ]
    97. #cat("Kept", length(utr_down_clean), "sequences after cleaning.\n")
    98. utr_3_background_clean <- utr_3_background[
    99. width(utr_3_background) >= min_pwm_length &
    100. !grepl("[^ACGT]", as.character(utr_3_background))
    101. ]
    103. # Install packages (if not already)
    104. if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
    105. BiocManager::install(c("PWMEnrich", "Biostrings", "PWMEnrich.Hsapiens.background", "MotifDb"))
    106. BiocManager::install("PWMEnrich.Hsapiens.background")
    108. # Load libraries
    109. library(PWMEnrich)
    110. library(Biostrings)
    111. library(PWMEnrich.Hsapiens.background)
    112. library(MotifDb)
    114. # ============================================
    115. # 3. Load the hg19 GC‑aware background model
    116. # ============================================
    117. data("PWMLogn.hg19.MotifDb.Hsap",
    118.     package = "PWMEnrich.Hsapiens.background")
    119. bg_model <- PWMLogn.hg19.MotifDb.Hsap
    121. # ============================================
    122. # 4. Define your list of RBP gene symbols
    123. # ============================================
    125. # Specific RBP gene symbols
    126. rbp_genes <- c("A1CF", "A2BP1", "AC004381.6", "AC008073.5", "ACO1", "AKAP1", "AKAP17A", "AL662825.2", "AL662825.3", "AL662849.4", "AL844853.7", "ALKBH8", "ANKHD1", "ANKRD17", "APLF", "ASCC1", "BICC1", "BOLL", "BRAP", "BRUNOL4", "BRUNOL5", "BRUNOL6", "BX000357.2", "BX005143.2", "BX119957.7", "BX248507.1", "BX248518.2", "BX511012.1", "BX908728.4", "BX927220.1", "C14orf156", "C14orf21", "CARHSP1", "CELF1", "CELF2", "CELF3", "CIRBP", "CNBP", "CNOT4", "CPEB1", "CPEB2", "CPEB3", "CPEB4", "CPSF4", "CPSF4L", "CPSF6", "CPSF7", "CR388219.1", "CR388372.1", "CR753328.1", "CR759778.5", "CR759782.5", "CR847863.1", "CSDA", "CSDC2", "CSDE1", "CSTF2", "CSTF2T", "DAZ1", "DAZ2", "DAZ3", "DAZ4", "DAZAP1", "DAZL", "DDX41", "DDX43", "DDX53", "DHX57", "DHX8", "DNAJC17", "DPPA5", "DUS3L", "EIF2S1", "EIF3B", "EIF3G", "EIF4B", "EIF4H", "ELAVL1", "ELAVL2", "ELAVL3", "ELAVL4", "ENOX1", "ENOX2", "ENSG00000180771", "ENSG00000183403", "ENSG00000185728", "ENSG00000213250", "ENSG00000215042", "ENSG00000215492", "ENSG00000231942", "ENSG00000248163", "ENSG00000249536", "ENSG00000249644", "ENSG00000250177", "ERAL1", "ESRP1", "ESRP2", "EWSR1", "FMR1", "FUBP1", "FUBP3", "FUS", "Fusip1", "FXR1", "FXR2", "G3BP1", "G3BP2", "GRSF1", "HDLBP", "HELZ", "HNRNPA0", "HNRNPA1", "HNRNPA1L2", "HNRNPA2B1", "HNRNPA3", "HNRNPAB", "HNRNPC", "HNRNPCL1", "HNRNPD", "HNRNPF", "HNRNPH1", "HNRNPH2", "HNRNPH3", "hnRNPK", "HNRNPL", "hnRNPLL", "HNRNPM", "HNRNPR", "HNRPDL", "HTATSF1", "IGF2BP1", "IGF2BP2", "IGF2BP3", "KHDRBS1", "KHDRBS2", "KHDRBS3", "KHSRP", "KIAA0430", "KRR1", "LARP1", "LARP1B", "LARP4", "LARP4B", "LARP6", "LARP7", "LENG9", "LIN28A", "LIN28B", "MATR3", "MBNL1", "MBNL2", "MBNL3", "MDM2", "MEX3B", "MEX3C", "MEX3D", "MKI67IP", "MKRN1", "MKRN2", "MKRN3", "MRPS28", "MSI1", "MSI2", "MTHFSD", "MYEF2", "NCBP2", "NCBP2L", "NCL", "NEIL3", "NOL8", "NONO", "NOVA1", "NOVA2", "NPLOC4", "NUP153", "NUPL2", "PABPC1", "PABPC1L", "PABPC1L2A", "PABPC1L2B", "PABPC3", "PABPC4", "PABPC5", "PABPN1", "PABPN1L", "PAN3", "PARP12", "PCBP1", "PCBP2", "PCBP3", "PCBP4", "PDCD11", "PEG10", "PNMA3", "PNO1", "PNPT1", "POLDIP3", "POLR2G", "PPARGC1A", "PPARGC1B", "PPIE", "PPIL4", "PPP1R10", "PPRC1", "PRR3", "PSPC1", "PTBP1", "PTBP2", "PUF60", "PUM1", "PUM2", "QKI", "RALY", "RANBP2", "RAVER1", "RAVER2", "RBBP6", "RBCK1", "RBFOX2", "RBFOX3", "RBM10", "RBM11", "RBM12", "RBM12B", "RBM14-RBM4", "RBM15", "RBM15B", "RBM17", "RBM18", "RBM19", "RBM22", "RBM23", "RBM24", "RBM25", "RBM26", "RBM27", "RBM28", "RBM3", "RBM33", "RBM34", "RBM38", "RBM39", "RBM4", "RBM41", "RBM42", "RBM44", "RBM45", "RBM46", "RBM47", "RBM4B", "RBM5", "RBM6", "RBM7", "RBM8A", "RBMS1", "RBMS2", "RBMS3", "RBMX", "RBMX2", "RBMXL1", "RBMXL2", "RBMXL3", "RBMY1A1", "RBMY1B", "RBMY1D", "RBMY1E", "RBMY1F", "RBMY1J", "RBPMS", "RBPMS2", "RBP_Name", "RC3H1", "RC3H2", "RCAN2", "RDBP", "RDM1", "RNF113A", "RNF113B", "RNF31", "RNPC3", "RNPS1", "ROD1", "RP11-1286E23.4", "RRP7A", "RYBP", "SAFB", "SAFB2", "SAMD4A", "SAMD4B", "SART3", "SCAF4", "SCAF8", "SETD1A", "SETD1B", "SF1", "SF3B4", "SFPQ", "SHARPIN", "SLTM", "SNRNP35", "SNRNP70", "SNRPA", "SNRPB2", "SOLH", "SPEN", "SRBD1", "SREK1", "SRRT", "SRSF1", "SRSF11", "SRSF12", "SRSF2", "SRSF3", "SRSF4", "SRSF5", "SRSF6", "SRSF7", "SRSF9", "SSB", "STAR-PAP", "SUPT6H", "SYNCRIP", "TAB2", "TAB3", "TAF15", "TARDBP", "TDRD10", "TDRKH", "TEX13A", "THOC4", "TIA1", "TIAL1", "TMEM63A", "TMEM63B", "TNRC6A", "TNRC6B", "TNRC6C", "TOE1", "TRA2A", "TRA2B", "TRMT1", "TRMT2A", "TRNAU1AP", "TTC14", "U2AF1", "U2AF1L4", "U2AF2", "U2SURP", "UHMK1", "UNK", "UNKL", "UPF3B", "YAF2", "YB-1", "YBX2", "YTHDC1", "YTHDC2", "YTHDF1", "YTHDF2", "ZC3H10", "ZC3H13", "ZC3H15", "ZC3H18", "ZC3H3", "ZC3H4", "ZC3H6", "ZC3H7A", "ZC3H7B", "ZC3H8", "ZCCHC11", "ZCCHC13", "ZCCHC14", "ZCCHC17", "ZCCHC2", "ZCCHC3", "ZCCHC5", "ZCCHC6", "ZCCHC7", "ZCCHC8", "ZCCHC9", "ZCRB1", "ZFP36", "ZFP36L1", "ZFP36L2", "ZGPAT", "ZMAT5", "ZNF638", "ZRANB1", "ZRANB2", "ZRANB3", "ZRSR1", "ZRSR2")
    128. # ============================================
    129. # 5. Extract all PWMs from background & filter to RBPs
    130. # ============================================
    131. #all_pwms <- bg_model@pwms  #2287
    132. #rbp_pwms <- all_pwms[
    133. #  vapply(names(all_pwms),
    134. #         function(nm) any(nm %in% rbp_genes),
    135. #         logical(1))
    136. #]
    138. # 如果你想做模糊匹配（基因名作为子串出现），可以用：
    139. rbp_pwms <- all_pwms[grepl(paste(rbp_genes, collapse="|"), names(all_pwms), ignore.case=TRUE)]
    141. cat("Number of RBP PWMs selected:", length(rbp_pwms), "\n")
    142. cat("Example PWM names:\n", head(names(rbp_pwms), 10), "\n")
    144. # ============================================
    145. # 6. Create a custom background model with only RBP PWMs
    146. # ============================================
    147. custom_bg <- bg_model
    148. custom_bg@pwms <- rbp_pwms
    150. # ============================================
    151. # 7. Run motif enrichment
    152. # ============================================
    153. enrichment_result <- motifEnrichment(utr_3_down_clean, custom_bg)
    155. # ============================================
    156. # 8. Generate report & convert to data.frame
    157. # ============================================
    158. report_list <- groupReport(enrichment_result)
    159. report_df   <- as.data.frame(report_list)
    161. # ============================================
    162. # 9. Apply BH correction and filter significant hits
    163. # ============================================
    164. report_df$adj_pval <- p.adjust(report_df$p.value, method = "BH")
    165. significant_hits <- subset(report_df, adj_pval < 0.05)
    167. # ============================================
    168. # 10. View results
    169. # ============================================
    170. # Display motif name, raw p-value, adjusted p-value, and enrichment score
    171. colnames(significant_hits)[
    172. match(c("target","id","raw.score"), colnames(significant_hits))
    173. ] <- c("Gene", "Motif_ID", "EnrichmentScore")
    175. print(
    176. head(
    177.     significant_hits[, c("Gene","Motif_ID","EnrichmentScore","p.value","adj_pval")],
    178.     6
    179. )
    180. )
    182. p <- ggplot(significant_hits,
    183.             aes(x = reorder(Motif_ID, adj_pval),
    184.                 y = -log10(adj_pval))) +
    185. geom_col() +
    186. coord_flip() +
    187. labs(
    188.     x = "Motif ID",
    189.     y = expression(-log[10]("FDR‑adjusted p‑value")),
    190.     title = "Significant RBP Motif Enrichment"
    191. ) +
    192. theme_minimal()
    193. ggsave(
    194. filename = "RBP_motif_enrichment.png",  # or any path you like
    195. plot     = p,
    196. width    = 8,     # inches
    197. height   = 6,     # inches
    198. dpi      = 300    # resolution
    199. )

  4. Perform similar RBP-mapping by myself

    1. # Load a set of RBP motifs from a database (e.g., RBPDB) RBPDB, ENCODE, JASPAR, or even UCL's RBP motifs.
    3. # Set path to your PWM directory
    4. pfm_dir <- "/media/jhuang/Elements/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/motif_analysis/PFMDir"
    6. # List all .pfm files
    7. pfm_files <- list.files(pfm_dir, pattern = "\\.pfm$", full.names = TRUE)
    9. # Create an empty list to store PWMs
    10. pfm_list <- list()
    12. # Loop through each file
    13. for (file in pfm_files) {
    14. # Read table
    15. pfm <- tryCatch({
    16.     mat <- read.table(file, header = FALSE)
    17.     mat <- as.matrix(mat)
    19.     # Optionally add column names (assuming 4 columns for A, C, G, T)
    20.     if (ncol(mat) == 4) {
    21.     colnames(mat) <- c("A", "C", "G", "T")
    22.     }
    24.     mat
    25. }, error = function(e) {
    26.     message("Failed to read: ", file)
    27.     NULL
    28. })
    30. # Use the filename (without extension) as list name
    31. name <- tools::file_path_sans_ext(basename(file))
    32. pfm_list[[name]] <- pfm
    33. }
    35. # Check one example
    36. str(pfm_list[[1]])
    37. #saveRDS(pfm_list, file = "compiled_pfm_list.rds")
    39. #The error you're seeing (invalid PFM 'x': not an integer matrix) indicates that the PWM function expects the input matrix to represent raw position frequency matrices (PFMs), which are expected to be integer counts, but the matrices you have are likely log-odds matrices (floating point values), which represent the log-transformed scores instead of raw counts.
    41. # convert to PFMatrix
    42. pfm_named_list <- lapply(pfm_list, function(mat) {
    43. mat <- as.matrix(mat)
    44. if (!all(rownames(mat) %in% c("A", "C", "G", "T"))) {
    45.     rownames(mat) <- c("A", "C", "G", "T")
    46. }
    47. PFMatrix(ID = "motif", name = "motif", profileMatrix = mat)
    48. })
    50. # Create an empty list to store results
    51. match_results_list <- list()
    53. # Iterate through each PFMatrix object in the list
    54. for (i in seq_along(pfm_named_list)) {
    55. # Get the current PFMatrix object
    56. current_pfm <- pfm_named_list[[i]]
    58. # Run matchMotifs for the current motif
    59. match_result <- tryCatch({
    60.     #matchMotifs(pfm_obj, subject = utr_3_down, genome = NULL, p.cutoff = 1e-3)
    61.     matchMotifs(current_pfm, subject = utr_3_down, genome = NULL, p.cutoff = 1e-3)
    62. }, error = function(e) {
    63.     message("Error with motif ", names(pfm_named_list)[i], ": ", e$message)
    64.     return(NULL)  # Return NULL if error occurs
    65. })
    67. # Store the result if no error
    68. if (!is.null(match_result)) {
    69.     match_results_list[[names(pfm_named_list)[i]]] <- match_result
    70. }
    71. }
    73. # Check the results
    74. match_results_list
    76. # Extract motif match results for one motif
    77. motif_matches_1004_8676391 <- assay(match_results_list[["1004_8676391"]], "motifMatches")
    79. # Convert sparse matrix to dense matrix
    80. dense_matrix <- as.matrix(motif_matches_1004_8676391)
    82. # Optional: Add rownames for better visualization (if available)
    83. if (!is.null(names(utr_3_down))) {
    84. rownames(dense_matrix) <- names(utr_3_down)
    85. }
    87. # Convert logical matrix to numeric (1 for TRUE, 0 for FALSE)
    88. dense_matrix_num <- matrix(as.numeric(dense_matrix),
    89.                         nrow = nrow(dense_matrix),
    90.                         dimnames = dimnames(dense_matrix))
    92. # Now plot the heatmap
    93. png("motif_match_1004_8676391.png", width = 800, height = 1600)
    94. pheatmap(dense_matrix_num,
    95.         cluster_rows = FALSE,
    96.         cluster_cols = FALSE,
    97.         color = c("grey", "steelblue"),
    98.         legend_breaks = c(0, 1),
    99.         legend_labels = c("No match", "Match"),
    100.         main = "Motif Match: 1004_8676391")
    101. dev.off()
    103. #✅ Output
    104. #You will see a binary heatmap, where:
    105. #   * Blue cells = motif matched
    106. #   * White cells = no match
    107. #Let me know if you want to loop this over all motifs or save the plot to a file.

  5. Does hg19 vs hg38 Position Matter in Motif Enrichment?

    1. ✅ If you're doing motif enrichment analysis only (like PWMEnrich):
    3. Position DOES NOT matter.
    4. Here's why:
    5. Aspect  Explanation
    6. 🔍 PWMEnrich Only scans sequences for motif matches based on nucleotide content (e.g., enrichment of AU-rich motifs in downregulated 3′UTRs).
    7. 📏 Coordinates   Genomic coordinates (hg19 vs hg38) are not used in enrichment. Only the sequence itself matters.
    8. ⚖️ Statistical Test The enrichment is based on observed motif frequency vs background expectations (GC-aware or custom). No genome alignment is done.
    10. So:
    12.     ✅ As long as your 3′ UTR sequences are correct and biologically relevant, using an hg19 background is safe — even if your gene models are from hg38.
    14. ❗️ When does genome version matter?
    16. If you're doing things like:
    17. Use case    Position-sensitive?
    18. 🔬 Variant effect on motifs (e.g., SNPs in motifs)   ✅ Yes
    19. 🧭 Mapping motifs to genome positions (e.g., ChIP-seq overlaps)  ✅ Yes
    20. 📦 Loading sequences from coordinates (not FASTA)    ✅ Yes
    21. 🧠 Expression + genomic track integration (e.g., motif + ATAC-seq)   ✅ Yes
    23. Then you must ensure genome versions match exactly.
    24. ✅ In Your Case: Motif Enrichment of Downregulated 3′UTRs
    26.     If you are:
    28.     Feeding real 3′UTR sequences (as DNAStringSet) for downregulated genes
    30.     Using PWMEnrich purely for motif enrichment
    32. Then:
    34.     ✅ hg19 vs hg38 does not matter
    35.     ❗ Only sequence content matters
    37. Let me know if you want to later map enriched motifs back to genome positions — that’s when we’d switch to position-aware tools like FIMO, MOODS, or rtracklayer.

  6. Custom RBP-database preparation (failed!)

    1. # Example: Create a PWM for 2 fake motifs (replace with your actual PWMs)
    2. #motif1 <- matrix(c(0.2, 0.3, 0.3, 0.2,
    3. #                   0.1, 0.4, 0.4, 0.1,
    4. #                   0.25, 0.25, 0.25, 0.25,
    5. #                   0.3, 0.2, 0.3, 0.2), nrow = 4, byrow = TRUE,
    6. #                 dimnames = list(c("A", "C", "G", "T"), NULL))
    7. #motif2 <- matrix(c(0.3, 0.2, 0.2, 0.3,
    8. #                   0.4, 0.1, 0.4, 0.1,
    9. #                   0.25, 0.25, 0.25, 0.25,
    10. #                   0.2, 0.3, 0.2, 0.3), nrow = 4, byrow = TRUE,
    11. #                 dimnames = list(c("A", "C", "G", "T"), NULL))
    12. #rbp_pwms <- list(RBP1 = motif1, RBP2 = motif2)
    14. pwm_dir <- "/media/jhuang/Elements/Data_Ute/Data_RNA-Seq_MKL-1+WaGa/motif_analysis/PWMDir"
    15. pwm_files <- list.files(pwm_dir, pattern = "\\.pwm$", full.names = TRUE)
    16. rbp_pwms <- list()
    17. for (file in pwm_files) {
    18. # Read table
    19. pwm <- tryCatch({
    20.     mat <- read.table(file, header = FALSE)
    21.     mat <- as.matrix(mat)
    23.     # Optionally add column names (assuming 4 columns for A, C, G, T)
    24.     if (ncol(mat) == 4) {
    25.     colnames(mat) <- c("A", "C", "G", "T")
    26.     }
    28.     mat
    29. }, error = function(e) {
    30.     message("Failed to read: ", file)
    31.     NULL
    32. })
    34. # Use the filename (without extension) as list name
    35. name <- tools::file_path_sans_ext(basename(file))
    36. rbp_pwms[[name]] <- pwm
    37. }
    39. # -- filtering (ERROR) --
    40. pwm_list <- motifs(bg_model)
    41. # Example: filter motifs containing 'ELAVL' or 'HNRNP'
    42. #rbp_pwms <- pwm_list[grep("ELAVL|HNRNP|IGF2BP|RBFOX|FMR1|PUM", names(pwm_list))]
    43. rbp_enrichment <- motifEnrichment(utr_3_down, bg_model, pwmList = rbp_pwms)
    44. rbp_report <- groupReport(rbp_enrichment)
    46. # -- from MotifDb --
    47. library(MotifDb)
    48. human_motifs <- query(MotifDb, "Hsapiens")
    49. rbp_motifs <- query(human_motifs, "RBP")
    51. # ---- try Build your background model directly from these PWMs ----
    53. # 1. List all PWM files
    54. files_all <- list.files("~/Downloads/cisbp-rna/pwms_all_motifs", pattern="\\.txt$", full.names=TRUE)
    56. # 2. Filter out empty files
    57. file_sizes <- file.info(files_all)$size
    58. files     <- files_all[file_sizes > 0]
    59. cat("Reading", length(files), "non‑empty PWM files out of", length(files_all), "\n")
    61. # 3. Define the reader (as before)
    62. readCisbpPfm <- function(path) {
    63. df <- read.table(path, header=TRUE, stringsAsFactors=FALSE)
    64. mat_rna <- as.matrix(df[, c("A","C","G","U")])
    65. mat_dna <- t(mat_rna)
    66. rownames(mat_dna) <- c("A","C","G","T")
    67. sweep(mat_dna, 2, colSums(mat_dna), FUN="/")
    68. }
    70. # 4. Read only the non‐empty files
    71. cisbp_mats <- lapply(files, readCisbpPfm)
    72. names(cisbp_mats) <- tools::file_path_sans_ext(basename(files))
    74. # 5. (Optional) Check the first few
    75. head(names(cisbp_mats))
    76. length(cisbp_mats)  #193
    78. # 6. Build your background model directly from these PWMs:
    79. bg_utrs <- utr_3_background_clean    #readDNAStringSet("run2/3utr_background.fasta")
    80. #bg_utrs <- bg_utrs[ width(bg_utrs)>=6 & !grepl("[^ACGT]", as.character(bg_utrs)) ]
    82. custom_bg <- makeBackground(
    83. seqs    = bg_utrs,
    84. motifs = cisbp_mats,
    85. type    = "logn",
    86. organism = "hg19",
    87. verbose  = TRUE
    88. )
    90. custom_bg <- makeBackground(
    91. sequences = bg_utrs,       # your DNAStringSet of background UTRs
    92. motifs    = cisbp_mats,    # your named list of PWM matrices
    93. type      = "logn",        # log‑normal GC‑aware model
    94. organism  = NULL,          # disable built‑in genomes
    95. verbose   = TRUE
    96. )
    98. custom_bg <- makeBackground(
    99. cisbp_mats,        # 1) motifs
    100. organism = "hg19",
    101. type    = "logn",  # 2) background model type
    102. bg.seq  = bg_utrs, # 3) override organism by supplying your own sequences
    103. quick   = FALSE    # (optional) fit to all sequences for best accuracy
    104. )
    106. # 1. Pick a library size to re‑scale your PWMs into PFMs
    107. library_size <- 1000   # you can also try 1000
    109. # 2. Convert each normalized matrix into integer counts
    110. cisbp_pfms <- lapply(cisbp_mats, function(pmat) {
    111. counts <- round(pmat * library_size)
    112. # ensure no column sums to zero (add pseudocount 1 if needed)
    113. zero_cols <- which(colSums(counts) == 0)
    114. if (length(zero_cols) > 0) {
    115.     counts[, zero_cols] <- 1
    116. }
    117. counts
    118. })
    120. # 3. Build a custom background (PFMs + your bg UTRs)
    121. custom_bg <- makeBackground(
    122. cisbp_pfms,      # 1) a list of integer PFMs
    123. type   = "logn", # 2) log‑normal GC‑aware model
    124. bg.seq = bg_utrs,# 3) your background DNAStringSet
    125. quick  = FALSE   # fit on all sequences
    126. )
    128. # … then continue with motifEnrichment() on your down‑regulated UTRs as before.
    130. # 5) Read & clean your down‑regulated 3′UTRs
    131. down_utrs <- readDNAStringSet("downregulated_utrs.fasta")
    132. down_utrs <- down_utrs[ width(down_utrs)>=6 & !grepl("[^ACGT]", as.character(down_utrs)) ]
    134. # 6) Run the enrichment
    135. res    <- motifEnrichment(down_utrs, custom_bg)
    136. report <- groupReport(res)
    137. df     <- as.data.frame(report)
    138. df$adj_pval <- p.adjust(df$p.value, method="BH")
    140. # 7) Subset significant
    141. sig <- subset(df, adj_pval < 0.05)
    142. print(sig[, c("motif.name","p.value","adj_pval","enrichmentScore")])

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. pheatmap vs heatmap.2
  6. Setup conda environments
  7. Should the inputs for GSVA be normalized or raw?
  8. File format for single channel analysis of Agilent microarray data with Limma?
  9. Kraken2 Installation and Usage Guide
  10. PiCRUST2 Pipeline for Functional Prediction and Pathway Analysis in Metagenomics

最新文章

最多评论文章

推荐相似文章

Workflow using QIIME2 for Data_Marius_16S_2025

Variant calling for Data_Pietschmann_229ECoronavirus_Mutations_2025 (via docker own_viral_ngs)

How to debug and construct the docker docker own_viral_ngs?

Visualization and Export of miRNA Expression Profiles Using Manhattan Plots in R

© 2023 XGenes.com Impressum