说实话我第一次用PHARP给猪芯片数据做基因型填充的时候真正卡住我的不是软件本身而是把数据喂进去之前的那一步“文件准备”。当时我手里是一份从测序公司导出的芯片Excel报告六万多SNP位点几百头猪的样本但PHARP一启动就报错报错信息还特别不友好根本看不出是哪里的问题。后来反复排查才发现问题出在最基础的地方染色体编号写得不统一、样本ID里有隐藏空格、等位基因编码方向反了。这些坑单独拎出来都不算大但堆在一起就足够让人抓狂一整周。这篇内容就是围绕“PHARP对猪芯片数据进行填充时的文件准备”展开的。我会把一份真实常见的芯片原始数据一步步改造成PHARP能直接吃进去的输入文件包括格式要求、字段含义、常见报错和对应的排查思路。不管你是第一次接触基因型填充还是已经被各种文件格式折磨过几轮这篇文章都值得你花十分钟读完大概率能帮你少走不少弯路。1. 为什么“填充”第一步卡在文件准备而不是软件本身1.1 一份真实猪芯片数据通常以什么形式存在猪的芯片数据尤其是从商业化检测平台拿回来的结果最常见的形态是Excel或CSV表格。每一行是一个样本在某一个SNP位点上的分型结果列通常包括样本编号、位点名称、染色体、物理位置、两个等位基因Allele1和Allele2。听起来挺规整对吧但实际拿到手你会发现不同平台导出的字段顺序、命名方式、编码风格差异非常大。有些用“A/G”这种双字母表示基因型有些用“AA/AG/GG”三字母还有直接用数值编码的。更麻烦的是很多表格里会夹带一些辅助列比如分型质量分数、GC值、重复样本的标记等。PHARP这种填充工具并不关心这些额外的花活它只认特定结构的最小文件集合。所以文件准备的第一件事就是从这一堆乱糟糟的原始数据里提取出PHARP真正需要的核心信息。1.2 “填充”到底在填什么缺失基因型填补的机制简述猪芯片数据填充本质上做的是基因型填补imputation。芯片本身覆盖的位点有限比如常见的GGP PorcineHD芯片大约有80K个标记但参考猪基因组里已知的变异位点数远不止这些。填充的目的就是利用参考群体的单倍型信息和样本间的连锁不平衡LD关系预测出样本在未被芯片直接检测到的位点上的基因型。我习惯用一个简单的类比来解释你有一幅拼图芯片相当于只给了你其中一部分碎片而填充算法则根据隔壁几幅完整拼图的图案规律推断出你手里缺的那些位置应该是什么颜色。PHARP做的就是这件事。但它需要三个基础输入一是已经分型好的样本基因型数据二是位点位置信息三是样本之间的对应关系也就是ID体系。这些信息如果文件里没写对填充算法再聪明也白搭。1.3 文件准备的本质把平台语言翻译成软件语言为什么单独把“文件准备”拎出来说因为填充工具本身往往是“死脑筋”的。PHARP读文件的时候不会像人一样去猜你的意思它读第一列就认为是样本家族ID读第二列就认为是样本个体ID读错位了不会警告你只会给出一个让人摸不着头脑的运行日志或者干脆静默输出一个全是Na的结果。所以文件准备的本质就是一次严格的“翻译”工作把你手里平台格式的数据翻译成PHARP定义的输入协议。这个环节做得越干净后面运行填充就越顺利。我见过太多人花大量时间去调参数、研究算法结果回头发现错误出在某个文件里一个不起眼的空格上这种低级错误是最亏的。2. PHARP能识别的输入长什么样字段顺序、编码与系谱文件2.1 PHARP输入文件的最小集合不同版本的PHARP输入文件的具体定义会有一些差异我用的是基于PED/MAP体系的版本这也是绝大多数填充软件通用的格式。最小集合包含两个文件基因型文件通常后缀是.ped样本信息和基因型数据混排。位点信息文件通常后缀是.map位点的染色体、名称和物理位置。如果你的数据里还涉及家系信息猪场的数据一般都有系谱那么还需要一个系谱文件把样本之间的亲缘关系补上。亲缘关系对填充准确率有很大影响尤其是在家系结构明显的群体里父本母本的信息能明显提升低密度芯片的填充效果。表里的内容比较枯燥但它值得你逐字检查文件必需列说明.ped前6列固定家族ID、个体ID、父本ID、母本ID、性别、表型.ped第7列起每两个等位基因表示一个位点顺序与.map严格一致.map4列固定染色体、位点ID、遗传距离、物理位置系谱文件3列个体ID、父本ID、母本ID2.2 基因型编码规则从双字母到数值编码PED文件里一个SNP位点的基因型用两个字符表示比如“AG”。两个字符相同就是纯合AA或GG不同就是杂合AG。缺失用“0 0”表示。这里有一个非常容易错的点字符之间到底有没有空格严格来说PED格式里每个等位基因之间可以用空格分隔但也可以连着写。不过为了兼容性和后续处理我强烈建议统一用“单个字符中间不加空格”的方式也就是一个样本在一个位点上的基因型整体写作“AA”“AG”“GG”“00”这样的两字符形式。如果你混用了“A G”和“AG”某些解析器会崩。PHARP内部处理时通常会把双字母基因型进一步转成数值编码比如AA记作0AG记作1GG记作2缺失记作-9或NA。这个转换不一定需要你手动做只要你提供的PED文件标准PHARP会自己处理。但如果你拿到的原始数据本身就是数值编码比如已经分好0/1/2了那你需要确认PHARP版本是否支持直接读数值矩阵支持的话反而省事直接跳过PED环节。2.3 系谱文件与表型文件怎么和基因型关联系谱和表型文件的核心作用是解释样本ID之间的关系。好多人在这一步翻车不是因为格式写得不对而是因为ID对不上。比如基因型文件里样本ID写的是耳标号“HN-2024-0135”系谱文件里却写成了“135”PHARP运行时会认为这是完全不同的个体所有亲缘关系全部断裂。关联逻辑很简单文件之间的ID必须“完全字符串匹配”。不存在模糊匹配这一说哪怕多一个空格、多一个字母都会被当成不同个体。我自己的做法是在做任何格式转换之前先单独写一个脚本把全部文件的样本ID提取出来做交集比对确保每一个要参与填充的样本ID在所有文件里都长得一模一样。表型文件对填充本身不是必需品但如果你后续要做全基因组关联分析GWAS那么表型文件建议和基因型文件一起准备。PED文件第6列就是表型列二分类性状写1和2连续性状写具体数值缺失写-9。3. 从芯片原始数据到PHARP输入文件的转换实操3.1 用PLINK做第一轮数据清洗拿到原始芯片数据后我一般先不急着转换而是用PLINK做第一轮清洗。PLINK是一个开源的全基因组关联分析工具集虽然它的主要用途是GWAS但它的数据清洗和格式转换能力非常实用处理PED/MAP和BED/BIM/FAM格式是它的老本行。第一步把原始数据整理成PLINK能识别的格式。如果公司给你的数据是Excel导出的CSV你需要先把它转成常规的PED/MAP这一步我用R来完成后面细说。如果公司直接给了你PLINK格式的文件那就省事了直接进入清洗环节。假设你已经有了一个名为“rawdata”的PED/MAP文件基础版常用的清洗命令如下# 1. 查看缺失率统计 plink --file rawdata --missing --out check_missing # 2. 剔除缺失率过高的SNP位点--geno 0.1 指位点缺失率超过10%则剔除 plink --file rawdata --geno 0.1 --recode --out clean_geno # 3. 剔除缺失率过高的样本--mind 0.1 指样本缺失率超过10%则剔除 plink --file clean_geno --mind 0.1 --recode --out clean_geno_mind # 4. 输出清洗后的PED和MAP作为PHARP的候选输入 plink --file clean_geno_mind --recode --out pharp_candidate这里要提醒一句PLINK处理猪的常染色体没问题但X、Y和线粒体基因组的处理逻辑和人类不完全一样。默认情况下PLINK会把X、Y、MT识别为自己的内部编码你的MAP文件里如果是“X”“Y”“MT”需要确认PLINK版本能正确识别或者干脆在清洗之前就把性染色体和常染色体分开处理常染色体走PLINK性染色体单独写脚本处理。3.2 用R脚本把Excel导出的数据改造成PED/MAP大部分公司的芯片交付文件是类似“样本ID、位点ID、Allele1、Allele2”的长表结构每一行是一个样本-位点组合你需要把它变成宽表每一行是一个样本每一列是一个位点。这个操作用R的data.table包最顺手。下面这段代码我一直在用改改路径就能跑library(data.table) # 读取原始交付文件 raw - fread(chip_raw.csv) # 假设列名SampleID, SNPID, Chr, Pos, Allele1, Allele2 # 生成基因型列按A/G双字母格式 raw[, Genotype : paste0(Allele1, Allele2)] # 缺失分型统一写成00 raw[is.na(Allele1) | is.na(Allele2), Genotype : 00] # 把长表转成宽表每行一个样本每列一个位点 geno_wide - dcast(raw, SampleID ~ SNPID, value.var Genotype) # 构建PED文件前6列 ped - data.table( FID 0, # 家族ID无家系信息时统一填0 IID geno_wide$SampleID, PID 0, # 父本ID后续用系谱补 MID 0, # 母本ID Sex 0, # 0未知1公2母 Pheno -9 # 表型缺失 ) # 合并基因型列 ped - cbind(ped, geno_wide[, -SampleID]) # 构建MAP文件 map - unique(raw[, .(SNPID, Chr, Pos)]) map[, GeneticDist : 0] # 很多芯片没有遗传距离填0即可 map - map[, .(Chr, SNPID, GeneticDist, Pos)] # 输出 write.table(ped, pharp_input.ped, row.names FALSE, col.names FALSE, quote FALSE) write.table(map, pharp_input.map, row.names FALSE, col.names FALSE, quote FALSE)注意事项dcast这一步如果同一个样本在同一个位点出现多次比如重复分型会报警告。如果数据里存在重复行一定要先处理掉否则宽表会多出很多带后缀的列污染后续所有分析。3.3 小样本试跑验证文件可用性的最快路径文件转换完成后不要急着把全部几百个样本直接灌给PHARP。我强烈建议先抽5到10个样本单独生成一份迷你输入文件试跑一遍完整流程。这一步能帮你过滤掉绝大部分格式错误还不会浪费太多时间。试跑的逻辑很简单如果小样本能顺利跑完那说明文件的整体框架没问题如果报错错误定位也容易因为数据量小你可以逐行检查。等小样本验证通过后再跑全量数据你就只需要担心计算资源和运行时间的问题了。# 用head命令取前8行PED文件的前几列这里假设你有8个样本 head -n 8 pharp_input.ped | cut -f1-10当然更稳妥的是在R里用sample函数随机抽样本ID再根据ID从全量数据里子集出来避免破坏文件结构。子集文件跑通PHARP全流程后你至少能确认三件事文件格式正确、行为标准、算法能启动。这三件事确认了文件准备这个阶段才算真正过关。4. 文件准备里容易翻车的四个细节ID不一致、等位基因链方向、染色体命名与缺失率4.1 样本ID不一致最容易出错的隐形地雷我见过的最常见错误不是格式错而是ID不一致。芯片数据文件里样本叫“HN-2024-0135”系谱文件里叫“20240135”Excel表格里则被自动改名成了“1.35E07”。这些ID看似“差不多”但在PHARP眼里完全不是同一个个体。最典型的坑是Excel自动转换你知道Excel会把长数字ID转成科学计数法吧十几位的个体编号或条码编号被Excel一保存就变了样。这个问题隐蔽性极高因为你在Excel里看到的可能还是原有显示格式但底层存储的数值已经变了。解决办法很简单拿到任何平台交付的Excel文件第一时间另存为CSV然后用文本编辑器或R、Python打开检查ID列确认没有科学计数法和隐藏字符。如果有条件建议从一开始就要求平台同时交付一份文本格式的ID清单这份清单就是你所有文件ID比对的基准。4.2 等位基因链方向A/T与G/C的“镜像陷阱”等位基因链方向这个问题做过人类芯片数据的人应该不陌生但在猪的数据里同样存在而且更容易被忽略。它的本质是同一个SNP位点正链forward strand和负链reverse strand的等位基因是互补的。比如正链上是A/G负链上就变成了T/C。如果你的数据来源不同比如一部分样本用的是GGP芯片另一部分用的是PorcineSNP80芯片或者跟参考面板的数据做了合并那么同一个位点可能出现“A/G”和“T/C”两种记录方式。如果直接把数据合并进填充流程PHARP会把这个位点当成完全不同的两个位点或者更糟糕把它们当成同一个位点但全是杂合结果完全不可用。解决办法是在文件准备阶段做一个链方向一致性检查。理想情况下你的参考面板文件里已经标注了链方向你需要把待填充数据的等位基因编码调整到与参考面板一致。如果一个位点在猪参考基因组Sscrofa11.1上能查到明确的正链坐标那直接以参考序列为准校正即可。实操中我建议把全部SNP位点的等位基因编码都转成与参考基因组正链一致后再进入填充流程这样最省心。# 示例把负链等位基因转成正链 # 互补映射 complement - c(A T, T A, G C, C G) # 假设SNP数据里有个allele列如果链方向标记为-则取互补 raw[Strand -, Allele1 : complement[Allele1]] raw[Strand -, Allele2 : complement[Allele2]]但要注意A/T和G/C这种纯互补的位点正负链无法简单区分这类位点在跨平台合并时一般建议直接剔除或者用额外的单倍型信息来判断。好在猪芯片数据里这类位点比例不高删掉对总体填充效果影响有限。4.3 染色体命名与物理位置的版本问题染色体命名不统一是另一个高频坑。猪有18条常染色体加X、Y和线粒体基因组。不同数据源里染色体可能写作“1”“chr1”“Sus scrofa chromosome 1”或者干脆用“1.1”“1.2”这种带组装版本后缀的格式。PHARP内部对染色体名称的处理是字符串级别的它不会去做智能归一化。所以你在准备MAP文件时必须把所有染色体统一成同一种写法。我建议统一写成数字1到18X写成“X”Y写成“Y”线粒体写成“MT”。当然有些版本的PHARP只接受数字染色体编号这时候就需要把X映射为19、Y映射为20、MT映射为21之类的固定编号。这个映射方式一定要写进你的内部说明文档里免得下次自己都忘了。物理位置也要注意版本问题。猪的参考基因组已经更新了好几版Sscrofa10.2和Sscrofa11.1之间的坐标差异不小。如果你的多个数据源使用了不同版本的基因组坐标必须在合并之前做坐标转换否则位点顺序会乱填充时LD关系也会算错。这个问题不那么明显因为文件不会报错但最后的填充结果准确率会莫名其妙地下降。4.4 缺失率阈值哪些位点值得填充哪些应该直接剔除文件准备阶段还要想清楚一个问题是不是所有缺失位点都要填充当然不是。一个在绝大多数样本里都缺失的位点填充出来的结果根本没有可靠性可言强行填进去只会给下游分析增加噪音。关于缺失率阈值我自己的经验是分两个维度看位点缺失率和样本缺失率。位点层面如果某个位点在超过10%的样本中缺失我会考虑在填充前剔除。样本层面如果某个样本的整体缺失率超过10%说明这个样本的DNA质量或芯片杂交过程可能有问题这类样本的填充结果也不可靠建议单独处理或剔除。# 查看每个样本缺失率 plink --file pharp_input --missing --out sample_missing # 输出文件里 F_MISS 列就是每个样本的缺失率这里有个个人心得阈值不是死的要看你的样本量和目标位点密度。如果样本特别少比如就四五十头剔除太多位点会导致填充参考信息不足反而得不偿失。这种情况下可以适当放宽到15%优先保证参与填充的位点数量。如果你的样本量充足严一点没问题。5. 提交前的自检清单与常见报错排查5.1 提交前十分钟逐项核对的检查清单每次准备完文件我不会立刻跑PHARP而是先过一遍自检清单。这套清单是我被坑过无数次之后总结出来的现在分享给你。第一检查文件行数。PED文件的行数应该等于样本数MAP文件的行数应该等于位点数两个文件的行数应该和你的预期完全一致。这一步用wc -l就能完成一旦发现行数不对说明数据里有重复行或缺失行必须马上排查。第二检查PED文件前6列。随机挑几行人工看一眼FID、IID、父本、母本、性别、表型这几列是否都是合法值。性别列只能是0/1/2出现其他数字就是有问题。表型列如果是连续性状不能有字母混入。第三检查基因型字符集。合法的两字符基因型只能是A、T、G、C、0这五种字符的组合如果出现“N”“R”“M”这些模糊碱基代码你需要决定是剔除还是转成缺失。PHARP默认读不懂含N的基因型。第四检查MAP文件排序。PHARP虽然不一定强制要求按染色体和物理位置排序但乱序会让后续的LD计算和填充效率大打折扣。建议按染色体编号排序同一条染色体内部按物理位置升序排列。第五检查文件编码。所有输入文件必须是纯文本格式不要是带BOM的UTF-8更不要是UTF-16。因为带BOM的文件第一行会多出一个不可见字符某些解析器会把第一个样本ID读错。5.2 常见报错信息与对应处理方案PHARP的报错信息在不同版本里五花八门但根因其实就那么几类。我挑了四个我实际遇到频率最高的来拆解一下。如果你看到“line X has fewer columns than expected”意思是PED文件某一行列数不够。最常见的原因是这个样本在某几个位点上缺失但你的转换脚本用了短格式而不是双字符“00”补位。PED文件不管有没有缺失每个位点都必须输出两个字符缺失就是“00”不能省略。如果你看到“unknown allele code”意思是基因型里出现了不认识的字符。排查方法很简单把你的PED文件里所有等位基因字符取一个unique集合看看有没有意外字符混进去。我遇到过的最离谱的情况是Excel转CSV时把“A”自动加上了引号引到了下一行。如果你看到“sample ID duplicated”意思是有重复样本ID。这个通常不是格式问题而是原始数据里同一个样本被检测了两次或者合并时没去重。解决办法是在转换之前对样本ID做去重保留分型完整度高的那一条记录如果两条完整度差不多保留第一条。如果你看到“position out of boundary”或者出现大量警告日志但没有直接中断大概率是MAP文件的物理位置越界或者染色体命名和软件内置的参考信息对不上。这时候就要回头检查MAP文件里的染色体写法和位置是否与参考基因组版本一致。5.3 我运行PHARP后还会顺手做的四件事文件成功跑进PHARP不代表万事大吉。填充完成后我一般还会做四个快速检查确认输出结果是真的可用而不是“看起来能用”。第一件事检查输出文件的缺失率。填充之后缺失率应该是显著下降的如果某个样本或某个位点的缺失率几乎没变说明这个样本或位点在填充过程中被跳过了需要回到文件准备阶段排查原因。第二件事检查等位基因频率的变化。把填充前后的等位基因频率做一个散点图大部分位点应该落在对角线附近。如果有一批位点明显偏离对角线通常是链方向搞反了或者填充参考面板本身有问题。第三件事检查填充质量分数如果PHARP输出了DR2或类似指标。不同软件的填充质量分数定义不同但逻辑一样分数越低填充结果越不可靠。我通常会把低质量位点过滤掉再做下游分析。第四件事按家系结构做一次抽样验证。选几个有完整系谱记录的家系检查后代样本是否在目标位点上符合孟德尔遗传规律。如果在家系里出现大量不符合孟德尔的位点说明填充的误差比预期大这种数据最好别急着用于后续选择或研究。文件准备是个不起眼的环节但它决定了后面所有工作的地基牢不牢。我见过太多人把时间花在调算法参数上却忽略了最基础的输入质量控制。希望这篇内容能帮你把“文件准备”这几个字背后的坑都提前绕过去。下次做猪芯片数据的填充时你大概率会发现最难的那关其实在你打开PHARP之前就已经过了。