做叶绿体基因组分析的人这两年越来越多，但很多朋友一上来就直接跑软件，跑到一半发现组装出来的基因组不对、注释结果IR区边界明显异常、建树支持率低得没法看，然后又回头重来。说实话，这套流程看着就那几个环节——质控、组装、注释、比较分析、建树，但每一步的坑都比想象中多。我前前后后跑过上百套叶绿体基因组，踩过的雷基本都踩遍了，这篇就把整套高级分析的流程、参数、避坑点一次说清楚。
这套流程能用、且适合复现的前提是：你手里有植物样本的二代测序数据（最好是全基因组鸟枪法测序，或者纯叶绿体富集测序也行），目标是拿到完整的环状叶绿体基因组，然后往下做注释、比较基因组学、重复序列分析、密码子偏好性分析和系统发育分析。无论你是研究系统发育、群体遗传，还是做DNA条形码筛选，这套流程都能覆盖。我会把每一环节“为什么这么做”“参数怎么定”“出问题了怎么排查”都写透。
1. 分析全貌与流程设计
1.1 一套完整的叶绿体基因组高级分析都包含什么
叶绿体基因组是植物分子系统学里最常用的标记之一，因为它母系遗传、单拷贝数量大、保守区与高变区交替分布，既有全基因组比较的价值，又能筛出可用于设计的条形码区段。一套完整的高级分析，我大概分成六个环节：
数据质控与预处理、基因组组装与环化验证、注释（含IR边界确认）、重复序列与密码子分析、选择压力检验、系统发育分析。有些追求深入的组还会加群体遗传多样性分析、ndh基因丢失检测、RNA编辑预测等。
这里我多说一句，很多教程把“组装”当成重头戏，实际上注释阶段花的时间一点不少。组装可能跑一晚上就出了结果，注释之后你会发现IR边界对不上、基因有假阳性、tRNA注释缺位，这些都需要手动校正。所以整体时间分配上，组装占两成，注释占三成，后面四成在比较分析和反复的合理性验证上。心里有数，就不会在某个环节上过度焦虑。
1.2 不同研究目的下的流程裁剪
不是所有研究都需要跑全套。如果你是做单物种的叶绿体基因组图谱报道，那核心在组装+注释+重复序列+IR边界分析，系统发育可以只做基础的最大似然树；如果你做的是属内系统发育，那比较基因组、选择压力、单倍型网络才是重点，IR边界的详细图可能反而不需要；如果做的是DNA条形码筛选，你得额外关注高变区间的筛选统计，Pi值和核苷酸多样性是必须输出的指标。
所以我的建议是，动手前先花半小时明确三个问题：物种数量是多少、研究目的是什么、后续投稿大概需要哪些图。这直接决定你选哪套工具链、要不要跑PAML这种耗时大户。否则一把梭跑完，结果发现参考文献里的分析你少了两个，补起来又要浪费几天。
1.3 环境准备与软件清单
先说硬件。叶绿体基因组不大（通常是120kb到180kb），单样本组装对内存要求不高，16G内存的机器足够跑GetOrganelle；但如果你要批量处理几十上百个样本的比对、建树，建议上服务器或者有32G以上内存的工作站。存储上，一个样本的中间文件加起来大约5到10G，批量处理要注意磁盘空间。
软件链我推荐这么一套，全部开源且经过大量文献检验：
数据质控用fastp或Trimmomatic，选fastp是因为速度快，一条命令同时完成过滤、去接头和质量校正。
组装首选GetOrganelle，备选NOVOPlasty。NOVOPlasty需要种子序列，适合近缘参考明确的场景；GetOrganelle不需要参考基因组，通过延伸种子序列来组装，对叶绿体这种多拷贝基因组非常友好。
注释用GeSeq做自动化批注，或者用CPGAVAS2.在线的GeSeq方便，但有时需要手动校正；本地的CPGAVAS2可以批量跑，但依赖关系有点老，装起来可能费点功夫。tRNA注释会用到tRNAscan-SE，GeSeq内部已集成。
重复序列分析用MISA（SSR）和REPuter（散在重复序列）。
密码子偏好用CodonW或自写脚本处理，但CodonW在老系统上编译可能有问题，我用Python脚本比较多。
选择压力分析用KaKs_Calculator或PAML的codeml分支模型。
多序列比对用MAFFT，比对后修剪用trimAl，保守区段提取可以用Gblock。建树用IQ-TREE或RAxML-NG，我平时用得最多的是IQ-TREE，因为自带ModelFinder模型选择，对非专门做系统发育方法研究的人更省心。
单倍型网络用DnaSP输出格式，再丢到PopART里画图。
2. 从测序数据到环状基因组：组装与验证要点
2.1 数据质控与参考序列选择
拿到测序数据第一件事就是质控。叶绿体基因组虽然在细胞内拷贝数高，但测序数据里依然混着大量的核基因组片段和线粒体片段，所以质控不是简单去掉低质量碱基就行，而是要尽力保留有效叶绿体reads的占比。
我用fastp的命令一般长这样：
fastp -i R1.fastq.gz -I R2.fastq.gz -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe -q 20 -u 30 -n 10 -l 100 \ --thread 8 --json sample_fastp.json --html sample_fastp.html参数里-q 20表示质量值低于20的碱基会被修剪，-u 30允许最多30%的碱基低质量（超过这条read直接丢弃），-l 100是过滤掉短于100bp的read。
加这一句--detect_adapter_for_pe很关键，现在很多测序平台的接头序列和标准接头不完全一致，让fastp自动检测比手动指定更稳。
质控之后，我习惯用K-mer分布或者直接比对到叶绿体参考基因组看一眼reads占比。比如用一个近缘物种的叶绿体完整基因组做bwa mem比对，统计比对率：
bwa index ref_chloroplast.fasta bwa mem -t 8 ref_chloroplast.fasta clean_R1.fastq.gz clean_R2.fastq.gz > aln.sam samtools view -bS aln.sam > aln.bam samtools flagstat aln.bam如果比对率低于5%，意味着叶绿体reads占比太低，组起来会比较费劲。这种情况下可以考虑用GetOrganelle的--reduce-reads-for-coverage参数处理，或者干脆重新评估测序数据有没有问题。
2.2 组装策略选型：参考引导还是从头组装
叶绿体基因组组装有两种主流策略。第一种是参考引导组装，代表工具是NOVOPlasty。这种方案需要你提供一个近缘物种的完整叶绿体基因组做种子，程序从种子序列出发向两端扩展，把reads慢慢延伸拼接起来。
NOVOPlasty适合单样本快速组装，给一段2kb到5kb的种子序列就能跑，速度非常快。但隐患在于，如果研究的物种和参考序列差异大（比如存在大片段倒位、IR区结构变异），组装结果容易被参考序列带偏。你最后得到的可能不是这个物种的真实基因组结构，而是“参考序列的变异版本”。这是很多人在比较基因组学阶段发现结构差异很大的根本原因——不是物种间差异真的那么大，而是组装策略引入的偏差。
第二种是GetOrganelle这种基于De Bruijn图的从头组装策略，不需要参考基因组。GetOrganelle的思路是先对reads做多个k-mer规模的组装，再筛选出与目标细胞器相关的连接分支进行延伸。因为叶绿体基因组在细胞内有很高的拷贝数，相应k-mer的覆盖深度会显著高于核基因组片段，所以可以基于深度把目标序列“捞”出来。
我日常用GetOrganelle比较多，命令示范：
get_organelle_from_reads.py -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ -o chloroplast_assembly \ -F embplant_pt \ -t 8 \ -k 21,45,65,85,105 \ -R 15 \ --max-ext 50-F embplant_pt是使用植物叶绿体的预设参数集，内部包含对植物叶绿体做了优化的种子库和过滤规则。-k指定多个k-mer大小，从小到大覆盖不同重复结构的组装需求。-R 15表示对延伸过程中的低覆盖分支做15次迭代修复，这个参数在遇到复杂重复区时很管用，我一般不会低于15。
这里有个经验：如果测序深度很高（比如超过200X的叶绿体reads），GetOrganelle的自动深度过滤可能会导致组装结果被截断，这时候可以加--coverage-cutoff手动指定深度阈值，或者提高-R次数。我遇到过几次低深度样本反而组装比高深度样本顺的情况，原因就是高深度样本的嵌合read干扰了图结构。
2.3 组装结果的验证与环化
GetOrganelle跑完，会输出complete_graph和assembly_graph等fastg格式文件。不要直接拿走，先用Bandage打开看一下图结构。这个步骤特别重要，它让你直观看到基因组是单一环状结构还是断成了几截。
一个完整的叶绿体基因组应该呈现一个典型的“四部分结构”的环状：大单拷贝区（LSC）、小单拷贝区（SSC）、两个反向重复区（IRa和IRb）。在Bandage里，你会看到一个大的主环，IR区表现为一对对称的分支结构。如果看到的是开放线状、多个互不相连的contig，或者有可疑的额外环，说明组装不完整或混入了线粒体片段。
确认主结构之后，需要提取完整基因组序列。GetOrganelle会自动识别环化的起点，但有时起点位置在基因中间，影响后续注释起点的一致性。我的习惯是，先把环状序列导出，然后人工旋转，让序列起点落在trnH-psbA基因间区附近。这是叶绿体基因组比较中的一个默认习惯，NCBI上大多数叶绿体基因组序列的起点都在这个区域。旋转操作用一个小Python脚本就能实现：找到序列里trnH-psbA区段的坐标，将起点调整到该区域上游20bp左右的位置。
环化验证上还可以做一步：将组装结果与原始reads重新比对，检查覆盖率是否均匀。我一般用bwa加上samtools depth统计每个位点的覆盖度，然后看有没有明显的断崖式下降。如果某一区域覆盖度骤降，往往是组装错误或者该区域存在PCR扩增偏好。
3. 基因组注释：几何结构、IR边界与基因结构校正
3.1 注释工具与流程
组装完成后进入注释，这也是新手最容易“跑完就发”然后被审稿人骂的部分。自动化注释推荐用GeSeq平台，网址是chlorobox.mpimp-golm.mpg.de，在线提交序列之后它会整合ChloroplastDB和RefSeq数据库进行基因预测，同时调用tRNAscan-SE注释tRNA。
GeSeq的操作界面不算复杂，需要注意几个选项：参考数据库选Chloroplast_refseq，基因搜索模式选BLAT加HMMER联合，ID覆盖率阈值建议设到70到80。如果采用MUSCLE的参考比对模式，结果里会标注每个基因的完整度，方便判断伪基因。
在线注释完之后，强烈建议下载GFF3格式的结果，然后在本地用Plann或Chloe再做一轮交叉验证。因为GeSeq偶尔会在IR边界附近的基因上犯错，比如把ndhF跨IR边界的部分注释错、把ycf1这种高变基因漏掉或错位。交叉验证不是多此一举，我在十几个样本里发现过GeSeq把ycf1的长度预测差了好几百bp的现象。
本地另一款好用工具是CPGAVAS2，适合批量处理。CPGAVAS2的输出包含蛋白编码基因、rRNA、tRNA以及重复序列信息，还支持这种图片输出：IR区边界缩略图、基因组环形图。但CPGAVAS2依赖的BLAST和HMMER版本较老，在新系统下安装可能报错，编译时需要有点耐心。
3.2 IR区边界的收缩与扩张分析
IR区收缩扩张是叶绿体基因组比较分析里几乎必做的内容。判断IR边界说白了就是精确找出LSC与IR、IR与SSC之间的四个衔接点，通常基因标记分别是rps19基因（LSC/IRb边界）、ndhF基因（IRb/SSC边界）、ycf1或rps15（SSC/IRa边界）、trnH（IRa/LSC边界）。
拿到四个衔接点之后，建议用IRscope这个在线工具生成IR边界比较图，或者自己写脚本绘制。IRscope使用很简单，上传你要比较的多个物种的GenBank文件即可一次出图。
不过我可以很肯定地说，自动注释的IR边界位置经常会偏移，尤其是SSC/IR边界区域。这里的基因（ndhF和ycf1）往往部分跨越边界，软件如果对基因结构预估不准，边界就能差出几十甚至几百bp。
实操中的校正方法：把注释输出的GFF文件导入Geneious或IGV，与参考物种的注释结果对齐，人工检查ndhF、ycf1的终止密码子和IR区起始位置。更稳的办法是用BLAST找到IR区序列与LSC、SSC区序列的连接点，IR区在“反向重复”特性上可以通过比对自身序列之间的正反链相似区域来自动定位。我写过一个小流程：将完整基因组以1bp滑窗切成多段，检测两段之间是否呈反向互补关系，从而精确定位IR的边界。这个方法比盯着注释文件猜要可靠得多。
3.3 手动校正与可视化
注释结果手动校正阶段，最常用的工具是Apollo或者Geneious。如果你没有商业授权，用UGENE免费版也行，它能加载GFF和FASTA，直接查看基因与序列比对结果。
校正的标准动作是：先看蛋白编码基因有没有提前出现终止密码子、有没有跨越IR边界却不合理的情况；再看tRNA注释，tRNAscan-SE输出里的score低于50的基本都是可疑注释，要人工比对近缘物种的tRNA序列确认；最后检查rRNA，rRNA基因在叶绿体基因组里高度保守，如果注释结果里16S或23S rRNA长度明显异常，基本就是注释错误。
手工校正这一步很耗时，但回报也很大。别嫌麻烦，在叶绿体基因组相关期刊投稿时，审稿人非常看重IR边界和基因结构的准确性，我遇到过直接被要求提供手动校正截图的情况。
校正完成后，用OGDRAW（OrganellarGenomeDRAW）生成环形基因组图。OGDRAW输入GenBank格式，可以自定义颜色、基因分组、GC含量显示。它输出的是svg格式，可以直接用AI或者Inkscape微调，投期刊时基本够用。
4. 比较基因组学：重复序列、密码子偏好与选择压力
4.1 SSR与散在重复序列分析
叶绿体基因组里的SSR（简单重复序列）是群体遗传学里重要的分子标记来源，散在重复序列则与基因组重排、结构变异密切相关。
SSR检测我用MISA（MicroSAtellite identification tool），命令格式很简单：
misa sample.fastaMISA默认的搜寻规则是：单核苷酸重复至少10次，二核苷酸重复至少6次，三核苷酸重复至少5次，四、五、六核苷酸重复至少3次。这个默认规则我基本不调，因为如果放低阈值会检测出大量无意义微卫星。
统计输出结果时，除了SSR数量，还要区分完全型、复合型和间隔型。复合型SSR在叶绿体基因组里比例不低，很多时候可以作为高多态标记。另外，A/T碱基重复在叶绿体基因组SSR中占绝对主导，这是正常的，不用惊讶。如果结果里出现大量C/G重复，反而要怀疑是不是序列本身组装错误或者注释对象不是叶绿体。
散在重复序列检测用REPuter，它能识别正向、反向、回文和互补四种重复类型。命令：
reputer -h 3 -p -l 30 -e 3 sample.fasta-h 3表示汉明距离为3（允许3个错配），-p是显示正向重复，-l 30是最小重复长度30bp，-e 3是最大E-value阈值。我一般把最小重复长度设置在30bp，太短检测出来的大多是随机分布序列，没有实际意义。
这里分享一个常见误区：很多人在统计重复序列时只数数量，不区分长度类别。但审稿人经常要求提供重复序列的长度分布区间（如30-39bp、40-49bp、>50bp），因为长重复序列在基因组重排中作用更大。建议在REPuter输出结果里加一列长度，做后续柱状图或堆叠图时直接用。
4.2 密码子使用偏好性
密码子使用偏好性反映的是突变压力和选择压力的平衡结果，叶绿体基因组里一般会分析RSCU（相对同义密码子使用度）、GC含量（全序列、CDS、三位密码子位置）和有效密码子数（ENC）。
RSCU的计算逻辑并不复杂：某个密码子的RSCU值等于它的实际出现频率除以同义密码子的平均期望频率。RSCU大于1表示该密码子使用偏多，小于1表示使用偏少。
现在跑这个分析最省事的方式是直接用Python脚本解析GenBank或者FASTA+表格式CDS信息。我提供一个简化思路：把注释得到的CDS序列逐一提取，合并成一个长FASTA，然后用CodonW跑，或者用biopython自带的CodonUsage类：
from Bio.Seq import Seq from Bio.Data import CodonTable codon_table = CodonTable.unambiguous_dna_by_name["Standard"] count = {} for seq in coding_sequences: for i in range(0, len(seq) - 2, 3): codon = seq[i:i+3] count[codon] = count.get(codon, 0) + 1然后按各同义密码子族计算RSCU。如果你不想写代码，用dambe软件也行，它提供窗口化的密码子偏好分析。
密码子分析结果里有几个点需要留意。第一，叶绿体基因组一般AT含量高，所以第三位密码子偏好A/T结尾，这是普遍现象。第二，如果某个样本的GC3值异常高，可能混入了核基因序列。第三，在做物种间密码子偏好比较时，必须统一使用相同的注释版本，因为注释基因数量不同会直接影响统计结果。
4.3 核苷酸多态性与选择压力检验
选择压力分析是叶绿体基因组比较中的核心环节，尤其是当研究对象包含多个物种或者多个个体时。最常用的是Ka/Ks比值（也叫dN/dS），即非同义替换率与同义替换率的比值。Ka/Ks > 1 表示正选择，= 1 表示中性进化，< 1 表示纯化选择。
叶绿体基因组蛋白编码基因的Ka/Ks大多远小于1，因为叶绿体基因受强纯化选择。如果你算出来某个管家基因（比如matK、rbcL）的Ka/Ks大于1，基本可以断定是比对错位或者注释基因结构有误，要回去查序列。
常用工具有KaKs_Calculator和PAML。KaKs_Calculator适合批量计算基因对之间的Ka/Ks，推荐用MA方法或YN方法：
KaKs_Calculator -i input.axt -o output.kaks -m MA这里的input.axt是两两物种间同源基因的密码子比对结果。构建axt文件的常规步骤是：先用MAFFT对CDS核苷酸序列做比对，再翻译成氨基酸比对来校正（用PAL2NAL这个在线工具或者本地脚本实现），确保比对过程中不引入移码。
PAML的codeml适合做更复杂的假设检验，比如检测某个分支是否发生了加速进化。叶绿体基因组研究中，用branch-site model检测特定谱系的正选择基因比较常见。PAML的配置文件写起来容易劝退新手，但核心参数就几个：seqtype = 1（密码子模型）、model = 2（支模型）、NSsites = 2（正选择位点模型）。输出文件中有BEB分析得到的位点后验概率，后验概率大于0.95的位点可以报告为候选正选择位点。
做选择压力分析时，基因数量多的情况下建议写一个循环批处理，把每个基因的比对结果单独建目录跑，避免输出文件互相覆盖。还有一个经验是：对长度差异大的基因（比如ycf1在属内差异可能很大）要格外小心，建议先用氨基酸比对严格校正，极端情况下直接剔除该基因不纳入统计，否则会得出不可靠的正选择信号。
5. 系统发育分析与群体遗传信号挖掘
5.1 序列比对与保守区段提取
系统发育分析的第一步是对所有物种的叶绿体基因组做多序列比对。这里有两个选择：全基因组比对和特定基因/片段比对。如果物种间结构高度保守，全基因组比对没有问题；但如果研究样本之间有IR区收缩扩张或者倒位，全基因组比对会有大量错配区域，反而影响建树。
推荐流程是：先对完整叶绿体基因组使用MAFFT的E-INS-i策略比对（适合存在大片段结构差异的序列）：
mafft --auto --maxiterate 1000 chloroplast_genomes.fasta > chloroplast_aln.fasta然后必须进行比对优化和修剪。叶绿体基因组比对后插入缺失区域会产生大片gap。这些区域有很多是比对噪声而非真实的系统发育信号。用trimAl自动化修剪：
trimal -in chloroplast_aln.fasta -out chloroplast_aln_trimmed.fasta \ -automated1-automated1是启发式方法，自动决定最优修剪参数。对高度可变区存在较多的情况，可以改用gt 0.8参数表示保留在80%以上序列中出现的列，会更保守。
如果你只是做基于特定DNA片段的系统发育分析，比如只分析matK+rbcL+trnH-psbA的组合，建议先分别比对再合并。比对时需要注意各个片段的方向一致性，所有序列都要调整成相同方向后再做合并。合并后要检查每个样本的片段长度，如果某个样本某个片段缺失，在矩阵里用-补全，并在后续建树中明确这一点。
5.2 建树策略与模型选择
叶绿体基因组系统发育的主流建树方法是最大似然法和贝叶斯法。最大似然法推荐IQ-TREE，贝叶斯法用MrBayes或BEAST。日常绝大多数场景下，IQ-TREE就够用，命令：
iqtree2 -s chloroplast_aln_trimmed.fasta \ -m MFP \ -B 1000 \ -alrt 1000 \ -T AUTO-m MFP让IQ-TREE自动跑ModelFinder选择最优替代模型，不用自己纠结GTR还是TVM；-B 1000做1000次超快自举（UFBoot2），-alrt 1000做1000次SH-aLRT检验。
UFBoot2和SH-aLRT的支持率参考习惯是：SH-aLRT >= 80% 且 UFBoot >= 95% 认为支持率较高。很多文献只报bootstrap值，如果你用了IQ-TREE，建议把两个支持率都报上去，审稿人会更满意。
如果样本量小（比如少于20个物种）且研究的是近缘物种，可以考虑用*BEAST或者BEAST 2做物种树联合估计，但它对数据的要求高很多，叶绿体基因组单亲遗传的特点也会影响估计结果。我的个人建议是，除非做分化时间估算，否则跑最大似然树就够了，把主要精力放在拓扑结构解读和分支支持率检验上。
建树时间和数据量、模型复杂度直接相关。叶绿体全基因组比对矩阵一般有12万到18万个位点，几十个样本用IQ-TREE跑，单线程可能要几个小时。-T AUTO会自动选择线程并估计最优运行时间配置。实测下来，40个样本的全叶绿体基因组矩阵，8线程大概一晚上能跑完千次自举，不用焦虑。
5.3 单倍型网络与遗传多样性
如果研究对象是同一物种或近缘种的多个个体，单倍型网络图比系统发育树更直观。构建单倍型网络的标准步骤：先用DnaSP读取比对好的序列，识别所有变异位点，划分单倍型，输出单倍型序列文件。然后在PopART里选择TCS网络或者Median-Joining网络，直接导入DnaSP的输出文件就能画图。
DnaSP这个软件界面看起来有点古旧，但功能很扎实。操作流程是：File → Open Data File载入比对序列，然后点“Haplotypes”按钮，它会自动列出单倍型及序列所属个体。再把“Haplotype Data File”导出，用PopART打开。
群体遗传分析的另一个核心指标是核苷酸多样性Pi和单倍型多样性Hd。这些在DnaSP里也能直接算。需要注意的是，如果样本来自不同的地理居群，最好根据居群分组计算组内和组间多样性，并在结果里做AMOVA分析，看变异主要发生在居群内还是居群间。
做群体遗传时，叶绿体基因组与SSR数据结合起来用效果会比单独任何一种都好。叶绿体基因组提供的是高分辨率单倍型信息，SSR提供的是多态性位点分布，两者互补。我在做药用植物群体研究时经常看到，叶绿体基因组分型区分度高但变异位点少，SSR则相反。放一起解读会让结论立体很多。
6. 常见问题与排查技巧实录
6.1 组装失败与结果异常的排查
组装没跑出完整环状结构，是最常见的问题。可以按顺序排查：先看测序数据量，看看测序深度够不够。如果叶绿体reads占比低于5%，组装结果会有很多缺口。这时候可以考虑提高测序量，或者用靶向富集方法单独测叶绿体。如果占比尚可但组装碎片化，尝试把GetOrganelle的-R参数提高到30或50，-k列表里加入更大的k-mer。
还有一个容易忽略的点：如果你的物种和数据库里的种子库差异很大（比如移栽后高度杂交的园艺物种），GetOrganelle可能因为种子匹配失败而无法启动。这种情况下可以手动提供一段保守序列做种子，比如用近缘物种的完整基因组拆成小片段，放到--seed里喂给它。
组装结果里如果出现两个长度相近的环，其中一个是线粒体基因组的可能性很高。线粒体基因组的重复序列和叶绿体有部分相似，容易被误捞。判断方法是：对两个环分别注释，看基因类型。叶绿体基因组含有rbcL、matK、ndh家族等标志基因，线粒体则以atp、cox、nad基因群为主，一眼就能区分。
6.2 IR区注释错乱的典型现象和处理
IR区错乱最常见的表现是：注释结果中IRa和IRb的长度完全一样但基因组成不镜像对称，或者ndhF/ycf1的坐标跨越了IR边界但BLAST比对显示同源区域并不完整。
解决思路很简单：通过序列自身反向比对来精确定位IR边界。具体来说，IR区两段是反向互补的，把一个假设的边界位置分别在基因组上从两边取等长片段进行反向比对，比对分最高的边界位置就是真实边界。这个小脚本拆分出来不到50行，但比人工盯着看强太多了。
还有一种情况是组装阶段把IR区组装成单拷贝，这样整体基因组长度会比预期短，注释也会跟着错。遇到这种结构，建议回头检查Bandage图里IR区对应分支的覆盖度，如果IR分支覆盖度接近其他区域的两倍，说明组装没问题，是注释阶段没有正确识别重复结构，你应该把基因组重新环化后再注释。
6.3 建树支持率低的排查思路
系统发育树跑出来，很多分支的支持率都很低，新手容易慌张。第一步先看数据量：如果你只用了两三个片段，总长度不到3000bp，那支持率低是正常的，要么增加片段，要么用全叶绿体基因组。第二步看比对质量：重新审视比对中的移动gap区域，这些区域可能导致系统发育信号被噪声淹没，用trimAl更严格地修剪后再建树。第三步看物种采样：如果有远缘物种混入，长枝吸引效应会让部分节点支持率异常高或异常低，建议先用单基因树检查是否有样本的长枝异常。
支持率低不等于结果不能用，如果研究的是近缘物种，共享的稀有变异少，支持率天然不高。这个时候在文章里如实报告，并结合形态学或其他证据讨论，反而比强行解读更有说服力。
6.4 流程加速与批量化处理的建议
最后说一点批处理效率的问题。如果你要分析十个以上的样本，手工跑GeSeq再下载结果会非常痛苦。我建议把整个分析流水线脚本化。组装和质控可以用automated的snakemake或nextflow管理，注释则可以考虑在同一批中统一跑CPGAVAS2本地版。
批处理有一点务必注意：生成环形图或IRscope图时，要统一所有样本的序列起点方向。如果把起点不一致的样本放在一起比较，IR边界图看起来会非常混乱，而且很难在后期统一修改。我的做法是，在组装完成后立刻统一调整序列起点到trnH-psbA区域，之后再进入所有下游步骤。
算力方面，如果样本量几十个上百个，建议多用并行而不是依赖单机。GetOrganelle和IQ-TREE都支持多线程，把-t参数拉满或者用-T AUTO。BLAST搜索类的步骤尽量用--num_threads开启并行。云服务器按小时计费时，选择带AVX指令集的CPU对BLAST这类计算有明显加速效果。
在我个人试过的所有工具里，这个流程里最值得反复推敲的永远是注释阶段。很多分析结果出问题，追根溯源都是组装序列本身有小问题，或者注释基因结构不对。如果你能做到“组装完先看图，注释完先查IR，建树完先看支持率分布”这三步，整套流程会顺畅非常多。跑完一两个项目的经验，基本就能覆盖后续80%的异常情况了。
