FigureYa86SNPmotif
Ying Ge, Taojun Ye
2025-5-20
需求描述
突变会影响转录因子结合吗？作出判断，并同时画出motif logo和SNP的图。
#Requirement description
Does mutation affect transcription factor binding? Make a judgment and draw the motif logo and SNP diagram simultaneously.
应用场景
在基因组上同时展示突变位点和motif，为突变影响转录因子结合提供量化（pvalue、score）和可视化的证据。
在“小丫画图”公众号回复“SNP”，查看这个图在更多paper里的应用。
运行下面这行，查看motifbreakR的官方手册
#Application scenarios
Simultaneously displaying mutation sites and motifs on the genome provides quantitative (p-value, score) and visual evidence for the impact of mutations on transcription factor binding.
Reply to “SNP” on the “Xiaoya Painting” official account to check the application of this picture in more papers.
Run the following line to view the official manual of motionbreakR
# 打开motifbreakR包的示例文档浏览器（中文注释）
# Open the vignette browser for the motifbreakR package (English comment)
browseVignettes("motifbreakR")motifbreakR还可以跟其他工具结合使用，一系列结果图作为证据，帮你充实文章。看这篇：Variant Annotation Workshop with FunciVAR, StateHub and MotifBreakR
MotifbreakR can also be used in conjunction with other tools, with a series of result graphs as evidence to help enrich your article. Read this article: [Variant Annotation Workshop with FunciVAR, StateHub, and MotifBreakR]（ https://www.simoncoetzee.com/bioc2017.html )
环境设置
使用国内镜像安装包
#Environment settings
Use domestic image installation package
# 设置R包安装源为清华大学镜像（中文注释）
# Set R package installation source to Tsinghua University mirror (English comment)
options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# 设置Bioconductor包镜像为清华大学镜像（中文注释）
# Set Bioconductor package mirror to Tsinghua University mirror (English comment)
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
# 检查是否安装了BiocManager包，如果没有则安装（中文注释）
# Check if BiocManager package is installed, install if not (English comment)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 使用BiocManager安装motifbreakR包（中文注释）
# Install the motifbreakR package using BiocManager (English comment)
BiocManager::install("motifbreakR")
# 如果你只提供rs ID，就需要安装这个包（中文注释）
# If you only provide rs IDs, you need to install this package (English comment)
# SNP locations and alleles for Homo sapiens extracted from NCBI dbSNP Build 151. The source data files used for this package were created by NCBI between February 16-22, 2018, and contain SNPs mapped to reference genome GRCh38.p7
# 这个版本的SNP文件480M，其他版本更大，建议下载后本地安装，<http://bioconductor.org/packages/3.8/data/annotation/src/contrib/SNPlocs.Hsapiens.dbSNP142.GRCh37_0.99.5.tar.gz>
# This version of the SNP file is 480M, other versions are even larger. It is recommended to download and install locally. <http://bioconductor.org/packages/3.8/data/annotation/src/contrib/SNPlocs.Hsapiens.dbSNP142.GRCh37_0.99.5.tar.gz>
BiocManager::install("SNPlocs.Hsapiens.dbSNP142.GRCh37")
# 如果你提供bed或vcf，有下面这个包就够了（中文注释）
# If you provide bed or vcf files, this package is sufficient (English comment)
# Full genome sequences for Homo sapiens (Human) as provided by UCSC (hg19, Feb. 2009) and stored in Biostrings objects.
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")其他物种到http://www.bioconductor.org/packages/3.8/data/annotation查找相应的包的名字。
你的电脑可能需要安装GhostScript，参考https://github.com/Simon-Coetzee/motifBreakR里的Prepairing to install
加载包
Other species to< http://www.bioconductor.org/packages/3.8/data/annotation >Search for the name of the corresponding package.
Your computer may need to install GhostScript, refer to< https://github.com/Simon-Coetzee/motifBreakR >Preparing to install in
Library Package
# 加载motifbreakR包 - 用于分析SNP对转录因子结合位点的影响（中文注释）
# Load the motifbreakR package - for analyzing the impact of SNPs on transcription factor binding sites (English comment)
library(motifbreakR)## 载入需要的程辑包：grid## 载入需要的程辑包：MotifDb## 载入需要的程辑包：BiocGenerics##
## 载入程辑包：'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## 载入需要的程辑包：S4Vectors## 载入需要的程辑包：stats4##
## 载入程辑包：'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## 载入需要的程辑包：IRanges## 载入需要的程辑包：GenomicRanges## 载入需要的程辑包：GenomeInfoDb## 载入需要的程辑包：Biostrings## 载入需要的程辑包：XVector##
## 载入程辑包：'Biostrings'## The following object is masked from 'package:grid':
##
## pattern## The following object is masked from 'package:base':
##
## strsplit## See system.file("LICENSE", package="MotifDb") for use restrictions.## # 加载人类SNP位点信息包（基于GRCh37/hg19基因组和dbSNP142数据库）（中文注释）
# Load the human SNP locations package (based on GRCh37/hg19 genome and dbSNP142 database) (English comment)
# 包含约8400万个SNPs的位置、参考等位基因和常见变异等位基因信息（中文注释）
# Contains location, reference alleles, and common variant alleles for ~84 million SNPs (English comment)
library(SNPlocs.Hsapiens.dbSNP142.GRCh37)## 载入需要的程辑包：BSgenome## 载入需要的程辑包：BiocIO## 载入需要的程辑包：rtracklayer##
## 载入程辑包：'rtracklayer'## The following object is masked from 'package:BiocIO':
##
## FileForFormat# 加载人类基因组序列包（基于UCSC hg19版本）（中文注释）
# Load the human genome sequence package (based on UCSC hg19 version) (English comment)
library(BSgenome.Hsapiens.UCSC.hg19)
# 设置环境变量，使R显示英文错误信息（便于查找解决方案）（中文注释）
# Set environment variable to display R error messages in English (easier for troubleshooting) (English comment)
Sys.setenv(LANGUAGE = "en")
# 禁止R自动将字符向量转换为因子类型（中文注释）
# Prevent R from automatically converting character vectors to factors (English comment)
# 因子类型在数据分析中可能导致意外结果，通常建议关闭此自动转换（中文注释）
# Factors can cause unexpected results in data analysis, it's generally recommended to disable this automatic conversion (English comment)
options(stringsAsFactors = FALSE)输入数据
motifbreakR可以用SNP的rs ID，或bed文件，或vcf文件作为输入。
#Input data
MotifbreakR can use SNP’s rs ID, bed file, or vcf file as input.
提供SNP的rs ID
输入文件可以只提供SNP的rsID，例如rs1006140
##Provide rs ID for SNP
The input file can only provide the rsID of SNP, for example rs1006140
#SNPID <- read.table("easy_input_rs.txt")$V1
#SNPID
# 从文件读取SNP ID列表（当前已注释掉，使用下面的单个SNP ID）（中文注释）
# Read SNP ID list from file (currently commented out, using single SNP ID below) (English comment)
#SNPID <- read.table("easy_input_rs.txt")$V1
#SNPID
# 直接指定单个SNP ID进行分析（中文注释）
# Directly specify a single SNP ID for analysis (English comment)
SNPID <- "rs1927914"
# 从dbSNP数据库中获取指定SNP的详细信息（中文注释）
# Retrieve detailed information for the specified SNP from dbSNP database (English comment)
# 参数说明（Parameter explanation）:
# - rsid: 要查询的SNP ID (The SNP ID to query)
# - dbSNP: 使用的dbSNP数据库版本 (The version of dbSNP database to use)
# - search.genome: 参考基因组 (Reference genome)
SNPinfo <- snps.from.rsid(rsid = SNPID,
dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,
search.genome = BSgenome.Hsapiens.UCSC.hg19)
# 显示获取的SNP信息（中文注释）
# Display the retrieved SNP information (English comment)
SNPinfo## GRanges object with 1 range and 3 metadata columns:
## seqnames ranges strand | SNP_id REF
## <Rle> <IRanges> <Rle> | <character> <DNAStringSet>
## rs1927914 chr9 120464725 * | rs1927914 G
## ALT
## <DNAStringSet>
## rs1927914 A
## -------
## seqinfo: 1 sequence from hg19 genome提供突变位点的bed或vcf文件
easy_input.bed，整理自例文里的Table S1
##Provide bed or vcf files for mutation sites
easy_input.bed， Organize Table S1 from the example text
# 从BED格式文件中批量读取SNP信息（中文注释）
# Read SNP information in bulk from a BED format file (English comment)
# 参数说明（Parameter explanation）:
# - file: 输入的BED格式文件路径 (Path to the input BED format file)
# - search.genome: 参考基因组 (Reference genome)
# - format: 输入文件格式，支持"BED"或"VCF" (Input file format, supports "BED" or "VCF")
# SNPinfo <- snps.from.file(file = "easy_input.bed", # 输入BED格式文件 (Input BED format file)
# search.genome = BSgenome.Hsapiens.UCSC.hg19,
# format = "bed") # 也可以选择"vcf"格式 (Can also choose "vcf" format)判断SNP对motif的影响
#Assessing the impact of SNPs on motif
# 运行motifbreakR分析SNP对转录因子结合位点的影响（中文注释）
# Run motifbreakR to analyze the impact of SNPs on transcription factor binding sites (English comment)
# 参数说明（Parameter explanation）:
# - snpList: 包含SNP信息的对象 (Object containing SNP information)
# - filterp: 是否过滤p值 (Whether to filter by p-value)
# - pwmList: 使用的位置权重矩阵库 (Position Weight Matrix library to use)
# - threshold: 显著性阈值 (Significance threshold)
# - method: 评分方法 (Scoring method)
# - bkg: 背景碱基频率 (Background base frequencies)
# - BPPARAM: 并行计算参数 (Parallel computation parameters)
hocomoco <- query(MotifDb, "HOCOMOCO")
results <- motifbreakR(snpList = SNPinfo, filterp = TRUE,
pwmList = hocomoco,
threshold = 1e-4,
method = "ic",
bkg = c(A=0.25, C=0.25, G=0.25, T=0.25),
BPPARAM = BiocParallel::bpparam())
# 提取特定SNP(rs1927914)的分析结果（中文注释）
# Extract analysis results for a specific SNP (rs1927914) (English comment)
result1 <- results[names(results) %in% "rs1927914"]
# 去除结果中的重复行（中文注释）
# Remove duplicate rows from the results (English comment)
result1 <- unique(result1)
# 计算每个motif预测的p值（中文注释）
# Calculate p-values for each motif prediction (English comment)
result1pval <- calculatePvalue(result1)
# 显示计算得到的p值（中文注释）
# Display the calculated p-values (English comment)
result1pval## GRanges object with 2 ranges and 23 metadata columns:
## seqnames ranges strand | SNP_id REF
## <Rle> <IRanges> <Rle> | <character> <DNAStringSet>
## rs1927914 chr9 120464725 + | rs1927914 G
## rs1927914 chr9 120464725 - | rs1927914 G
## ALT varType motifPos geneSymbol dataSource
## <DNAStringSet> <character> <list> <character> <character>
## rs1927914 A SNV -10, 5 PO5F1 HOCOMOCOv11-core-A
## rs1927914 A SNV -5, 7 PO6F1 HOCOMOCOv10
## providerName providerId seqMatch
## <character> <character> <character>
## rs1927914 PO5F1_HUMAN.H11MO.0.A PO5F1_HUMAN.H11MO.0.A ctaggacttagcatGcataa..
## rs1927914 PO6F1_MOUSE.H10MO.D PO6F1_MOUSE.H10MO.D aggacttagcatAcataa..
## pctRef pctAlt scoreRef scoreAlt Refpvalue Altpvalue
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## rs1927914 0.848018 0.747277 9.30192 8.24358 5.56898e-05 9.32275e-04
## rs1927914 0.741396 0.854286 10.30735 11.83275 6.12885e-04 2.73734e-05
## snpPos alleleRef alleleAlt effect altPos alleleDiff
## <integer> <numeric> <numeric> <character> <integer> <numeric>
## rs1927914 <NA> NA NA strong 1 -1.05834
## rs1927914 <NA> NA NA strong 1 1.52540
## alleleEffectSize
## <numeric>
## rs1927914 -0.097108
## rs1927914 0.110523
## -------
## seqinfo: 1 sequence from hg19 genome# 计算突变型与野生型结合分数的差异（Altscore-Refscore）（中文注释）
# Calculate the difference in binding scores between mutant and wild-type (Altscore-Refscore) (English comment)
result1pval$score <- result1pval$scoreAlt-result1pval$scoreRef
# 将最终分析结果保存为CSV文件（中文注释）
# Save the final analysis results to a CSV file (English comment)
# write.csv(result1pval, "output.csv", quote = F)开始画图
#Start drawing
# 保存结果为PDF文件（当前已注释掉，直接显示图形）（中文注释）
# Save results as PDF file (currently commented out, displaying graph directly) (English comment)
#pdf("SNPmotif.pdf")
#grid.newpage()
# 可视化特定SNP(rs1927914)对motif结合的影响（中文注释）
# Visualize the impact of a specific SNP (rs1927914) on motif binding (English comment)
# 参数说明（Parameter explanation）:
# - results: 包含motif分析结果的对象 (Object containing motif analysis results)
# - rsid: 要可视化的SNP ID (SNP ID to visualize)
# - effect: 显示的效应强度阈值，"strong"显示强效应，"weak"显示所有效应 (Threshold for effect strength to display, "strong" shows strong effects, "weak" shows all effects)
plotMB(results = result1, rsid = "rs1927914",
effect = "strong") # 也可以选择"weak"显示所有预测结果 (Can also choose "weak" to show all predictions)## Loading required namespace: Cairo
# 关闭PDF设备（当前已注释掉）（中文注释）
# Close the PDF device (currently commented out) (English comment)
#dev.off()后期加工
根据output.csv中的最后一列score，向图中添加每个转录因子的Altscore-Refscore
#Post processing
Add the Altscore Refscore for each transcription factor to the graph based on the last column score in the output. csv
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8
## [5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8
## [7] LC_PAPER=zh_CN.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Asia/Shanghai
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3
## [2] SNPlocs.Hsapiens.dbSNP142.GRCh37_0.99.5
## [3] BSgenome_1.70.2
## [4] rtracklayer_1.62.0
## [5] BiocIO_1.12.0
## [6] motifbreakR_2.16.0
## [7] MotifDb_1.44.0
## [8] Biostrings_2.70.3
## [9] XVector_0.42.0
## [10] GenomicRanges_1.54.1
## [11] GenomeInfoDb_1.38.8
## [12] IRanges_2.36.0
## [13] S4Vectors_0.40.2
## [14] BiocGenerics_0.48.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.17.1
## [3] jsonlite_2.0.0 magrittr_2.0.3
## [5] GenomicFeatures_1.54.4 farver_2.1.2
## [7] rmarkdown_2.29 zlibbioc_1.48.2
## [9] vctrs_0.6.5 Cairo_1.6-2
## [11] memoise_2.0.1 Rsamtools_2.18.0
## [13] RCurl_1.98-1.17 base64enc_0.1-3
## [15] htmltools_0.5.8.1 S4Arrays_1.2.1
## [17] progress_1.2.3 curl_6.2.3
## [19] CNEr_1.38.0 SparseArray_1.2.4
## [21] Formula_1.2-5 sass_0.4.10
## [23] bslib_0.9.0 htmlwidgets_1.6.4
## [25] plyr_1.8.9 Gviz_1.46.1
## [27] cachem_1.1.0 GenomicAlignments_1.38.2
## [29] lifecycle_1.0.4 pkgconfig_2.0.3
## [31] Matrix_1.6-5 R6_2.6.1
## [33] fastmap_1.2.0 GenomeInfoDbData_1.2.11
## [35] MatrixGenerics_1.14.0 digest_0.6.37
## [37] colorspace_2.1-1 TFMPvalue_0.0.9
## [39] AnnotationDbi_1.64.1 Hmisc_5.2-3
## [41] RSQLite_2.3.11 seqLogo_1.68.0
## [43] filelock_1.0.3 httr_1.4.7
## [45] abind_1.4-8 compiler_4.3.0
## [47] bit64_4.6.0-1 htmlTable_2.4.3
## [49] backports_1.5.0 BiocParallel_1.36.0
## [51] DBI_1.2.3 R.utils_2.13.0
## [53] biomaRt_2.58.2 MASS_7.3-59
## [55] poweRlaw_1.0.0 rappdirs_0.3.3
## [57] DelayedArray_0.28.0 rjson_0.2.23
## [59] gtools_3.9.5 caTools_1.18.3
## [61] tools_4.3.0 splitstackshape_1.4.8
## [63] foreign_0.8-82 nnet_7.3-18
## [65] R.oo_1.27.1 glue_1.8.0
## [67] restfulr_0.0.15 checkmate_2.3.2
## [69] reshape2_1.4.4 cluster_2.1.4
## [71] ade4_1.7-23 TFBSTools_1.40.0
## [73] generics_0.1.4 gtable_0.3.6
## [75] tzdb_0.5.0 R.methodsS3_1.8.2
## [77] ensembldb_2.26.0 data.table_1.17.4
## [79] hms_1.1.3 xml2_1.3.8
## [81] pillar_1.10.2 stringr_1.5.1
## [83] dplyr_1.1.4 BiocFileCache_2.10.2
## [85] lattice_0.21-8 bit_4.6.0
## [87] deldir_2.0-4 annotate_1.80.0
## [89] biovizBase_1.50.0 DirichletMultinomial_1.44.0
## [91] tidyselect_1.2.1 GO.db_3.18.0
## [93] knitr_1.50 gridExtra_2.3
## [95] grImport2_0.3-3 ProtGenerics_1.34.0
## [97] SummarizedExperiment_1.32.0 xfun_0.52
## [99] Biobase_2.62.0 matrixStats_1.5.0
## [101] stringi_1.8.7 lazyeval_0.2.2
## [103] yaml_2.3.10 evaluate_1.0.3
## [105] codetools_0.2-19 interp_1.1-6
## [107] tibble_3.2.1 cli_3.6.5
## [109] rpart_4.1.19 xtable_1.8-4
## [111] jquerylib_0.1.4 dichromat_2.0-0.1
## [113] Rcpp_1.0.14 dbplyr_2.5.0
## [115] png_0.1-8 XML_3.99-0.18
## [117] parallel_4.3.0 readr_2.1.5
## [119] ggplot2_3.5.2 blob_1.2.4
## [121] prettyunits_1.2.0 latticeExtra_0.6-30
## [123] jpeg_0.1-11 AnnotationFilter_1.26.0
## [125] bitops_1.0-9 VariantAnnotation_1.48.1
## [127] scales_1.4.0 motifStack_1.46.0
## [129] crayon_1.5.3 rlang_1.1.6
## [131] KEGGREST_1.42.0