一、研究背景与整体思路Rice3K3000 份水稻基因组计划3K Rice Genomes Project是水稻遗传研究、群体遗传分析和基因组选择中非常常用的公开资源。本文基于服务器目录/share/x-tech/OGR-rice/3K下的实际数据按照四个步骤逐块探索先看清表型压缩包的内容和格式再确认基因型三件套的规模与样本 ID 体系然后重点检查每个性状的表型 ID 与 fam ID 对齐率最后可选地做一轮 MAF 谱与缺失率初筛。整套流程的目的不是单纯统计文件而是回答三个后续建模必须解决的问题表型能不能直接与基因型 join位点规模是否符合预期ID 映射是否存在丢失或错配风险。二、第一步表型 zip 包结构与格式观察进入数据目录逐个解压查看每个*.zip里的文件清单、行数和头部几行内容。cd /share/x-tech/OGR-rice/3K for z in *.zip; do echo $z for f in $(unzip -Z1 $z); do echo -- $f ($(unzip -p $z $f | wc -l) 行) unzip -p $z $f | head -3 done done 21 | less culm_length.zip -- culm_length.txt (1765 行) IRIS_313-9814 27 IRIS_313-10617 35 IRIS_313-11980 36 culm_number.zip -- culm_number.txt (1766 行) IRIS_313-10796 5 IRIS_313-10917 5 IRIS_313-11127 6 days_to_heading.zip -- days_to_heading_2018HN.txt (2799 行) IRIS_313-10349 95 IRIS_313-7885 114 IRIS_313-10228 112 flag_leaf_length.zip -- flag_leaf_length_2018HN.txt (2792 行) CX2 26.49 CX3 28.65 CX4 24.67 flag_leaf_width.zip -- flag_leaf_width_2018HN.txt (2792 行) CX2 1.48 CX3 1.59 CX4 1.27 grain_length_width_ratio.zip -- grain_length_width_ratio.txt (2013 行) IRIS_313-10000 2.7 IRIS_313-10001 4.0 IRIS_313-10007 2.7 grain_length.zip -- grain_length.txt (2013 行) IRIS_313-10000 8.2 IRIS_313-10001 9.6 IRIS_313-10007 9 grain_width.zip -- grain_width.txt (2013 行) IRIS_313-10000 3.1 IRIS_313-10001 2.4 IRIS_313-10007 3.3 leaf_rolling_index.zip -- leaf_rolling_index_2015SZ.txt (1138 行) IRIS_313-8113 0 IRIS_313-8065 0 IRIS_313-9996 19.5 ligule_length.zip -- ligule_length.txt (1744 行) IRIS_313-10301 4 IRIS_313-8339 5 IRIS_313-10092 5 panicle_length.zip -- panicle_length.txt (1763 行) IRIS_313-10092 13 IRIS_313-10617 13 IRIS_313-10084 13 panicle_number.zip -- panicle_number_2018HN.txt (2771 行) CX2 10.1 CX3 7.9 CX4 9 plant_height.zip -- plant_height_2018HN.txt (2666 行) CX2 90.5 CX3 126.8 CX4 80.1 seedling_height.zip -- seedling_height.txt (1744 行) IRIS_313-10988 12 IRIS_313-10956 12 IRIS_313-11586 13 thousand_grain_weight.zip -- thousand_grain_weight.txt (1847 行) B001 24.3 B002 23.3 B003 19.8这一步重点看四件事每个 zip 里是一个还是多个 txt例如 days_to_heading 可能存在 2018HN 等多个年份地点版本而本地整合时只取了其中一个要确认最终采用的版本是否明确。列结构一般是 ID 与值两列以制表符分隔。是否有表头这直接决定第三步要不要使用NR1跳过表头。缺失值如何表示是 NA、空字符串还是 -9后续清洗时必须统一处理。三、第二步基因型三件套的规模与 ID 体系基因型数据以 PLINK 二进制格式存在fam、bim、bed 三件套分别描述样本、位点和基因型矩阵。先看样本规模和 ID 前缀构成。# 样本数量 ID 前缀构成3K 里有 IRIS_313-xxx 和 B/CX 等多种编号体系 wc -l NB_final_snp.fam 3024 NB_final_snp.fam head -3 NB_final_snp.fam B001 B001 0 0 0 -9 B002 B002 0 0 0 -9 B003 B003 0 0 0 -9 awk {print $2} NB_final_snp.fam | sed s/[0-9_-].*// | sort | uniq -c 246 B 312 CX 2466 IRIS head -3 NB_final_snp.bim 1 1026 0 1026 G A 1 1033 0 1033 G A 1 1047 0 1047 T A 1 1026 0 1026 G A │ │ │ │ │ │ 列1 列2 列3 列4 列5 列6 ┌─────┬────────────────────────────┬────────────────────────────────┐ │ 列 │ 含义 │ 这行的值 │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 1 │ 染色体号 │ 1 号染色体 │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 2 │ SNP 的 ID(名字,任意字符串) │ 1026(建库时直接拿位置当名字) │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 3 │ 遗传距离(cM) │ 0 占位符,没填,分析用不到 │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 4 │ 物理位置(第几个 bp) │ 1026 │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 5 │ 等位基因 1 │ G │ ├─────┼────────────────────────────┼────────────────────────────────┤ │ 6 │ 等位基因 2 │ A │ └─────┴────────────────────────────┴────────────────────────────────┘ awk {c[$1]} END{for(k in c) print k, c[k]} NB_final_snp.bim | sort -k1,1n 1 3057565 #1 号染色体上有 305.7 万个 SNP 2 2548707 3 2460262 4 2936357 5 2274771 6 2468410 7 2363916 8 2459990 9 1898922 10 2002536 11 2683870 12 2479918 #12 行加起来应该正好等于总数 29,635,224 awk length($5)1 || length($6)1 {n} END{print indel/多碱基位点:, n0} NB_final_snp.bim indel/多碱基位点: 0 其中: length($5)1 || length($6)1 —— 对每一行问:第 5 列或第 6 列的等位基因,字符串长度是否超过 1 个字母 {n} —— 满足条件的行,计数器加 1 END{print ..., n0} —— 读完全部行后打印计数(n0 是个小技巧:如果一行都没匹配到,n 是空的,加 0 强制显示成数字 0); 输出 0 —— 2963 万行里,没有任何一行的等位基因超过 1 个字母。这一步的核对标准是fam 应为 3024 行bim 约 2960 万行与本地3k_for_gs.log的记录保持一致。ID 前缀分布尤其重要3K 数据中并存 IRIS_313-xxx 和 B/CX 等多种编号体系前缀的构成决定了表型 ID 能否直接 join 到 fam。四、第三步表型 ID 与 fam ID 的对齐率重点这是整个探索中最建议先跑的一块。先生成 fam 样本 ID 集合再对每个性状 zip 的第一个文件逐性状计算交集。FAMIDS$(mktemp); awk {print $2} NB_final_snp.fam | sort -u $FAMIDS for z in *.zip; do f$(unzip -Z1 $z | head -1) n$(unzip -p $z $f | awk NR1{print $1} | sort -u | wc -l) hit$(comm -12 (unzip -p $z $f | awk NR1{print $1} | sort -u) $FAMIDS | wc -l) echo $z: 表型ID $n 个, 与fam直接匹配 $hit done rm -f $FAMIDS culm_length.zip: 表型ID 1764 个, 与fam直接匹配 1764 culm_number.zip: 表型ID 1765 个, 与fam直接匹配 1765 days_to_heading.zip: 表型ID 2798 个, 与fam直接匹配 2798 flag_leaf_length.zip: 表型ID 2791 个, 与fam直接匹配 2791 flag_leaf_width.zip: 表型ID 2791 个, 与fam直接匹配 2791 grain_length_width_ratio.zip: 表型ID 2012 个, 与fam直接匹配 2012 grain_length.zip: 表型ID 2012 个, 与fam直接匹配 2012 grain_width.zip: 表型ID 2012 个, 与fam直接匹配 2012 leaf_rolling_index.zip: 表型ID 1134 个, 与fam直接匹配 991 ligule_length.zip: 表型ID 1743 个, 与fam直接匹配 1743 panicle_length.zip: 表型ID 1762 个, 与fam直接匹配 1762 panicle_number.zip: 表型ID 2770 个, 与fam直接匹配 2770 plant_height.zip: 表型ID 2665 个, 与fam直接匹配 2665 seedling_height.zip: 表型ID 1743 个, 与fam直接匹配 1743 thousand_grain_weight.zip: 表型ID 1846 个, 与fam直接匹配 18463k_phenotype_wide.csv(表型宽表) —— 把集群上那 15 个 zip 里的性状 txt 整合成的一张总表:每行一个样本、每列一个性状,3105 行 × 15 性状。sample_id plant_height grain_length thousand_grain_weight ...B001 102.0 8.4 24.1IRIS_313-7684 NaN 8.9 20.73k_for_gs.fam(训练集的样本名单) —— 训练用基因型三件套的行名文件,2221 行,每行一个进入训练的样本 ID。就是参加考试的学生名单。模型训练时,程序拿 fam 里的每个 ID 去宽表里找同名行,配对成基因型 → 表型的训练样本。这一步的观察要点本地3k_phenotype_wide.csv的 sample_id 是 B001 体系而3k_for_gs.fam是 IRIS_313- 体系两者之间经过一次3k_rice_accession_mapping.csv映射。如果某些性状的原始 ID 匹配率偏低那么表型与基因型的错配或丢样就是效果差的最直接候选原因。需要注意脚本里的NR1假设文件有表头应先用第一步的结果确认如果没有表头要去掉NR1。五、第四步可选 MAF 谱与缺失率初筛如果时间允许可以用 plink2 快速出一份频率与样本缺失率报告结果只写到/tmp不污染原始目录。第一个染色体举例如下plink2 --bfile NB_final_snp --chr 1 --to-bp 50000 --export A --out /tmp/peekcut -f1,2,7-12 /tmp/peek.raw | head -5--chr 1 只要 1 号染色体 (看 bim 第 1 列)--to-bp 50000 只要位置 ≤ 50000 bp 的 (看 bim 第 4 列)#读取NB_final_snpbed/.bim/.fam三件套计算每个 SNP 的等位基因频率计算每个样本的基因型缺失率(34 秒,3024 样本 × 2963 万位点) plink2 --bfile NB_final_snp --freq --missing sample-only --out /tmp/3k_qc --threads 16 PLINK v2.0.0-a.7LM 64-bit Intel (25 Aug 2025) cog-genomics.org/plink/2.0/ (C) 2005-2025 Shaun Purcell, Christopher Chang GNU General Public License v3 Logging to /tmp/3k_qc.log. Options in effect: --bfile NB_final_snp --freq --missing sample-only --out /tmp/3k_qc --threads 16 Start time: Thu Aug 27 17:17:56 2026 257218 MiB RAM detected, ~227650 available; reserving 128609 MiB for main workspace. Using up to 16 threads (change this with --threads). 3024 samples (0 females, 0 males, 3024 ambiguous; 3024 founders) loaded from NB_final_snp.fam. 29635224 variants loaded from NB_final_snp.bim. Note: No phenotype data present. Calculating sample missingness rates... done. Calculating allele frequencies... done. --freq: Allele frequencies (founders only) written to /tmp/3k_qc.afreq . --missing: Sample missing data report written to /tmp/3k_qc.smiss . End time: Thu Aug 27 17:18:30 2026 head -2 /tmp/3k_qc.afreq #CHROM ID REF ALT PROVISIONAL_REF? ALT_FREQS OBS_CT 1 1026 A G Y 0.000625391 3198 #查看一共多少列 head -1 /tmp/peek.raw | awk {print NF} 2882 # 按表头自动找 ALT_FREQS 列重算: awk NR1{for(i1;iNF;i) if($iALT_FREQS) ci; next} {m($c0.5?$c:1-$c); n; if(m0) mono; else if(m0.01) r1; else if(m0.05) r5; else common} END{printf 总位点 %d\n单态: %d (%.1f%%)\nMAF0.01: %d (%.1f%%)\n0.01-0.05: %d (%.1f%%)\nMAF0.05: %d (%.1f%%)\n, n, mono,100*mono/n, r1,100*r1/n, r5,100*r5/n, common,100*common/n} /tmp/3k_qc.afreq 总位点 29635224 单态: 0 (0.0%) MAF0.01: 19921753 (67.2%) 0.01-0.05: 3557543 (12.0%) MAF0.05: 6155928 (20.8%) # 分布概况 最差 10 个样本 awk NR1{print $NF} /tmp/3k_qc.smiss | sort -n | \ awk {a[NR]$1} END{print 样本缺失率 中位a[int(NR*.5)], 90%a[int(NR*.9)], 99%a[int(NR*.99)], 最大a[NR]} sort -k5,5gr /tmp/3k_qc.smiss | head -11 IRIS_313-8921 IRIS_313-8921 29294555 29635224 0.988505 B203 B203 8916374 29635224 0.300871 B181 B181 8419324 29635224 0.284099 CX9 CX9 8392578 29635224 0.283196 IRIS_313-8859 IRIS_313-8859 8379822 29635224 0.282766 IRIS_313-8457 IRIS_313-8457 8248512 29635224 0.278335 IRIS_313-10539 IRIS_313-10539 8231076 29635224 0.277746 IRIS_313-11168 IRIS_313-11168 8195858 29635224 0.276558 IRIS_313-10729 IRIS_313-10729 7952722 29635224 0.268354 IRIS_313-9575 IRIS_313-9575 7916370 29635224 0.267127 IRIS_313-11174 IRIS_313-11174 7866043 29635224 0.265429 # 缺失率10% 的样本数(这些样本如果在 2221 训练集里,就是嫌疑) awk NR1 $NF0.10 {n} END{print 缺失率10% 的样本:, n0} /tmp/3k_qc.smiss 缺失率10% 的样本: 2576- .afreq → 全集 MAF 谱(67% 罕见变异、616 万常见位点池)→ 结论:位点筛选合理;- .smiss → 逐样本缺失率 → 揪出幽灵样本 IRIS_313-8921(98.9%)和 42 个高缺失样本 → 发现 --mind 缺口。这里主要观察 MAF 分位数判断低 MAF 位点的占比。如果 5% 分位已经非常接近 0说明有大量极低频位点后续基因组选择或关联分析前可能需要先做 MAF 过滤。六、四步检查总结步骤检查内容关键判断① 表型 zip文件数、列结构、表头、缺失值编码是否多版本 txt、缺失值如何表示② 基因型三件套fam/bim/bed 规模、ID 前缀fam 是否 3024 行、bim 是否约 2960 万行③ ID 对齐率逐性状 zip 与 fam 的交集低匹配率即错配或丢样线索④ MAF 谱plink2 频率与缺失率低 MAF 占比、是否需要过滤如果第三步各性状对齐率整体较高可以放心地按统一的样本 ID 体系建设表型-基因型对应表若个别性状匹配率明显偏低则需要回到3k_rice_accession_mapping.csv检查映射是否完整并对该性状单独复核而不是直接带入建模流程。