从PLINK到CMplot：三步绘制高颜值SNP密度图
1. 从PLINK数据到SNP密度图：为什么需要可视化
做基因组分析的朋友都知道，拿到原始数据后的第一件事就是检查数据质量。我刚开始做GWAS研究时，导师问的第一个问题就是："你的SNP在染色体上分布均匀吗？"当时我就懵了——难道SNP不是随机分布的吗？后来才发现，很多测序平台会有偏好性，某些染色体区域可能出现SNP堆积或缺失，这会直接影响后续分析的可靠性。
这时候SNP密度图就派上用场了。它能直观展示每条染色体上每1Mb区间内的SNP数量，用颜色梯度反映密度高低。就像给基因组做"CT扫描"，哪里密集哪里稀疏一目了然。我后来养成了习惯，拿到新数据必先画密度图，曾经还因此发现过样本标签错配的问题（某些样本的SNP密度模式明显异常）。
传统方法需要写几十行代码，直到发现了CMplot这个R包。它把整个流程简化到了三步：读数据、调参数、出图。最让我惊喜的是，它可以直接处理PLINK的.map文件格式——这是大多数基因型分析产出的标准格式，省去了繁琐的格式转换步骤。
2. 数据准备：从PLINK到R的完美衔接
2.1 理解PLINK的.map文件
PLINK是遗传学研究的"瑞士军刀"，它的.map文件包含四个核心字段：
染色体编号 SNP名称 遗传距离 物理位置比如：
1 rs12345 0 1000000 1 rs67890 0 2000000 2 rs11111 0 1500000重点在于第四列的物理位置（单位是碱基对），这是画密度图的关键。有时候数据可能缺少遗传距离列（第三列全为0），这完全不影响密度图绘制，因为SNP密度只与物理位置有关。
2.2 用data.table高效读取大数据
我测试过多种读取方式，当.map文件超过1GB时，基础R的read.table()会慢到怀疑人生。这里推荐data.table包的fread()函数，速度能快10倍以上：
library(data.table) map_data <- fread("genotype.map", header=FALSE)如果数据没有表头，一定要设置header=FALSE，否则第一行数据会被误认为列名。曾经有同事因为这个参数没设置，导致后续分析全部错位，白白浪费了一周时间。
2.3 数据格式的精简转换
CMplot只需要三列数据：SNP名称、染色体编号、物理位置。用dplyr的select操作可以轻松完成：
library(dplyr) clean_data <- map_data %>% select(SNP=V2, Chromosome=V1, Position=V4)注意染色体编号最好是数字格式。如果用的是"chr1"这种格式，需要先用正则表达式处理：
clean_data$Chromosome <- gsub("chr", "", clean_data$Chromosome)3. CMplot的核心参数详解
3.1 基础绘图：一行代码出图
最简单的密度图只需要一行代码：
CMplot(clean_data, plot.type="d", bin.size=1e6)但这样出来的图可能不太美观。让我分享几个实战中总结的参数技巧：
3.2 颜色与分箱的艺术
- bin.size：控制统计窗口大小，1e6表示1Mb。对于高密度芯片（如WGS），可以尝试500kb；对于低密度芯片（如GWAS芯片），2Mb可能更合适
- col：设置颜色梯度，从左到右对应低密度到高密度。推荐几组实测好看的配色：
- 保守风格：
c("darkgreen", "yellow", "red") - 印刷友好：
c("grey90", "grey50", "black") - 色盲友好：
c("#E69F00", "#56B4E9", "#009E73")
- 保守风格：
CMplot(clean_data, plot.type="d", bin.size=1e6, col=c("darkgreen", "yellow", "red"))3.3 输出设置与排版
- file.output：TRUE时自动保存图片，FALSE时在R中显示
- file：格式支持"jpg","pdf","tiff"等。投稿论文推荐tiff格式
- dpi：印刷质量建议300以上，屏幕展示72就够了
- memo：可以在文件名后添加备注，比如"final_version"
CMplot(clean_data, plot.type="d", bin.size=1e6, file="tiff", memo="density", dpi=300)4. 高级技巧与问题排查
4.1 处理特殊染色体
有些数据包含性染色体或线粒体DNA：
- 将X染色体编码为23，Y为24，MT为25
- 或者先过滤掉这些染色体：
filter(clean_data, Chromosome %in% 1:22)
4.2 超大数据的处理技巧
当数据超过100万SNP时：
- 先随机抽样10%画图看趋势：
sample_n(clean_data, nrow(clean_data)*0.1) - 或者先用PLINK做区域过滤：
plink --bfile data --chr 1-22 --make-bed
4.3 常见报错解决
- "Error in plot.window(...)"：通常是因为染色体编号包含字符，确保转换为数值型
- "Need at least 3 SNPs to plot"：检查数据是否成功读入
- 图形元素重叠：调整
width和height参数，或减小bin.size
5. 解读密度图：从图案发现隐藏问题
一张好的密度图不仅能展示数据，还能反映潜在问题：
- 全局密度不均：可能提示测序深度不一致或样本污染
- 局部尖峰：可能是重复区域或比对错误
- 染色体末端缺失：常见于某些芯片设计缺陷
曾经有个案例：某项目的1号染色体出现周期性密度波动，后来发现是DNA提取时存在技术重复。这种问题在传统QC指标中很难发现，却在密度图上显露无遗。
最后分享一个我常用的完整代码模板，保存为.R文件随时调用：
library(data.table) library(dplyr) library(CMplot) # 数据读取 map_data <- fread("your_data.map", header=F) # 数据清洗 plot_data <- map_data %>% select(SNP=V2, Chromosome=V1, Position=V4) %>% mutate(Chromosome=as.numeric(gsub("chr", "", Chromosome))) # 绘图输出 CMplot(plot_data, plot.type="d", bin.size=1e6, col=c("grey90", "grey50", "black"), file="pdf", memo="final", dpi=300, width=20, height=12, file.output=TRUE)