使用deepTools绘制基因分布图：从BED文件到出版级可视化

1. 从一张图说起：为什么我们需要基因的“全局视野”？

在基因组学研究中，我们常常会问：我感兴趣的这些基因，它们在染色体上是随机分布的吗？还是倾向于聚集在特定的区域？比如，研究某个转录因子调控的靶基因，我们想知道它们是否富集在转录起始位点附近；或者分析一组在特定条件下差异表达的基因，看它们是否在染色体上成簇出现，这可能暗示着染色质高级结构或表观遗传调控的区域性影响。回答这些问题，一张直观的基因在基因组上的分布图（Genomic Distribution Plot）是必不可少的。

然而，从原始的基因列表（比如一个包含基因名或基因组坐标的文本文件）到一张信息丰富、可发表的分布图，中间往往隔着数据处理、格式转换、统计计算和可视化绘图等多个步骤。手动操作不仅繁琐，而且容易出错，重现性也差。这时候，一个强大、高效且被广泛认可的工具链就显得尤为重要。在生物信息学领域，deepTools正是为此而生的瑞士军刀之一。它并非专门用于绘制基因分布图，但其核心功能——将高通量测序数据（如ChIP-seq, ATAC-seq, RNA-seq）转化为各种汇总图（summary plots）和热图（heatmaps）——经过巧妙的“改装”，完全可以用来优雅地解决我们的问题。

简单来说，我们可以把每个基因看作一个“区间”（interval），其分布特征（如相对于转录起始位点TSS的上下游分布、在染色体上的绝对位置分布）就是我们需要可视化的“信号”。deepTools的computeMatrix和plotProfile/plotHeatmap组合，正是计算和绘制这种区间相关信号的利器。结合BED格式的基因区间文件和适当的参数设置，我们就能生成揭示基因分布模式的精美图表。本文将手把手带你走通这条从基因列表到出版级分布图的完整路径，并分享我在实际分析中积累的参数调优心得和避坑指南。

2. 核心工具链解析：deepTools如何为基因分布图赋能

deepTools是一套用 Python 编写的工具，主要用于处理高通量测序数据，其设计哲学是“将大数据转化为可解释的图”。对于绘制基因分布图，我们主要用到其中两个核心模块：computeMatrix和plotProfile/plotHeatmap。理解它们的工作原理，是灵活运用和排错的基础。

2.1 computeMatrix：从坐标到信号矩阵的引擎

computeMatrix是整个流程的计算核心。它的任务可以概括为：根据你提供的基因组区间列表（例如基因的TSS区域），从一个大范围的信号文件（例如全基因组的覆盖度文件）中，提取每个区间及其周边指定范围内的信号值，并将所有区间的信号对齐、缩放、平均，最终计算生成一个数值矩阵。

这个过程的输入和输出至关重要：

• 输入1：区间文件 (Regions File)：通常是一个BED格式的文件。每一行定义了一个基因组区间，例如chr1 1000 1500 geneA。对于基因分布图，这个文件通常包含你感兴趣的所有基因的坐标。一个关键技巧是，我们通常不直接用基因的整个body区域，而是用基因的转录起始位点 (TSS)作为代表点。因为许多调控事件（如转录因子结合、组蛋白修饰）都集中在TSS附近。我们可以用awk或专门的脚本从基因注释文件（如GTF）中提取TSS坐标并生成BED文件。
• 输入2：信号文件 (Score File)：这是一个描述全基因组范围内“信号强度”的文件。最常用的格式是bigWig (.bw)。bigWig文件是一种索引化的、压缩的二进制格式，可以高效地查询任意基因组区间的信号平均值。对于基因分布图，这里的“信号”需要根据你的科学问题来定义：
  • 如果你想看基因的绝对位置分布：你可以创建一个“虚拟”的均匀信号文件。例如，用bedtools genomecov生成一个全基因组每个碱基覆盖度为1的bedGraph文件，再转换为bigWig。这样，computeMatrix计算出的“信号”实际上就是区间的密度，经过后续绘图，就能反映出基因在基因组上的富集情况。
  • 如果你想看基因相对于某个表观标记的分布：那么信号文件就应该是相应的ChIP-seq或ATAC-seq数据的bigWig文件（经过标准化，如RPKM或CPM）。
• 核心参数与逻辑：
  • -b和-a: 分别定义每个区间上游 (upstream)和下游 (downstream)延伸多少碱基对(bp)来截取信号。例如，-b 3000 -a 3000会以每个区间的中心（或起点，取决于--referencePoint）为基准，向两侧各取3000bp。
  • --referencePoint: 定义区间的哪个位置作为对齐的参考点。可选center(区间中心)、TSS(转录起始位点，即BED的起点，需确保BED是TSS)、TES(转录终止位点)。对于基因分布图，最常用的是TSS，因为它是一个明确的、功能相关的位点。
  • --binSize: 将每个区间（包括上下游延伸区）划分成多少个小窗口（bins）来计算平均信号。较小的binSize(如10bp) 分辨率高，但数据量大；较大的binSize(如50bp) 更平滑，计算更快。通常50-100bp是一个不错的平衡点。
  • --missingDataAsZero: 如何处理信号文件中没有覆盖的区间？如果设为zero，则未覆盖区域信号计为0。这通常是你想要的，特别是使用虚拟均匀信号文件时。
  • --sortRegions和--sortUsing: 输出矩阵前如何对区间进行排序？descend降序，ascend升序，no不排序。排序可以让我们在热图中看到清晰的模式梯度。

computeMatrix的输出是一个二进制的.gz矩阵文件，它包含了所有区间、所有bin的信号值，以及相关的元数据（如区间名、坐标、排序信息等）。这个文件是下游绘图模块的输入。

2.2 plotProfile 与 plotHeatmap：从矩阵到洞察的画笔

拿到computeMatrix生成的矩阵文件后，我们可以用两个工具来可视化：

• plotProfile: 绘制线图 (line plot)。它会将所有区间的信号在每个bin上的值进行平均（或中位数等统计），然后绘制一条平均信号曲线，并通常带有阴影表示标准差或标准误。这张图能最清晰地展示信号的整体趋势。例如，你的基因集合是否在TSS上游表现出明显的信号峰？
• plotHeatmap: 绘制热图 (heatmap)。矩阵中的每一个值（区间x bin）对应热图中的一个色块。行代表一个基因区间，列代表基因组位置（从上游到下游）。这张图能同时展示整体趋势和个体差异。你可以看到是否所有基因都遵循同一模式，还是存在不同的亚群。结合排序，模式会更加明显。

两个绘图工具共享许多美化参数，如颜色 (--colors)、图例 (--legendLocation)、坐标轴标签 (--xAxisLabel,--yAxisLabel)、采样显示 (--plotType) 等。plotHeatmap还有额外的参数控制聚类 (--kmeans)、颜色标度 (--colorMap)、是否显示每行的基因标签 (--geneLabels) 等。

注意：deepTools的绘图是基于matplotlib的，其默认样式可能比较基础。为了得到出版级的图片，我们通常需要在生成图片后，用Inkscape、Adobe Illustrator或 Python 的matplotlib库直接进行二次美化（调整字体、线宽、图例位置等）。deepTools也提供--plotFileFormat参数输出svg或pdf矢量图，方便后期编辑。

3. 实战演练：从基因列表到分布图的完整流程

理论讲完，我们进入实战。假设我们有一个基因列表my_genes.txt，里面每行是一个基因名（如TP53,BRCA1）。我们的目标是看这些基因在基因组上的分布是否在TSS上游有特殊模式。

3.1 第一步：准备输入文件——BED与BigWig

1. 生成基因TSS的BED文件：我们需要一个基因名到基因组坐标的映射。最常用的来源是基因组注释文件，例如从GENCODE或Ensembl下载的GTF文件。

# 假设我们使用人类的GENCODE v44注释文件 (gencode.v44.annotation.gtf) # 使用 awk 提取所有基因的TSS坐标。注意：GTF中一个基因可能有多个转录本，这里取每个基因所有转录本中最小的起始位置作为代表TSS（对于+链基因，起始位置是start；对于-链基因，起始位置是end）。 awk 'BEGIN{OFS="\t"} $3=="gene" {gene_id=$10; gsub(/[";]/,"",gene_id); gene_name=$14; gsub(/[";]/,"",gene_name); if ($7=="+") {print $1, $4, $4+1, gene_name, ".", $7} else if ($7=="-") {print $1, $5-1, $5, gene_name, ".", $7}}' gencode.v44.annotation.gtf > all_genes_tss.bed

这个命令会生成一个包含所有基因TSS（1bp宽）的BED6文件（染色体，起始，终止，基因名，得分，链）。

接下来，从所有基因中筛选出我们感兴趣的基因：

# 假设 my_genes.txt 每行是基因名 grep -w -f my_genes.txt all_genes_tss.bed > my_genes_tss.bed

如果my_genes.txt里是其他标识符（如Ensembl ID），则需要根据GTF中的对应字段进行筛选。

2. 创建虚拟均匀信号BigWig文件：为了看基因的绝对分布，我们需要一个全基因组范围内信号恒定为1的bigWig文件。首先需要该基因组的染色体大小文件（例如hg38.chrom.sizes）。

# 方法1: 使用 bedtools genomecov (推荐) bedtools genomecov -bg -i my_genes_tss.bed -g hg38.chrom.sizes > uniform_coverage.bedgraph # 解释：-bg 输出bedGraph格式，-i 输入BED，-g 染色体大小文件。但这样生成的bedGraph只在我们有基因的位置有覆盖。 # 我们需要的是全基因组覆盖。一个技巧是创建一个包含所有位置的BED文件，但这不现实。 # 方法2: 更直接的方法，使用 wigToBigWig 工具包中的一个技巧，或者使用 ucsc 的 kent 工具。 # 这里介绍一个实用但取巧的方法：我们其实不需要一个真正的“均匀”信号，因为 computeMatrix 的 --missingDataAsZero 参数。 # 我们可以直接使用一个空的或无关的bigWig文件，并设置 --missingDataAsZero。但为了概念清晰，我们可以创建一个代表“基因存在”的信号。 # 更常见的做法是：我们直接分析基因的密度分布，这可以通过 deepTools 的 `computeMatrix` 对 BED 文件本身进行操作来实现，但需要另一种模式。 # 实际上，对于“基因在基因组上的分布图”，更常见的需求是看密度（每Mb有多少个基因），而不是信号强度。 # 因此，一个更合理的流程是：将基因组分成连续的窗口（如 1 Mb），计算每个窗口内我们目标基因的数量，然后用其他工具（如 R 的 ggplot2）做条形图或线图。 # 但如果我们坚持用 deepTools 的“信号”思路来模拟，可以这样做： # 用 awk 从染色体大小文件生成一个每1bp一个记录、得分全为1的庞大bedGraph，但这文件会巨大无比，不现实。 # 因此，我们必须重新审视目标：如果我们想看基因在染色体上的密度分布，deepTools 的 computeMatrix/plotProfile 并不是最直接的工具。 # 它更适合看“相对于某个点（如TSS）的信号分布”。对于绝对位置分布，应该用 bedtools 的 coverage 或 map。

看来这里遇到了一个关键概念区分。让我们回到原点：标题“基因在genome上的分布图”可能有两种理解：

1. 相对分布：基因集合的信号在某个参考点（如TSS）上下游的分布模式。这完美契合deepTools的computeMatrix模式。
2. 绝对分布/密度分布：基因在染色体不同区域（如着丝粒、端粒、染色体臂）的富集程度。这通常需要将基因组分箱 (binning)。

为了覆盖更常见的需求，我们调整示例：假设我们有一个组蛋白修饰（如H3K4me3）的ChIP-seq数据，我们已经有了其bigWig文件 (H3K4me3.bw)。我们想看看我们感兴趣的基因的启动子区域（TSS附近）的H3K4me3信号模式。这是一个非常经典且适合deepTools的分析。

那么，我们的输入文件就是：

• my_genes_tss.bed：我们目标基因的TSS坐标。
• H3K4me3.bw：全基因组H3K4me3 ChIP-seq信号文件（已标准化，如RPKM）。

3.2 第二步：运行computeMatrix计算信号矩阵

现在我们针对“相对分布”场景进行操作。我们想查看每个基因TSS上游3kb到下游3kb范围内的H3K4me3信号。

computeMatrix reference-point \ --referencePoint TSS \ -b 3000 -a 3000 \ -R my_genes_tss.bed \ -S H3K4me3.bw \ --binSize 50 \ --missingDataAsZero \ --sortRegions descend \ --sortUsing mean \ -o matrix_genes_H3K4me3_TSS.gz \ --outFileNameMatrix matrix_genes_H3K4me3_TSS.tab \ # 输出纯文本矩阵，可选，用于其他分析 --outFileSortedRegions sorted_regions_genes_H3K4me3.bed # 输出排序后的区间文件，可选

参数解释：

• reference-point: 模式，表示以参考点为中心进行分析。
• --referencePoint TSS: 以BED文件的起点作为参考点（我们的BED文件是1bp的TSS，所以正好）。
• -b 3000 -a 3000: 参考点上游和下游各取3000bp。
• -R: 输入区间BED文件。
• -S: 输入信号bigWig文件。
• --binSize 50: 将总共6000bp的区域分成 6000/50 = 120 个bins。
• --missingDataAsZero: 没有信号的地方记为0。
• --sortRegions descend --sortUsing mean: 按照所有bins的平均信号值从高到低对基因进行排序。
• -o: 输出压缩矩阵文件（主要输出）。
• --outFileNameMatrix和--outFileSortedRegions: 输出附加的文本文件，方便后续用其他工具处理或检查。

运行完成后，会生成matrix_genes_H3K4me3_TSS.gz。

3.3 第三步：绘制profile图与heatmap图

绘制平均信号曲线图 (Profile Plot):

plotProfile -m matrix_genes_H3K4me3_TSS.gz \ -o profile_plot.png \ --plotFileFormat png \ --perGroup \ # 如果有多组信号/多个样本，按组分别绘图。我们只有一组。 --colors blue \ --yAxisLabel "H3K4me3 signal (RPKM)" \ --refPointLabel "TSS" \ --plotTitle "H3K4me3 signal around TSS of target genes"

这会生成一张PNG图片，X轴是基因组位置（从-3kb到TSS再到+3kb），Y轴是平均信号强度。你通常会看到在TSS位置有一个尖锐的信号峰，这是活性基因启动子区域H3K4me3的典型特征。

绘制热图 (Heatmap):

plotHeatmap -m matrix_genes_H3K4me3_TSS.gz \ -o heatmap_plot.png \ --plotFileFormat png \ --colorMap RdBu_r \ # 使用红蓝渐变色系，_r表示反转 --yAxisLabel "Genes (sorted by mean signal)" \ --xAxisLabel "Distance from TSS (bp)" \ --refPointLabel "TSS" \ --legendLocation upper-right \ --dpi 300 # 提高分辨率

热图能更细致地展示每个基因的信号模式。排序后，信号强的基因在上方，信号弱的在下方。你可以清晰地看到信号模式的异质性。

实操心得：--colorMap的选择很有讲究。RdBu_r(红蓝)、viridis(黄-绿-蓝)、plasma(紫-黄) 等都是科学绘图常用的、对色盲友好的渐变色。避免使用jet，因为它虽然鲜艳但可能误导对数据相对大小的判断。

4. 高级技巧与深度参数调优

掌握了基础流程后，一些高级技巧和参数调优能让你的图更具洞察力和美感。

4.1 处理多个样本或条件组

如果你想比较不同条件下（如对照组 vs 处理组）同一组基因的信号分布，可以将多个bigWig文件同时输入给computeMatrix。

computeMatrix reference-point \ --referencePoint TSS \ -b 3000 -a 3000 \ -R my_genes_tss.bed \ -S control_H3K4me3.bw treated_H3K4me3.bw \ # 多个信号文件，用空格分隔 --binSize 50 \ --missingDataAsZero \ -o matrix_multi_sample.gz

然后在plotProfile中使用--perGroup参数，它会为每个样本画一条线，并用不同颜色区分。在plotHeatmap中，多个样本的信号会并排显示（默认是上下堆叠，可以用--regionsLabel和--samplesLabel来标注）。

4.2 使用scale-regions模式分析基因体信号

上面的reference-point模式专注于一个点（TSS）周围。如果你想分析整个基因体（从TSS到TES）以及上下游一定范围内的信号，可以使用scale-regions模式。

computeMatrix scale-regions \ -R my_genes_body.bed \ # 这个BED文件需要是基因的完整区域，而不仅仅是TSS -S H3K27ac.bw \ # 例如分析增强子标记H3K27ac在基因体的分布 -b 2000 -a 2000 \ # 基因体上下游额外延伸2kb --regionBodyLength 5000 \ # 将基因体本身缩放到一个固定长度（如5000bp），便于不同长度基因的比较 --binSize 50 \ --missingDataAsZero \ -o matrix_genes_body.gz

在这种模式下，X轴通常被分为三部分：上游区、基因体（缩放后）、下游区。这对于研究像RNA聚合酶II（Pol II）或某些组蛋白修饰在转录单元内的分布非常有用。

4.3 绘图美化与输出控制

• 输出矢量图：将--plotFileFormat设为svg或pdf，方便在Adobe Illustrator或Inkscape中无损编辑和组合。
• 调整图片尺寸：使用--plotWidth和--plotHeight参数（单位是英寸）。发表文章时通常需要特定宽度（如单栏8cm，双栏17cm）。
• 自定义颜色：--colors参数接受颜色名称（如red,blue）或十六进制码（如#FF0000,#0000FF）。对于多组数据，按顺序指定，如--colors red blue green。
• 修改坐标轴：--yMin和--yMax可以固定Y轴范围，使得多张图之间可以比较。--xAxisLabel和--yAxisLabel用于设置轴标签。
• 关闭默认图例或标题：使用--legendLocation none关闭图例，--plotTitle ""清除标题，以便在后期软件中添加更统一的格式。

4.4 性能优化与大数据集处理

当基因列表很大（>10,000）或信号文件很多时，computeMatrix可能会消耗大量内存和时间。

• 增大--binSize：从50bp增加到100bp或200bp，可以显著减少计算量和矩阵大小。
• 使用--smartLabels：当区间文件有重复名时，自动处理标签。
• 分而治之：如果内存不足，可以考虑将大的BED文件拆分成多个小文件，分别运行computeMatrix，然后使用plotProfile和plotHeatmap的-m参数同时接受多个矩阵文件进行绘图（但需要注意样本顺序一致）。
• 利用多核：computeMatrix支持--numberOfProcessors或-p参数来指定使用的CPU核心数，可以加速计算。

5. 常见问题排查与避坑指南

即使按照流程操作，也可能会遇到各种问题。以下是一些常见坑点及其解决方案。

5.1 错误：Error: The bigWig file appears to be malformed!或Error: Received error code -1

• 可能原因1：bigWig文件索引损坏或缺失。bigWig文件需要配套的.bwi索引文件。确保两者在同一目录下，且文件名正确（例如file.bw和file.bw.bwi或file.bwi）。你可以用bigWigInfo工具检查文件。
• 可能原因2：bigWig文件的染色体命名与BED文件不匹配。这是最常见的问题。BED文件中用的是chr1，而bigWig文件中可能用的是1（无chr前缀），反之亦然。使用head命令查看BED文件的前几行，用bigWigInfo查看bigWig文件的染色体列表，确保一致。

  # 查看BED文件染色体命名 head -n 5 my_genes_tss.bed # 查看bigWig文件染色体列表 (需要 ucsc-kent 工具包中的 bigWigInfo) bigWigInfo H3K4me3.bw | head -20

  解决方案：统一命名规范。可以使用sed命令为BED文件添加或去除chr前缀。

  # 为BED文件添加chr前缀（如果bigWig有chr） sed -i 's/^/chr/' my_genes_tss.bed # 或者，如果bigWig没有chr，而BED有，则去除chr前缀 sed -i 's/^chr//' my_genes_tss.bed

  注意：-i参数会直接修改原文件，操作前建议备份。也可以使用sed 's/^chr//' input.bed > output.bed生成新文件。

5.2 绘图时Y轴范围不合理或图形扭曲

• 现象：Profile图的Y轴从0开始，但信号值都是几十上百，导致曲线挤在顶部；或者热图的颜色条范围不合适，使得对比度很差。
• 解决方案：
  • 对于plotProfile，使用--yMin和--yMax手动设置Y轴范围。可以先从输出的矩阵文本文件 (matrix.tab) 或plotProfile的标准输出中查看信号的大致范围。
  • 对于plotHeatmap，使用--zMin和--zMax来设置颜色映射的值范围。例如，--zMin 0 --zMax 10会将所有低于0的值映射为最小颜色，高于10的值映射为最大颜色，0-10之间线性映射。这能有效增强对比度，突出差异。也可以使用--whatToShow调整显示内容（如heatmap and colorbar）。

5.3 热图中基因标签重叠或无法显示

• 现象：当基因数量很多时（比如超过100个），在热图左侧显示所有基因名会导致文字重叠，无法辨认。
• 解决方案：
  • 使用--geneLabels none关闭基因标签显示。
  • 如果必须显示，可以尝试增大图片高度 (--plotHeight)，但效果有限。
  • 更实用的做法：热图主要用于观察整体模式，而非识别单个基因。如果需要识别特定基因，可以在排序后的区间文件 (sorted_regions.bed) 中找到其排名，或者将热图与后续的基因功能分析结合。在论文中，通常只展示具有代表性的部分基因或聚类。

5.4 computeMatrix运行缓慢或内存不足

• 原因：区间太多、信号文件太大、binSize太小、上下游范围太大。
• 优化策略：
  1. 过滤区间：如果基因列表很大，考虑根据表达量、显著性等指标筛选出最感兴趣的子集进行分析。
  2. 调整参数：增大--binSize(如从10调到50)，减小-b和-a的范围（如果不是必须分析很远的区域）。
  3. 使用--blackListFileName：如果分析中需要排除某些区域（如高重复序列区、黑名单区域），提前指定黑名单文件，deepTools会在计算时跳过这些区域，有时能减少计算量。
  4. 增加内存和CPU：在计算集群上提交任务，指定更多的内存（如--mem-per-cpu=10G）和使用多核 (-p 8)。

5.5 生成的图片风格不符合期刊要求

• 问题：deepTools默认的字体、线宽、图例样式可能比较简陋。
• 终极解决方案：输出矢量图 (svg/pdf)，然后用专业矢量图形软件（如Adobe Illustrator、Inkscape、Affinity Designer）进行美化。这是发表前几乎必须的一步。你可以在这些软件中轻松修改字体为 Arial 或 Helvetica，调整字号，加粗线条，移动图例，添加面板标签（如 A, B, C）等。
• 进阶方案：deepTools的绘图函数是基于matplotlib的。你可以通过创建自定义的matplotlib样式文件 (rcParams)，并在plotProfile或plotHeatmap中通过--plotFile参数指定一个自定义的 Python 绘图脚本来实现更精细的控制，但这需要一定的 Python 编程能力。

6. 超越deepTools：其他可视化思路与工具

虽然deepTools非常强大，但它主要擅长“相对分布”。对于“绝对分布”或更复杂的可视化需求，可能需要结合其他工具。

• 染色体圈图 (Circos Plot)：如果你想展示基因在多条染色体上的分布，以及它们之间的关联（如共表达、互作），Circos 圈图非常强大，但学习曲线陡峭。R的circlize包是一个不错的替代。
• 曼哈顿图 (Manhattan Plot)：通常用于全基因组关联分析 (GWAS)，但也可以用来展示基因（或其他特征）在染色体上的分布密度。每个点代表一个基因组窗口（如 1 Mb），Y轴是该窗口内目标基因的数量或密度。可以用R的qqman或ggplot2包绘制。
• 基因密度条形图：使用bedtools的makewindows和coverage功能，将基因组划分为固定大小的窗口，计算每个窗口内目标基因的个数，然后用R/Python绘制条形图或折线图。这是查看基因在染色体上宏观分布的最直接方法。

  # 示例：计算目标基因在1Mb窗口内的覆盖度 bedtools makewindows -g hg38.chrom.sizes -w 1000000 > genome_1mb_windows.bed bedtools coverage -a genome_1mb_windows.bed -b my_genes_tss.bed > gene_coverage_per_1mb.bed

  得到的文件包含每个1Mb窗口内目标基因的计数，接下来可以用任何绘图工具进行可视化。

我个人在项目中的体会是，没有一种工具是万能的。deepTools在解决“围绕某个特征点的信号分布”这类问题上效率极高，且能产生可直接用于发表的中间结果图。但对于更全局的、绝对位置的分布观察，通常需要结合bedtools进行前期数据处理，再用R或Python进行灵活的可视化。理解每种工具的设计初衷和边界，根据具体的科学问题选择合适的工具链组合，才是高效生物信息学分析的关键。最后，无论使用什么工具，确保你的分析流程清晰、可重复，并且每一步的参数和结果都有据可查，这才是产出可靠科研结果的基石。