从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, headerFALSE)如果数据没有表头一定要设置headerFALSE否则第一行数据会被误认为列名。曾经有同事因为这个参数没设置导致后续分析全部错位白白浪费了一周时间。2.3 数据格式的精简转换CMplot只需要三列数据SNP名称、染色体编号、物理位置。用dplyr的select操作可以轻松完成library(dplyr) clean_data - map_data %% select(SNPV2, ChromosomeV1, PositionV4)注意染色体编号最好是数字格式。如果用的是chr1这种格式需要先用正则表达式处理clean_data$Chromosome - gsub(chr, , clean_data$Chromosome)3. CMplot的核心参数详解3.1 基础绘图一行代码出图最简单的密度图只需要一行代码CMplot(clean_data, plot.typed, bin.size1e6)但这样出来的图可能不太美观。让我分享几个实战中总结的参数技巧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.typed, bin.size1e6, colc(darkgreen, yellow, red))3.3 输出设置与排版file.outputTRUE时自动保存图片FALSE时在R中显示file格式支持jpg,pdf,tiff等。投稿论文推荐tiff格式dpi印刷质量建议300以上屏幕展示72就够了memo可以在文件名后添加备注比如final_versionCMplot(clean_data, plot.typed, bin.size1e6, filetiff, memodensity, dpi300)4. 高级技巧与问题排查4.1 处理特殊染色体有些数据包含性染色体或线粒体DNA将X染色体编码为23Y为24MT为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-bed4.3 常见报错解决Error in plot.window(...)通常是因为染色体编号包含字符确保转换为数值型Need at least 3 SNPs to plot检查数据是否成功读入图形元素重叠调整width和height参数或减小bin.size5. 解读密度图从图案发现隐藏问题一张好的密度图不仅能展示数据还能反映潜在问题全局密度不均可能提示测序深度不一致或样本污染局部尖峰可能是重复区域或比对错误染色体末端缺失常见于某些芯片设计缺陷曾经有个案例某项目的1号染色体出现周期性密度波动后来发现是DNA提取时存在技术重复。这种问题在传统QC指标中很难发现却在密度图上显露无遗。最后分享一个我常用的完整代码模板保存为.R文件随时调用library(data.table) library(dplyr) library(CMplot) # 数据读取 map_data - fread(your_data.map, headerF) # 数据清洗 plot_data - map_data %% select(SNPV2, ChromosomeV1, PositionV4) %% mutate(Chromosomeas.numeric(gsub(chr, , Chromosome))) # 绘图输出 CMplot(plot_data, plot.typed, bin.size1e6, colc(grey90, grey50, black), filepdf, memofinal, dpi300, width20, height12, file.outputTRUE)
)
