1. 为什么绵羊需要定制单细胞参考基因组第一次接触绵羊单细胞转录组分析时我也曾天真地以为直接用人或小鼠的参考基因组就能搞定。结果发现这个想法就像用汽车钥匙开防盗门——完全不匹配。单细胞测序对参考基因组的要求比传统转录组严格得多特别是对于绵羊这类非模式生物官方预制参考基因组往往不可用。这里有个关键区别常规转录组分析使用的参考基因组包含所有染色体序列而单细胞测序需要特别处理过的参考基因组。10X Genomics的cellranger流程要求参考基因组必须包含精确的基因边界注释特别是外显子区域。我在实际操作中就遇到过因为注释文件版本不对导致30%的reads无法比对的情况。Ensembl数据库提供了绵羊的参考基因组Ovis_aries_rambouillet.Oar_rambouillet_v1.0但直接下载的原始文件需要经过三步关键处理筛选仅保留蛋白质编码基因的GTF注释验证线粒体基因完整性重构适合STAR比对器的索引结构最坑的是不同版本的注释文件质量差异很大。有次我用了Ensembl release-103的*.chr.gtf文件结果发现竟然漏掉了全部线粒体基因注释导致后续细胞质控完全失效。后来换用非chr版本才解决问题这个坑我踩了整整两天才爬出来。2. 实战第一步获取正确的基因组与注释文件2.1 下载基因组FASTA文件从Ensembl下载绵羊基因组时会发现有几种不同版本dna.toplevel.fa包含所有contig序列dna.primary_assembly.fa仅包含主要染色体dna.nonchromosomal.fa额外序列对于单细胞分析首选primary_assembly版本。但绵羊的primary_assembly版本缺失时可以用toplevel替代。我通常用wget后台下载大文件wget -b -c http://ftp.ensembl.org/pub/release-103/fasta/ovis_aries_rambouillet/dna/Ovis_aries_rambouillet.Oar_rambouillet_v1.0.dna.toplevel.fa.gz gunzip Ovis_aries_rambouillet.Oar_rambouillet_v1.0.dna.toplevel.fa.gz下载后务必检查是否包含线粒体序列grep Ovis_aries_rambouillet.Oar_rambouillet_v1.0.dna.toplevel.fa | grep MT2.2 获取GTF注释文件的陷阱cellranger只接受GTF格式的注释文件而Ensembl提供多种GTF变体.gtf完整注释.chr.gtf仅染色体注释.abinitio.gtf预测基因这里有个大坑某些版本的.chr.gtf会漏掉线粒体基因我第一次分析时就中了招。建议同时下载两个版本比对wget http://ftp.ensembl.org/pub/release-103/gtf/ovis_aries_rambouillet/Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.chr.gtf.gz wget http://ftp.ensembl.org/pub/release-103/gtf/ovis_aries_rambouillet/Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.gtf.gz # 检查MT基因是否存在 awk -F \t {print $1} *.gtf | grep -i MT | sort | uniq3. 注释文件过滤的关键技巧3.1 为什么必须过滤GTF文件原始GTF包含各种非编码RNA、假基因等注释。这些冗余信息会导致reads被错误标记为多重比对multi-mapped降低UMI计数准确性增加计算资源消耗cellranger提供的mkgtf工具可以过滤保留protein_coding基因cellranger mkgtf Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.gtf \ Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.filtered.gtf \ --attributegene_biotype:protein_coding但要注意某些GTF文件可能缺少gene_biotype属性。这时需要手动检查GTF结构# 查看GTF包含哪些属性 head -n 100 *.gtf | awk -F \t {print $9} | tr ; \n | awk -F {print $1} | sort | uniq3.2 处理特殊基因案例绵羊线粒体基因命名可能与人/小鼠不同。例如人线粒体基因MT-ND1, MT-COX1绵羊可能是ND1, COX1建议预先提取线粒体基因列表grep -i MT Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.filtered.gtf | awk -F \t {print $9} | awk -F ; {print $1,$3} | sort | uniq4. 构建参考基因组的完整流程4.1 运行cellranger mkref准备好过滤后的GTF和FASTA后使用mkref命令构建参考基因组nohup cellranger mkref \ --genomeovis_aries \ --fastaOvis_aries_rambouillet.Oar_rambouillet_v1.0.dna.toplevel.fa \ --genesOvis_aries_rambouillet.Oar_rambouillet_v1.0.103.filtered.gtf \ --memgb32 \ mkref.log 关键参数说明--memgb根据基因组大小调整绵羊基因组约2.8Gb建议32GB内存--nthreads当前版本有bug无法指定线程数4.2 监控运行过程构建过程通常需要2-4小时可以通过日志监控进度tail -f mkref.log # 预期看到类似输出 # Apr 15 16:56:08 ..... finished successfully # Reference successfully created! 4.3 验证输出结构成功的参考基因组应包含以下目录结构ovis_aries/ ├── fasta/ │ ├── genome.fa │ └── genome.fa.fai ├── genes/ │ └── genes.gtf.gz ├── reference.json └── star/ # STAR索引文件特别要检查reference.json中的内容cat ovis_aries/reference.json | python -m json.tool # 确认fasta_hash和gtf_hash不为空5. 常见问题排查指南5.1 线粒体基因缺失问题症状下游分析时发现所有细胞的MT基因占比为0% 解决方法确认FASTA包含MT序列检查GTF是否有MT基因注释必要时手动添加缺失的MT基因注释5.2 版本兼容性问题不同cellranger版本对参考基因组要求不同v6.x 需要GTF包含gene_version和transcript_version属性v5.x 可以接受简化版GTF遇到版本问题时可以尝试# 为旧版GTF添加必要属性 awk -F \t BEGIN {OFS\t} $3gene {$9$9; gene_version \1\;} $3transcript {$9$9; transcript_version \1\;} {print} input.gtf output.gtf5.3 外源基因添加技巧如果需要添加报告基因如GFP需要获取基因序列FASTA创建对应的GTF条目合并到原有文件中示例添加GFP# 准备GFP序列 echo -e GFP\nATGAGTAAAGGAGAAGAACTTTTCACTGGAGTTGTCCCAATTCTTGTTGAATTAGATGGT GFP.fa # 创建GTF条目 echo -e GFP\tartificial\texon\t1\t717\t.\t\t.\tgene_id GFP; transcript_id GFP; gene_name GFP; gene_biotype protein_coding; GFP.gtf # 合并到参考基因组 cat Ovis_aries_rambouillet.Oar_rambouillet_v1.0.dna.toplevel.fa GFP.fa combined.fa cat Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.filtered.gtf GFP.gtf combined.gtf6. 参考基因组优化建议经过多次实践我总结了几个优化点内存分配对于大型基因组建议分配至少1.5倍基因组大小的内存。绵羊基因组约2.8Gb设置--memgb32比较稳妥并行处理虽然cellranger不能直接控制线程数但可以通过Linux环境变量间接优化export STAR_NUM_THREADS16 export STAR_SORT_THREADS8临时目录指定高速存储作为临时目录能显著提升速度export TMPDIR/path/to/fast/ssd版本选择Ensembl的release版本不是越新越好。建议先检查该物种的golden path版本查看版本更新日志确认没有重大注释变更最后提醒一点每次构建参考基因组后建议保存完整的日志文件和软件版本信息。我曾经因为忘记记录cellranger版本号导致三个月后无法复分析结果。现在我的习惯是创建一个version.txt文件{ echo cellranger $(cellranger --version); echo Build date: $(date); echo FASTA md5: $(md5sum *.fa); echo GTF md5: $(md5sum *.gtf); } version.txt