可视化
使用基因组浏览器来可视化区域的染色质景观。Integrative Genomic Viewer查看器提供了易于使用的web应用程序版本和本地桌面版本。UCSC基因组浏览器提供了最全面的补充基因组信息
以IGV为例
去官网下载带java的IGV版本

igv下载

首先load genome，点击Genomes----> load genomes from sever... ---->选择小鼠参考基因组最新版本mm39

load genome

load file

特定基因的峰图

可以在红框内输入感兴趣的基因名称，从而可视化特定基因的峰图。左边的框可以选择感兴趣的染色体。同时可以通过右键调节峰图的类型（如点图，线图等），颜色，峰图高度，设置数据范围等，从而绘制出文章中展示的图

igv image

在CUT&Tag peaks上的热图：在 CUT&Tag 实验识别的峰值上生成的热图，以可视化这些峰值区域的信号强度。SEACR（Sparse Enrichment Analysis for CUT&RUN）用于识别 CUT&Tag 数据中的峰值。SEACR 的输出文件中，第六列包含峰值区域的起始和结束位置，格式为 chr:start-end，表示该区域内信号最强的位置。为了在热图中对齐信号，我们需要计算每个峰值区域的中点。这有助于在生成热图时，使不同峰值的信号在同一位置对齐。如，给定一个峰值区域 chr1:1000-1500，中点的位置为 1000 + (1500 - 1000) / 2 = 1250。创建一个新的 bed 文件，其中第六列包含峰值区域的中点信息。这将用于后续的热图可视化。使用deeptools 的 computeMatrix 和 plotHeatmap 命令来生成和绘制热图。computeMatrix 命令用于计算矩阵数据，这些数据将用于热图可视化。plotHeatmap 命令用于绘制热图，显示信号强度在不同峰值区域的分布。

使用DiffBind进行差异peaks分析

setwd("/home/zhangyina/zbw_cut-tag/7.peakcalling_last")
getwd()
library(DiffBind)
packageVersion("DiffBind")

#首先构建一个metadata.csv文件，bamReads是去除重复、过滤好的bam文件，Peaks是calling peak步骤生成的，Peacaller是识别peak用的工具，ControlID，bamControl填空即可，如图所示
#读入文件。保证路径是字符型型数据，不是因子数据，否则没办法识别路径
#peak caller 指的是用来进行peak calling的工具。macs：表示使用MACS2 工具进行峰值调用，默认情况下会使用 .xls 文件 作为输入格式，但是会有注释信息，会报错，输入.broadpeak文件
metadata <- read.csv("metadata_cuttag.csv", header=TRUE,stringsAsFactors=FALSE)
str(metadata)
#broadPeak 文件的得分列位于第五列。
#还是要指定scoreCol=5,指定得分列
dba_obj <- dba(sampleSheet = metadata,scoreCol = 5)
dba_obj

#dba.count() 函数的主要作用是根据样本的 .bam 文件和peak文件，计算在每个峰值区域内有多少读取（reads）与峰值区域重叠。这一步是差异分析的基础。
#bUseSummarizeOverlaps=FALSE（默认）：DiffBind 使用 BEDTools 通过简单的区域重叠来计算读取计数。这是 DiffBind 默认的计算方法，它依赖于峰值区域与读取重叠的区域计算。
#TRUE：使用summarizeOverlaps 函数来计算读取计数，是一种更严谨的方法。能够处理更复杂的对齐情况（例如处理重复区域、偏移问题等）。这种方法一般更适合处理复杂的数据集。
dba_obj <- dba.count(dba_obj,bUseSummarizeOverlaps=TRUE)
View(dba_obj)

#attributes=DBA_FACTOR会自动使用Factor属性，而label=DBA_ID会自动用样本ID来标记点。
dba.plotPCA(dba_obj, attributes=DBA_CONDITION, label=DBA_ID)
dev.off()  # 关闭当前打开的图形设备（如果有）

# Establishing a contrast，用于为差异分析设置对比条件
#DBA_CONDITION、DBA_TISSUE、DBA_TREATMENT、和 DBA_FACTOR 这些常量用来根据样本的不同属性（metadata）来区分样本
#DBA_CONDITION识别dba_obj$samples 中的 Condition 列中的信息，对应的是样本的实验条件。DBA_TISSUE识别的是 dba_obj$samples 中的 Tissue 列。DBA_TREATMENT识别 dba_obj$samples 中的 Treatment 列， 识别的是样本的处理类型。DBA_FACTOR 对应的是实验中的因子，通常表示的是不同的调控因子，比如转录因子（Transcription Factor）或者其他实验中研究的因素。识别 dba_obj$samples 中的 Factor 列。
dba_obj <- dba.contrast(dba_obj, categories=DBA_CONDITION, minMembers=2)#每个组要求至少有两个样本，否则无统计意义
dba.show(dba_obj, bContrasts=T) #结果显示在DBA 对象 dba_obj 中，定义了一个对比 (contrast=1)，在 组 E（2 个样本） 与 组 C（2 个样本） 之间进行的

#dba.analyze用于在指定的对比条件下执行差异结合分析（简称DBA）。其主要作用是根据之前通过 dba.contrast() 设置好的对比条件，对不同组样本的DNA结合位点（peaks）进行统计学分析，找出在这些条件下具有显著差异的结合位点。
#支持多种差异结合分析方法,DBA_ALL_METHODS：表示同时使用多个方法（如DESeq2和edgeR），便于比较不同方法得到的结果。
dba_obj <- dba.analyze(dba_obj,method=DBA_ALL_METHODS)
#dba_obj_edger <-  dba.analyze(dba_obj, method=DBA_EDGER)
#dba_obj_deseq2 <-  dba.analyze(dba_obj, method=DBA_DESEQ2)

#总结结果
# 
dba.plotVenn(dba_obj,contrast=1,method=DBA_ALL_METHODS) #Venn 图将展示 E 组与 C 组之间的 peaks 重叠和独有情况。
dba.plotVenn(dba_obj, contrast=1, method=DBA_EDGER,th=1)
dba.plotVenn(dba_obj, contrast=1, method=DBA_DESEQ2,th=1)

#th=1 表示将 p 值或 FDR 阈值设置为 1，即 不进行任何显著性过滤，所有检测到的 peaks 都将包含在 dba.report() 的输出结果中。
#默认情况下，th 通常设定为 0.05。
dba.edgeR <- dba.report(dba_obj, method=DBA_EDGER, contrast = 1, th=1)
dba.deseq <- dba.report(dba_obj, method=DBA_DESEQ2, contrast = 1, th=1)
dba.deseq <- as.data.frame(dba.deseq)
dba.edgeR <- as.data.frame(dba.edgeR)
View(dba.edgeR)
View(dba.deseq)
dba.edgeR_DE <-  subset(dba.edgeR, abs(Fold) > 2 & p.value < 0.05)
dba.edgeR.up<- dba.edgeR_DE[which(dba.edgeR_DE$Fold >= 2 ),]      # 显著上升的peak
dba.edgeR.down<- dba.edgeR_DE[which(dba.edgeR_DE$Fold <= 2),]    # 显著下降的peak
View(dba.edgeR.up)
View(dba.edgeR.down)

write.csv(dba.edgeR.up,file = "edgeR_up.csv",row.names = T)
write.csv(dba.edgeR.down,file = "edgeR_down.csv",row.names = T)

metadata.csv

使用ChIPseeker对peaks进行注释和可视化

###########R###############
#使用ChIPseeker需要准备两个文件：一个就是要注释的peaks的文件，需满足BED格式。另一个就是注释参考文件，即需要一个包含注释信息的TxDb对象
library(ChIPseeker)
library(GenomicRanges)
library(clusterProfiler) #功能富集分析
library(org.Mm.eg.db) #基因注释数据库
library(TxDb.Mmusculus.UCSC.mm39.refGene)
#TxDb.Mmusculus.UCSC.mm39.refGene基于 UCSC mm39 小鼠基因组构建的 RefSeq 基因注释信息。具体来说，这个包包含了基因组中的转录本、外显子、内含子、启动子等结构的注释信息，特别是与 RefSeq 数据库中基因相关的注释。
txdb <- TxDb.Mmusculus.UCSC.mm39.refGene

#getPromoters 用于获取基因组中的启动子区域，因为许多转录因子会结合在这些区域。
#getTagMatrix 函数用于生成Tag Matrix，这是一个表示ChIP-seq信号在特定基因组区域（如启动子区域）中的覆盖情况的矩阵。Tag Matrix通常用于后续的可视化分析，比如绘制热图或信号分布曲线。
# 定义TSS上下游的距离
#promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
#tagMatrix <- getTagMatrix(GC_peak0.01, windows=promoter)

#注释所有区域的峰值，无需指定tssRegionTSS区域附近
#TxDb提供峰值注释所需的基因组参考信息。annotatePeak 将使用这个数据库来确定峰值在基因组中的位置，并将它们注释为可能的基因组区域（如启动子、外显子、内含子等）
#annoDb 参数指定了一个基因注释数据库org.Mm.eg.db,数据库包含了基因的各种注释信息，如Entrez Gene ID、基因符号（Symbol）、基因名称（Gene Name）等。它将 annotatePeak 注释到的基因信息进一步转换或关联到特定的基因ID、名称或其他注释信息中。TxDb主要指定了基因组结构注释信息

#Conc：这是总体的峰值浓度（concentration）。表示在所有样本中结合峰的总体信号强度。它是所有样本在该峰位置的归一化信号值的汇总，或者是取了均值。
#Conc_E：这是条件 E 中的峰值浓度。表示在 E 条件下，该峰的信号强度。这通常是在 E 组样本中检测到的该区域的归一化信号值。
#Conc_C：这是条件 C 中的峰值浓度。表示在 C 条件下，该峰的信号强度。通常是 C 组样本中该区域的归一化信号值。

#将前面的表格转换为GRanges 对象
#GRanges 对象是 R 中的内存数据结构，专门用于存储和操作基因组区域数据。它不仅包含和 BED 文件类似的信息（染色体、起始位点、终止位点、链信息等），还可以附加其他元数据，并且能够直接与 R 语言的基因组分析工具进行无缝交互。
edgeR_down <- GRanges(
  seqnames = dba.edgeR.down$seqnames,  # 染色体列
  ranges = IRanges(start = dba.edgeR.down$start, end = dba.edgeR.down$end),  # 起始和结束位置
  strand = dba.edgeR.down$strand,  # 链向信息
  scoreall=dba.edgeR.down$Conc,
  scoreC=dba.edgeR.down$Conc_C,
  scoreE=dba.edgeR.down$Conc_E,
  Fold = dba.edgeR.down$Fold,  # 添加 Fold 变化作为元数据
  p.value = dba.edgeR.down$p.value,  # 添加 p 值作为元数据
  FDR = dba.edgeR.down$FDR  # 添加 FDR 值作为元数据
)
edgeR_down

edgeR_up <- GRanges(
  seqnames = dba.edgeR.up$seqnames,  # 染色体列
  ranges = IRanges(start = dba.edgeR.up$start, end = dba.edgeR.up$end),  # 起始和结束位置
  strand = dba.edgeR.up$strand,  # 链向信息
  scoreall=dba.edgeR.up$Conc,
  scoreC=dba.edgeR.up$Conc_C,
  scoreE=dba.edgeR.up$Conc_E,
  Fold = dba.edgeR.up$Fold,  # 添加 Fold 变化作为元数据
  p.value = dba.edgeR.up$p.value,  # 添加 p 值作为元数据
  FDR = dba.edgeR.up$FDR  # 添加 FDR 值作为元数据
)
edgeR_up

peak0.01_edgeR_down <- annotatePeak(edgeR_down,TxDb = txdb,annoDb = "org.Mm.eg.db")
peak0.01_edgeR_up <- annotatePeak(edgeR_up,TxDb = txdb,annoDb = "org.Mm.eg.db")

peak0.01_edgeR_down <- as.data.frame(peak0.01_edgeR_down)
View(peak0.01_edgeR_down)
peak0.01_edgeR_up <- as.data.frame(peak0.01_edgeR_up)
View(peak0.01_edgeR_up)


write.csv(peak0.01_edgeR_down,file = "peak0.01_edgeR_down.csv",quote = F,row.names = T)
write.csv(peak0.01_edgeR_up,file = "peak0.01_edgeR_up.csv",quote = F,row.names = T)

Cut&Tag与RNA-seq结合分析
通常前面已经得到差异表达的基因（RNA-seq）以及差异修饰区域对应的基因（Cut&Tag）。分别将两种分析的差异上下调基因做overlap，得到4个交集。对这些交集基因分别进行GO注释或KEGG分析，以期找到下游机制。

RNA_up <- read.table("./Enew_DEG_up.xlsx")
RNA_up$SYMBOL <- row.names(RNA_up)
View(RNA_up)
RNA_down <- read.table("./Enew_DEG_down.xlsx")
View(RNA_down)
RNA_down$SYMBOL <- row.names(RNA_down)
cut_edgeRup_RNAdown <- merge(peak0.01_edgeR_up, RNA_down, by = "SYMBOL")
View(cut_edgeRup_RNAdown)
# 对cuttag_down和RNA_up做overlap
cut_edgeRdown_RNAup <- merge(peak0.01_edgeR_down, RNA_up, by = "SYMBOL")
# 查看重叠的结果
View(cut_edgeRdown_RNAup)

##绘制韦恩图##
library(VennDiagram)
# 提取基因集合
# 从 cut_edgeRdown_RNAup 和 cut_edgeRup_RNAdown 中提取 SYMBOL 列
cut_down_genes <- unique(peak0.01_edgeR_down$SYMBOL)
RNA_up_genes <- unique(RNA_up$SYMBOL)

cut_up_genes <- unique(peak0.01_edgeR_up$SYMBOL)
RNA_down_genes <- unique(RNA_down$SYMBOL)

# 绘制 cut_down 和 RNA_up 的 Venn 图
venn_cutdown_RNAup <- venn.diagram(
  x = list(
    cut_down = cut_down_genes,
    RNA_up = RNA_up_genes
  ),
  category.names = c("CUT&Tag Downregulated", "RNA-seq Upregulated"),
  filename = NULL,#filename = NULL，表示不会将图像保存到文件中
  output = TRUE, #指定是否输出图像对象。这个对象可以后续使用 grid.draw() 绘制出来
  main = "Overlap between CUT&Tag downregulated genes and RNA-seq upregulated genes",
  fill = c("#0072B5FF","#BC3C29FF"),
  alpha = 0.7,#透明度为 50%，可以让重叠区域更加清晰地显示。
  cex = 1.5,#集合标签和重叠区域数字的文本大小
  fontface = "bold",
  cat.fontface = "bold",
  cat.pos = c(0, 0) # 将两个标签移动到圆的底部
)

venn2_cutdown_RNAup <- venn.diagram(
  x = list(
    cut_down = cut_down_genes,
    RNA_up = RNA_up_genes
  ),
  category.names = c("CUT&Tag Downregulated", "RNA-seq Upregulated"),
  filename = NULL,#filename = NULL，表示不会将图像保存到文件中
  output = TRUE, #指定是否输出图像对象。这个对象可以后续使用 grid.draw() 绘制出来
  main = "Overlap between CUT&Tag downregulated genes and RNA-seq upregulated genes",
  fill = c("#0072B5FF","#BC3C29FF"),
  alpha = 0.7,#透明度为 50%，可以让重叠区域更加清晰地显示。
  cex = 1.5,#集合标签和重叠区域数字的文本大小
  fontface = "bold",
  cat.fontface = "bold",
  cat.pos = c(180, 180), # 将两个标签移动到圆的底部
  scaled = FALSE # 关闭自动缩放，使两个圆大小相同
)

# 显示 Venn 图
grid.newpage()
grid.draw(venn_cutdown_RNAup)
grid.draw(venn2_cutdown_RNAup)
ggsave("venn_cutdown_RNAup.pdf",plot=venn_cutdown_RNAup,width = 9,height = 8)
ggsave("venn2_cutdown_RNAup.pdf",plot=venn2_cutdown_RNAup,width = 9,height = 8)
#对于其他情况同理

####DAVID通路注释分析
library(stringr)
cuttag_down_RNAseq_up_GO <- read.csv("./cuttag_RNAseq_GO.csv",header = T,stringsAsFactors = FALSE)
View(cuttag_down_RNAseq_up_GO)
cuttag_down_RNAseq_up_GO$Term <- gsub("GO:[0-9]+~","",cuttag_down_RNAseq_up_GO$Term)
cuttag_down_RNAseq_up_GO <- arrange(cuttag_down_RNAseq_up_GO,cuttag_down_RNAseq_up_GO[,10])
cuttag_down_RNAseq_up_GO$Term <- factor(cuttag_down_RNAseq_up_GO$Term)
#数据框已经排序，并且Description已经转换为因子，但为了确保ggplot2使用的是正确的因子顺序，我们仍需显式地设置因子的排列顺序。排序数据框并将列转换为因子并不能保证因子的级别按照排序顺序排列。
cuttag_down_RNAseq_up_GO$Term <- factor(cuttag_down_RNAseq_up_GO$Term, 
                                levels = cuttag_down_RNAseq_up_GO$Term[order(cuttag_down_RNAseq_up_GO$Fold.Enrichment)])
cuttag_down_RNAseq_up_GO$Term <- str_wrap(cuttag_down_RNAseq_up_GO$Term, width =30)
#要再重新排序
cuttag_down_RNAseq_up_GO$Term <- factor(cuttag_down_RNAseq_up_GO$Term, 
                                        levels = cuttag_down_RNAseq_up_GO$Term[order(cuttag_down_RNAseq_up_GO$Fold.Enrichment)])
#气泡图可以展示4个变量，横坐标是富集倍数，富集倍数高意味着输入基因集在特定条目中的基因数相对背景基因集的比例更大，表明该条目在输入基因集中具有更高的显著性或重要性。纵坐标是通路名称，颜色根据显著性来定义，气泡大小是下调基因中属于该通路的基因数目
#太长的术语换行
cutdown_RNAup_GO <- ggplot(cuttag_down_RNAseq_up_GO,aes(x=Fold.Enrichment,y = Term,fill=-1*log10(PValue),size=Count))+
  geom_point(shape = 21, stroke = 0.5) +  # shape = 21 代表空心圆，并设置边框粗细
  #scale_size(range=c(2, 8))+
  scale_size_continuous(range = c(8, 18))+
  #颜色梯度：显著性由低到高为由蓝到红
  #scale_colour_gradient(low = "#0072B5FF",high = "#BC3C29FF")+
  scale_fill_gradient(low = "#0072B5FF",high = "#BC3C29FF")+
  # 设置图形主题为白底
  theme_bw()+
  ylab("GO Terms")+
  xlab("Fold Enrichment")+
  #指定了颜色轴的标签。使用expression 函数来创建一个数学表达式，-log[10](PValue)
  labs(fill=expression(-log[10](PValue)))+
  #设置图例标题和图例文本的大小为14
  theme(legend.title=element_text(size=18), legend.text = element_text(size=16))+
  #y轴标题的右侧边距为50个像素，x轴标题的顶部边距为20个像素
  theme(axis.title.y = element_text(size=20,margin = margin(r = 50)),axis.title.x = element_text(size=20,margin = margin(t = 20)))+
  #设置x轴文本的样式。face = "bold" 将文本设置为粗体。color = "black" 将文本颜色设置为黑色。angle = 0 将文本角度设置为0度（水平显示）。vjust = 1 垂直调整设置为1，通常指定文本在其分配空间中的垂直位置。
  theme(axis.text.y = element_text(size = 17, face = "bold", color = "black"))+
  theme(axis.text.x = element_text(size=17,face ="bold",color="black",angle=0,vjust=1))
cutdown_RNAup_GO
ggsave("cutdown_RNAup_GO.pdf",plot=cutdown_RNAup_GO,width = 13,height = 12) 

绘制实验组和对照组的Cut&Tag修饰的平均信号曲线

#########linux命令##########
nohup samtools merge C_merge.bam C_bowtie2_redup.sort.bam C2_bowtie2_redup.sort.bam &
nohup samtools merge E_merge.bam E_bowtie2_redup.sort.bam E2_bowtie2_redup.sort.bam &
#使用bamCoverage进行格式转换，将bam文件转换成bw文件
nohup samtools sort -@ 10 C_merge.bam -o C_merge_sort.bam &
nohup samtools sort -@ 10 E_merge.bam -o E_merge_sort.bam &
samtools index C_merge_sort.bam
samtools index E_merge_sort.bam
###binsize用默认的
nohup bamCoverage -b /home/zhangyina/zbw_cut-tag/03.analysis/sam2bam/C_merge_sort.bam --normalizeUsing RPKM --extendReads 300 -o ~/zbw_cut-tag/03.analysis/bam2bw2/C_merge_sort.RPKM.bw &
nohup bamCoverage -b /home/zhangyina/zbw_cut-tag/03.analysis/sam2bam/E_merge_sort.bam --normalizeUsing RPKM --extendReads 300 -o ~/zbw_cut-tag/03.analysis/bam2bw2/E_merge_sort.RPKM.bw &

#reference-point 模式：基于一个参考点来计算信号分布。由于所有区域长度相同，可以简单地选择基因区域的中心（或起始点/终止点）作为参考点。
#--referencePoint center：指定参考点为基因区域的中心点
computeMatrix reference-point --referencePoint center -b 1000 -a 1000 -R /home/zhangyina/zbw_cut-tag/IGV/E_down_FC2_p0.01.bed -S ~/zbw_cut-tag/03.analysis/bam2bw2/C_merge_sort.RPKM.bw --skipZeros -o /home/zhangyina/zbw_cut-tag/03.analysis/TSS/matrix_C_down_merge_gene_3.mat.gz -p 20
computeMatrix reference-point --referencePoint center -b 1000 -a 1000 -R /home/zhangyina/zbw_cut-tag/IGV/E_down_FC2_p0.01.bed -S ~/zbw_cut-tag/03.analysis/bam2bw2/E_merge_sort.RPKM.bw --skipZeros -o /home/zhangyina/zbw_cut-tag/03.analysis/TSS/matrix_E_down_merge_gene_3.mat.gz -p 20
gunzip -k ./matrix_C_down_merge_gene_3.mat.gz
gunzip -k ./matrix_E_down_merge_gene_3.mat.gz

#############################平均信号曲线#################################
library(tidyverse)
library(stringr)
library(reshape2)

# 输入, sample C 和 E 的 output_denisty 路径
data_path_C <- "./matrix_C_down_merge_gene_3.mat"
data_path_E <- "./matrix_E_down_merge_gene_3.mat"

# 读取合并后的信号矩阵
#read_delim用于读取分隔符格式（delimited）的文件（如以逗号、制表符等分隔的文件）。delim = "\t"：指定文件的分隔符为制表符（\t）。col_names = FALSE：指定不自动将第一行作为列名。skip = 1：跳过文件的第一行。na.omit() 函数用于移除 data_agg 中的任何缺失值（NA）所在的行
# 整合A B的信号值
data_C <- read_delim(data_path_C,
                     delim = "\t",
                     skip = 1,
                     col_names = F) %>%
  na.omit()
data_E <- read_delim(data_path_E,
                     delim = "\t",
                     skip = 1,
                     col_names = F) %>%
  na.omit()

# 保留区间id和信号值列
data_agg_C <- data_C[,-c(1:3,5,6)]
## 计算每个区间的信号平均值
data_agg_C_mean <- colMeans(data_agg_C[,-1]) %>% as.data.frame()
# 添加bin的顺序，用于画图。c(1:nrow(data_agg_A_mean)) 生成一个从 1 到 data_agg_A_mean 行数的序列，标记每个 bin 的位置。
data_agg_C_mean$bin_number <- c(1:nrow(data_agg_C_mean))
# 添加样本标签，用于画图
data_agg_C_mean$type <- "Control"
# 添加列名
colnames(data_agg_C_mean) <- c("signal", "bin_number","type")

############# 样本E
data_agg_E <- data_E[,-c(1:3,5,6)]
# 计算每个区间的信号平均值
data_agg_E_mean <- colMeans(data_agg_E[,-1]) %>%
  as.data.frame()
# 添加bin的顺序，用于画图
data_agg_E_mean$bin_number <- c(1:nrow(data_agg_E_mean))
# 添加样本标签，用于画图
data_agg_E_mean$type <- " cKO"
# 添加列名
colnames(data_agg_E_mean) <- c("signal", "bin_number","type")

# 合并 A and B信号值
data_agg_mean <-  rbind(data_agg_C_mean, data_agg_E_mean)


data_plot <- data_agg_mean
#将 data_plot 中的 type 列转化为因子，levels = unique(data_plot$type) 指定因子的级别顺序，以确保图例按照 type 的顺序显示
data_plot$type <- factor(data_plot$type, levels  = unique(data_plot$type))

y_label <- "signal (RPKM)"
title_label <- ""

#使用 round() 函数取 y 值的最小和最大值，并向上调整，以确保显示完整的信号范围，round() 函数用于对数值进行四舍五入取整
y_limit <- c(round(min(data_plot$signal)),
                          round(max(data_plot$signal)) + 1)
y_limit <- c(2, round(max(data_plot$signal))+0.5)
#y 轴的刻度。seq() 从最小值到最大值按步长 2 生成刻度，确保 y 轴标记间隔为 2。
y_axis <- seq(round(min(data_plot$signal)),
              round(max(data_plot$signal))+1,
              2)

x_label <- "Distance to Center (kb)"
legend_label <- ""
fill_value <- c("#0072B5FF","#BC3C29FF")
size <- 2
y <- data_plot$signal
x <- data_plot$bin_number
x_limit <- c(0,max(data_plot$bin_number))
x_axis <- seq(0, max(data_plot$bin_number),
              max(data_plot$bin_number) / 2)

#图例的位置。c(0.8, 0.9) 表示图例位于绘图区域的右上角。
#legend_direction：图例排列方向，"vertical" 表示垂直排
legend_position <- c(0.8,0.9)
legend_direction <- "vertical"


center <- ggplot(data = data_plot, aes(x = x, y = y, color = type)) +
  geom_line(size = 2) +
  #breaks设置 y 轴刻度为 y_axis。limits设置 y 轴范围。expand控制 y 轴边缘留白，c(0,0) 表示无边缘留白
  scale_y_continuous(breaks = y_axis,
                     limits = y_limit,
                     expand = c(0,0)) +
  #labels = c("-5","-2.5","Center","2.5","5")：指定 x 轴标签，分别表示 -5kb 到 5kb 范围。
  scale_x_continuous(breaks = x_axis,
                     labels = c("-1.0Kb","Center","1.0Kb")) +
  labs(title = title_label, caption = "", y = y_label, x = x_label) +
  theme_classic(base_family = "sans",
                base_line_size = 0.3) +
  scale_color_manual(values = fill_value,
                     name = legend_label) +
  theme(legend.title=element_blank(),
        legend.position = legend_position,
        legend.direction = legend_direction,
        legend.text = element_text(size = 14, colour = "black"),
        legend.key.size = unit(0.1, 'inch'),
        legend.key.width = unit(1, "inch"),
        legend.background = element_rect(fill = "transparent"),
        panel.border = element_rect(color = "black", fill = NA, size = 0.8),  # 添加黑色边框
        axis.text = element_text(size = 14),        # 增大坐标轴标签的大小
        axis.title = element_text(size = 16),        # 增大坐标轴标题的大
        )
center
ggsave("DE_regions_signal2.pdf",plot=center,width = 10,height = 7)

Average signal

对目的基因区域和基因组上其他基因区域的已知motif分布做统计检验，并且绘制已知的motif在目的基因区域上的分布情况。

###############linux命令#########
#提取出目的基因注释信息
grep -F -f select.txt GRCm38_genecode_m25.gtf > select_genes.gtf
#-v 表示反转匹配（invert match）。它会输出那些 不匹配 指定模式的行。
grep -v -F -f select.txt GRCm38_genecode_m25.gtf > other_genes.gtf

#从中提取出所需的列
awk '$3 == "gene"' select_genes.gtf | awk '{print $1, $4-1000, $5, $9,$10,$13,$14}' OFS="\t" > select_promoters.bed
#在生成 BED 文件时，确保起始坐标不小于 0。可以通过 awk 添加一个判断条件来处理负值
awk '$3 == "gene"' other_genes.gtf | awk '{print $1, ($4-1000 < 0 ? 0 : $4-1000), $5, $9, $10,$13,$14}' OFS="\t" > other_promoters.bed

#提取出fasta序列
bedtools getfasta -fi GRCm38.primary_assembly.genome.fa -bed selcet_promoters.bed -fo select_promoters.fa
bedtools getfasta -fi GRCm38.primary_assembly.genome.fa -bed other_promoters.bed -fo other_promoters.fa

#统计motif数量
seqkit locate -p "ATCGGGT" select_promoters.fa | awk 'NR>1 {print $1}' | sort | uniq -c > select_motif_counts.txt
seqkit locate -p "ATCGGGT" other_promoters.fa | awk 'NR>1 {print $1}' | sort | uniq -c > other_motif_counts.txt

#需要将count=0保留
seqkit seq -n other_promoters.fa > all_sequences.txt

#######################R###################
#进行统计检验
library(dplyr)
# 读取文件
other_data <- read.table("other_motif_positions_counts.txt", header = FALSE, stringsAsFactors = FALSE)
all_sequences <- read.table("all_sequences_sorted.txt", header = FALSE, stringsAsFactors = FALSE)

# 为数据命名列（假设other_motif_counts.txt有两列：基因区域和计数）
colnames(other_data) <- c("region", "Count")
colnames(all_sequences) <- c("region")

# 将all_sequences添加到数据框中，确保包括所有基因区域
all_merged_data <- all_sequences %>%
  left_join(other_data, by = "region") %>%
  mutate(Count = ifelse(is.na(Count), 0, Count))

# 将合并后的数据写回文件
#write.table(merged_data, "updated_other_motif_counts.txt", 
            #row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")

# 读取数据
select_motif <- read.table("select_motif_counts.txt", header = FALSE, stringsAsFactors = FALSE)


# 给数据列命名
colnames(select_motif) <- c("Count", "Location")


# 转换 Count 列为数值
select_motif$Count <- as.numeric(seclect_motif$Count)
all_merged_data$Count <- as.numeric(all_merged_data$Count)
# 描述性统计
summary(select_motif$Count)
summary(all_merged_data$Count)

# 标准差
sd_select <- sd(select_motif$Count)
sd_non_all <- sd(all_merged_data$Count)

cat("Select Motif - Std Dev:", sd_select, "\n")
cat("Other Motif - Std Dev:", sd_non_all, "\n")

#Mann-Whitney U 检验（非参数检验）。适合小样本且不依赖正态分布的情况。比较两组数据的中位数或秩分布差异。
wilcox_test_all <- wilcox.test(select_motif$Count, all_merged_data$Count, alternative = "two.sided")
print(wilcox_test_all)
#单尾检验，检验 select_motif$Count 是否显著大于 other_motif$Count
wilcox_test_greater_all <- wilcox.test(select_motif$Count, all_merged_data$Count, alternative = "greater")
print(wilcox_test_greater_all)

#反向，是否小于
wilcox_test_less_all <- wilcox.test(select_motif$Count, all_merged_data$Count, alternative = "less")
print(wilcox_test_less_all)

library(ggplot2)

# 对数据进行对数变换，避免极端值影响
select_motif_log <- log1p(select_motif$Count)  # log1p是对数变换，log(1+x)
all_merged_data_log <- log1p(all_merged_data$Count)

# 创建数据框，合并两组数据
data_combined <- data.frame(
  Count = c(select_motif_log, all_merged_data_log),
  Group = c(rep("select Motif", length(select_motif_log)),
            rep("Other Motif", length(all_merged_data_log)))
)

# 绘制箱线图
ggplot(data_combined, aes(x = Group, y = Count, fill = Group)) +
  geom_boxplot(outlier.shape = 16, outlier.colour = "red", outlier.size = 3) +
  labs(title = "Boxplot of Select Motif vs Other Motif Counts (Log Scale)",
       x = "Group",
       y = "Log-Transformed Count") +
  theme_minimal() +
  scale_fill_manual(values = c("skyblue", "orange")) +
  theme(text = element_text(size = 12))
# 计算 Wilcoxon 检验 p 值（您可以直接使用检验结果中的 p-value）
p_value <- 0.001876
# 创建一个数据框，用于标注显著性
significance_label <- ifelse(p_value < 0.05, "*", "ns")  # "ns" 表示不显著
# 小提琴
violin_motif <-ggplot(data_combined, aes(x = Group, y = Count, fill = Group)) +
  geom_violin(trim = TRUE, alpha = 0.5) +  # 仅小提琴图
  labs(title = "Violin Plot of Select Motif vs Other Motif Counts (Log Scale)",
       x = "Group",
       y = "Log-Transformed Count") +
  scale_fill_manual(values = c("#BC3C29FF","#0072B5FF")) +
  theme_minimal() +
  theme(  # 覆盖字体大小
    plot.title = element_text(size = 20, face = "bold"), 
        axis.title.x = element_text(size = 20, face = "bold"),      # x 轴标题字体大小
        axis.title.y = element_text(size = 20, face = "bold"),      # y 轴标题字体大小
        axis.text.x = element_text(size = 18, face = "bold"),       # x 轴刻度字体大小
        axis.text.y = element_text(size = 18, face = "bold"),       # y 轴刻度字体大小
        legend.title = element_text(size = 16, face = "bold"),      # 图例标题字体大小
        legend.text = element_text(size = 16, face = "bold")       # 图例标签字体         
      )+
  # 在图表中添加 p 值文本
  annotate("text", x = 1.5, y = max(data_combined$Count), label = paste("p-value =", round(p_value, 4)), size = 10, color = "red")
  # 添加显著性标记
  #annotate("text", x = 1.5, y = max(data_combined$Count), label = significance_label, size = 6, color = "red")


# 绘制密度图，添加 p 值注释
density_motif <- ggplot(data_combined, aes(x = Count, fill = Group, color = Group)) +
  geom_density(alpha = 0.5) +  # 密度图，alpha 设置透明度
  labs(title = "Density Plot of Seclect Motif vs Other Motif Counts (Log Scale)",
       x = "Log-Transformed Count",
       y = "Density") +
  scale_fill_manual(values = c("#BC3C29FF","#0072B5FF")) +  # 设置颜色
  scale_color_manual(values = c("#BC3C29FF","#0072B5FF")) +  # 设置曲线颜色
  theme_minimal() +
  theme(  # 覆盖字体大小
    plot.title = element_text(size = 20, face = "bold"), 
    axis.title.x = element_text(size = 20, face = "bold"),      # x 轴标题字体大小
    axis.title.y = element_text(size = 20, face = "bold"),      # y 轴标题字体大小
    axis.text.x = element_text(size = 18, face = "bold"),       # x 轴刻度字体大小
    axis.text.y = element_text(size = 18, face = "bold"),       # y 轴刻度字体大小
    legend.title = element_text(size = 16, face = "bold"),      # 图例标题字体大小
    legend.text = element_text(size = 16, face = "bold")       # 图例标签字体         
  )+
  # 在图表中添加 p 值文本
  annotate("text", x = max(data_combined$Count) - 1, 
           y = max(density(data_combined$Count)$y), 
           label = paste("p-value =", round(p_value, 4)), 
           size = 9, color = "red")
ggsave("select vs other violin_motif.pdf",plot=violin_motif,width = 12,height = 10)
ggsave("select vs other density_motif.pdf",plot=density_motif,width = 12,height = 10)

####################linux命令
seqkit locate -p "ATCGGGT" select_promoters.fa select_motif_position.txt
###################R###################
library(ggplot2)
library(gggenes)
# 加载数据
select_motif<- read.table("./select_motif_positions.txt", header = TRUE, stringsAsFactors = FALSE)
# 转换列名为 gggenes 所需格式
colnames(select_motif) <- c("gene", "patternName", "pattern","strand", "start", "end")
select_motif <- select_motif[,-7]
table(select_motif$gene)

select_mapping <- data.frame(
  region = c("chr1:83228741-83546044", "chr13:69760497-69905109", 
             "chr12:52168996-52226814", "chr9:89289174-89610503",
             "chr5:14803289-14839788", "chr8:157856149-158906138"),
  gene_name = c("gene1", "gene2", "gene3", "gene4", "gene5", "gene6")
)
# 假设 select_motif 数据框已有 gene 列.似乎不可以两个data.frame的都存在用于作为键的gene列，by="gene",不可以
select_motif <- merge(
  select_motif, 
  gene_mapping, 
  by.x = "gene", # 用 gene 列作为键
  by.y ="region",
  all.x = TRUE  # 确保所有 motif 都保留
)

motif_select_3 <- ggplot(select_motif, aes(xmin = start, xmax = end, y = gene_name, fill = pattern)) +
  geom_gene_arrow(color = "#BC3C29FF", size = 3) +  # 基因箭头表示 motif 范围
  facet_wrap(~ gene_name, scales = "free", ncol = 1) +
  theme_genes() +  # 使用 gggenes 的默认主题
  labs(
    fill = "Motif Pattern",
    title = "Motif Distribution in Genes",
    x = "Position",
    y = "Genes"
  ) +
  scale_fill_manual(values = c("AAAGTTT" = "#BC3C29FF")) +  # 填充颜色
  theme(  # 覆盖字体大小
    plot.title = element_text(size = 22, face = "bold"),
    axis.title.x = element_text(size = 20, face = "bold"),      # x 轴标题字体大小
    axis.title.y = element_text(size = 20, face = "bold"),      # y 轴标题字体大小
    axis.text.x = element_text(size = 18, face = "bold"),       # x 轴刻度字体大小
    axis.text.y = element_text(size = 18, face = "bold"),       # y 轴刻度字体大小
    legend.title = element_text(size = 16, face = "bold"),      # 图例标题字体大小
    legend.text = element_text(size = 16, face = "bold")       # 图例标签字体         
  )
#报错Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  : 
#Viewport has zero dimension(s)。其实是图显示面积不够
ggsave("select_motif.pdf",plot=motif_select_3,width = 12,height = 10)
motif_select_3

Violin plot

Density plot

motif position