内皮细胞中由H3K9乳酸化和HDAC2驱动的反馈回路调节vegf诱导的血管生成
对该文章中的Fig.3进行尝试复现。
文章简要讲解
研究背景
血管内皮生长因子（VEGF）是调控血管生成的重要因子，其机制涉及代谢重编程，尤其是糖酵解的增强。糖酵解增加导致乳酸积累，而乳酸作为新型表观遗传修饰——组蛋白乳酸化（如H3K9la）的底物，可能在血管生成中发挥关键作用。然而，乳酸如何通过组蛋白乳酸化影响内皮细胞行为，以及该修饰在血管生成中的具体机制尚不清楚。
研究方法与结果
在体外实验中，VEGF刺激显著提高了人视网膜微血管内皮细胞（HRMECs）中的H3K9la水平，并增强其增殖、迁移和成管能力。体内实验采用氧诱导视网膜病变（OIR）小鼠模型，结果显示H3K9la水平随病程进展上升，并在新生血管高峰期达到最高。使用糖酵解抑制剂可显著降低H3K9la水平并抑制血管生成。
CUT&Tag分析发现H3K9la富集于多个血管生成相关基因（如NECTIN1、TGFBR2、ABL1等）的启动子区域，ChIP-qPCR和qPCR进一步验证其对基因表达的促进作用。此外，研究发现VEGF下调HDAC2表达，而HDAC2过表达可抑制H3K9la水平和血管生成，提示H3K9la与HDAC2之间存在反馈调控机制。
简要介绍一下CUT&Tag、Chipseq、ATAC-seq的区别和联系
CUT&Tag
是一种靶向蛋白质-DNA相互作用的研究技术，与ChIP-seq类似，但不需要交联和免疫沉淀，而是通过抗体引导的切割实现DNA片段化，具有更高的分辨率和更低的背景信号。
[图片上传失败...(image-933f9c-1751184905263)]
ChIP-seq
染色质免疫共沉淀测序，用于研究特定蛋白（如组蛋白修饰、转录因子）在全基因组范围内的结合位点，需交联、沉淀和大规模测序。
[图片上传失败...(image-55daaa-1751184905263)]
ATAC-seq
转座酶可及性染色质测序，用于检测开放染色质区域，反映基因调控元件的可接近性，不依赖特异性抗体。
[图片上传失败...(image-8945aa-1751184905263)]
-
研究蛋白质-DNA相互作用
- 推荐使用 CUT&Tag：适用于微量样本，具有高信噪比，适合对样本量有限的实验。
- 或选择 ChIP-seq：具有广泛的适用性，适合多种实验条件和样本类型。
-
研究染色质开放性
- 推荐使用 ATAC-seq：可高效检测全基因组范围内的染色质开放区域，操作简便且灵敏度高。
-
探索未知调控因子
- 首先通过 ATAC-seq 筛查潜在的开放染色质区域，识别可能的功能性调控区域。
- 随后结合 CUT&Tag 或 ChIP-seq 对目标蛋白进行验证，明确其在开放区域的结合情况。
尝试复现
数据下载
创建项目目录并进入该目录：
mkdir -p Project/GSE248643 && cd Project/GSE248643
下载数据：
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE248nnn/GSE248643/suppl/GSE248643_RAW.tar -O GSE248643_RAW.tar
解压数据：
tar -xvf GSE248643_RAW.tar
bigWig（文件扩展名为.bw或.bigWig）是一种用于存储基因组连续数值数据的二进制格式。
下载参考文件并生成 gene.bed：
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_43/gencode.v43.annotation.gtf.gz
gunzip gencode.v43.annotation.gtf.gz
awk '$3 == "gene"' gencode.v43.annotation.gtf | \
awk 'BEGIN{FS="\t| |;";OFS="\t"} {
# 提取基因名和链信息
gene_name = "";
strand = $7;
for (i=9; i<=NF; i++) {
if ($i ~ /gene_name/) gene_name = $(i+1);
gsub(/"/, "", gene_name); # 去除引号
}
# 输出BED格式：chr, start, end, gene_name, score(0), strand
print $1, $4-1, $5, gene_name, "0", strand;
}' > genes.bed
head -n 3 genes.bed
输出示例：
chr1 11868 14409 DDX11L2 0 +
chr1 12009 13670 DDX11L1 0 +
chr1 14403 29570 WASH7P 0 -
或生成上下游3000bp的启动子区。
awk '$3 == "gene"' gencode.v43.annotation.gtf | \
awk 'BEGIN{FS="\t"; OFS="\t"} {
split($9, attrs, "; ");
gene_name = "";
for (i in attrs) {
if (attrs[i] ~ /gene_name/) {
split(attrs[i], tmp, "\"");
gene_name = tmp[2];
}
}
tss = ($7 == "+" ? $4 : $5);
print $1, tss-3000, tss+3000, gene_name, "0", $7;
}' > gene_promoters.bed
head -n 3 gene_promoters.bed
输出示例：
chr1 8869 14869 DDX11L2 0 +
chr1 9010 15010 DDX11L1 0 +
chr1 26570 32570 WASH7P 0 -
分析环境构建
创建 Conda 环境并指定 Python 版本：
conda create -n cuttag python=3.8
激活环境：
conda activate cuttag
安装 deeptools：
conda install -y -c bioconda deeptools
验证安装：
deeptools --version
输出示例：
deeptools 3.5.5
当前目录结构（示例）：
./
├── [1.5G] gencode.v43.annotation.gtf
├── [2.3M] gene_promoters.bed
├── [2.3M] genes.bed
├── [224M] GSE248643_RAW.tar
├── [113M] GSM7918343_Control.bw
├── [111M] GSM7918344_VEGF.bw
computeMatrix：基因组信号分布分析核心工具
computeMatrix 是 deepTools 工具包的核心模块，用于计算基因组信号在特定区域的分布，并生成标准化矩阵文件，为后续可视化（如热图、趋势图）提供数据支持。它支持两种主要计算模式：
| 模式 | 核心功能 | 适用场景 |
|---|---|---|
| reference-point | 以特定参考点（如 TSS、TES、peak 中心）对齐，计算上下游固定范围内的信号分布。 | 用于研究功能位点附近的信号富集，例如： - 转录起始位点（TSS）的染色质开放性 - ChIP-seq 中抗体富集效果评估 - Peak 中心的蛋白结合强度 |
| scale-regions | 将不同长度的区域（如基因体）缩放到相同长度，计算信号在缩放区域内的分布。 | 用于分析基因整体结构的信号趋势，例如： - 从 TSS 到 TES 的染色质状态变化 - 基因体区域的全局开放性或修饰模式 - 跨样本的基因结构信号比较 |
示例命令
1. reference-point 模式计算与绘图
computeMatrix reference-point \
--referencePoint TSS \
-b 3000 -a 3000 \
-S GSM7918343_Control.bw GSM7918344_VEGF.bw \
-p 35 \
-R genes.bed -o matrix.gz
plotHeatmap -m matrix.gz \
-out ExampleHeatmap1.png \
--whatToShow "plot, heatmap and colorbar" \
--colorMap Blues \
--zMin 0 --zMax 0.8 \
--heatmapHeight 15 \
--heatmapWidth 8 \
--dpi 300 \
--refPointLabel "summit"
2. scale-regions 模式计算与绘图
computeMatrix scale-regions \
-S GSM7918343_Control.bw GSM7918344_VEGF.bw \
-R genes.bed \
--regionBodyLength 6000 \
--startLabel "TSS" \
--endLabel "TES" \
--beforeRegionStartLength 3000 \
--afterRegionStartLength 3000 \
--binSize 50 \
-p 35 \
-o matrix_scale.gz
plotHeatmap -m matrix_scale.gz \
-out ExampleHeatmap2.png \
--whatToShow "plot, heatmap and colorbar" \
# 去除上方则线图
#--whatToShow "heatmap and colorbar" \
--colorMap Blues \
--zMin 0 --zMax 0.8 \
--heatmapHeight 15 \
--heatmapWidth 8 \
--averageType median \
--dpi 300 \
--refPointLabel "summit"
3. 折线图绘制
plotProfile -m matrix.gz \
-out ProfilePlot.png \
--perGroup \
--colors "#1f77b4" "#2ca02c" \
--yMin 0 --yMax 0.25 \
--refPointLabel "TSS" \
--startLabel "-3.0" --endLabel "3.0Kb" \
--plotHeight 10 --plotWidth 12 \
--dpi 300
数据预处理与峰注释分析
将 bigWig 文件转换为 BEDGraph 并提取峰
conda install -c bioconda ucsc-bigwigtobedgraph
bigWigToBedGraph GSM7918343_Control.bw Control.bedGraph
awk '$4 > 5 {print $1"\t"$2"\t"$3}' Control.bedGraph > Control_peaks.bed
bigWigToBedGraph GSM7918344_VEGF.bw VEGF.bedGraph
awk '$4 > 5 {print $1"\t"$2"\t"$3}' VEGF.bedGraph > VEGF_peaks.bed
R语言中进行峰注释与可视化
library(ChIPseeker)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
control_peaks <- readPeakFile("Control_peaks.bed")
vegf_peaks <- readPeakFile("VEGF_peaks.bed")
seqlevelsStyle(control_peaks) <- "UCSC"
seqlevelsStyle(vegf_peaks) <- "UCSC"
peakAnno_control <- annotatePeak(control_peaks, tssRegion = c(-3000, 3000),
TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)
peakAnno_vegf <- annotatePeak(vegf_peaks, tssRegion = c(-3000, 3000),
TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)
绘制韦恩图与峰注释饼图
library(VennDiagram)
venn.diagram(x = list(Control = control_ids, VEGF = test_ids),
filename = "H3K9la_venn.png",
fill = c("#5399bd", "#d7aa91"), cat.col = c("#5399bd", "#d7aa91"),
main = "H3K9la Binding Peaks Overlap")
# 绘制饼图
ppie = ggplot(combined_dist, aes(x = "", y = Frequency, fill = Feature)) +
geom_bar(stat = "identity", width = 1) +
coord_polar("y", start = 0) +
facet_wrap(~Group, ncol = 2) +
theme_void() +
scale_fill_brewer(palette = "Paired") +
geom_text(aes(label = paste0(round(Frequency, 1), "%")),
position = position_stack(vjust = 0.5), size = 3, color = "white")
ggsave("test.pie.png", plot = ppie, height = 4, width = 9, bg = "white")
差异峰富集
MACS2 bdgdiff：直接比较两组样本的peak信号强度，输出差异peak。
DiffBind（R包）：适用于有生物学重复的数据，通过统计检验（如DESeq2或edgeR）识别差异结合位点。
ChIPComp：整合多个样本的peak信息，计算差异结合显著性。
作者提供的文件有些少，提供的SRA文件是可以下载的，建议还是下载原始文件从上游开始处理
结果展示
- 热图与折线图成功复现文章趋势，与原文结论一致。
- 韦恩图展示了 Control 与 VEGF 组之间峰的重叠情况。
- 饼图展示了不同组中峰在基因组区域（如启动子、基因体、远端区域）的分布比例。
- 差异结合峰，还是从原始文件开始处理可能会更好
