输入数据与依赖环境
1.1 输入文件
文件 | 作用 |
*.vcf.gz | 群体变异信息 |
*.fa | 参考基因组 |
*.gff | 基因注释信息 |
*.csv | 需要分析的基因对ID信息 |
VCF 文件 (Variant Call Format)：主要记录个体或群体中的基因变异信息（单核苷酸多态性（SNP）、插入/缺失（Indel）以及结构变异（SV）等）。
GFF文件（Generic Feature Format）：描述基因组特征的注释文件，记录了基因、转录本、外显子等元件在染色体上的位置信息（如基因 ID、起始/终止坐标、链的方向 (+/-) 等）。通过结合 GFF 的坐标信息与参考基因组（FASTA 文件），可以精准地提取出目标基因的完整序列。
选择*.csv而非*.xlxs形式：*.csv 是一种纯文本格式，仅存储原始数据值，*.xlxs 是一种复杂的二进制/XML 压缩包，包含了字体、颜色、单元格样式、公式和宏等大量冗余信息，对于脚本程序而言，解析纯文本比解压并读取复杂的*.xlxs 结构要快得多。
1.2 软件依赖
bcftools
vcftools
mafft
python >= 3.8
pyfaidx
pandasbcftools：构建Consensus Sequence
MAFFT：多序列比对工具
详细步骤
2.1 第一步：根据样本名自动区分雌雄样本
在VCF文件header中解析雌雄样本：
样本名以 F 开头为雌性，以 M 开头为雄性（排除 MF龙凤榧）
VCF 的第 10 列开始是样本
def parse_sex_samples(vcf_file):
with gzip.open(vcf_file, "rt") as f:
for line in f:
if line.startswith("#CHROM"):
header = line.strip().split("\t")
samples = header[9:]
break
female = [s for s in samples if s.startswith("F")]
male = [s for s in samples if s.startswith("M") and not s.startswith("MF")]
return female, male
female_samples, male_samples = parse_sex_samples(VCF_FILE)2.2 第二步：从 GFF 中提取基因坐标信息
*.csv文件中只包含了旁系同源基因对的ID和location位置，并不清楚起始/终止坐标、链的方向 (+/-) 等信息，此时，就需要根据基因对ID在GFF文件中提取基因序列的起始/终止坐标等信息。
在*.csv文件中，基因对名称使用的是PTG000013L.28，与GFF文件中，ID为ID=evm.TU.PTG000013L.28
#脚本策略：
- 过滤与定位 (过滤非基因行)
跳过注释：忽略以 # 开头的行
筛选特征：通过 if feature_type != "gene": continue，代码会跳过所有外显子（exon）、内含子（intron）和 mRNA，只处理 gene 这一行
坐标转换：GFF 文件是基于 1 的坐标系统，代码通过 int(start) - 1 将其转换为 Python 字符串切片使用的 0-based 坐标系统
- 属性解析 (拆解第 9 列)
GFF 的第 9 列通常是一长串用分号 ; 隔开的字符串（如 ID=gene1;Name=...）
代码将其拆分为字典 attr_dict，方便后续根据键（Key）直接获取值（Value）
- ID 提取
多重检索：按顺序检查 ID、Name、gene_id 这三个字段
正则匹配：使用 re.search(r"(PTG\d+L\.\d+)", ...) 提取 ID
PTG：必须以这三个字母开头
\d+：匹配后续的一串数字
L\.：匹配字母 L 和一个点号 .
\d+：匹配点号后的数字
def load_geneinfo(gff_file):
gene_coords = {}
with open(gff_file) as f:
for line in f:
if line.startswith("#"):
continue
fields = line.strip().split("\t")
if len(fields) < 9:
continue
chrom, source, feature_type, start, end, score, strand, phase, attributes = fields
if feature_type != "gene":
continue
start = int(start) - 1
end = int(end)
attr_dict = {}
for item in attributes.split(";"):
if "=" in item:
k, v = item.split("=", 1)
attr_dict[k] = v
gene_id = None
for key in ["ID", "Name", "gene_id"]:
if key in attr_dict:
m = re.search(r"(PTG\d+L\.\d+)", attr_dict[key])
if m:
gene_id = m.group(1)
break
if gene_id:
gene_coords[gene_id] = (chrom, start, end)
return gene_coords
gene_coords = load_geneinfo(GFF_FILE)
print(f"[INFO] Loaded {len(gene_coords)} genes from GFF")2.3 第三步：构建样本特异的全基因组 consensus
使用bcftools构建Consensus Sequence：通过分析比对文件 (BAM/CRAM) 生成变异 (VCF) 文件，然后根据变异信息或直接基于覆盖度生成包含所有样本（或指定样本）的单一共有序列。
输入：标准参考基因组和包含样本变异的 VCF 文件
处理逻辑：脚本遍历 VCF 中的记录，如果发现某个位置有变异（如 SNP 或 Indel），就把参考基因组上对应的碱基替换成该样本特有的碱基
输出：生成一个全新的 FASTA 文件，这个文件反映了该样本真实的 DNA 序列
cmd = f"{BCFTOOLS} consensus -s {sample} -f {FASTA_FILE} -o {out_fa} {VCF_FILE}"-s {sample}：只提取某一特定样本的变异信息
-f {FASTA_FILE}：指定原始的参考基因组
-o {out_fa}：指定输出文件名
{VCF_FILE}：变异源文件
#脚本策略：“先整体，后局部”
传统做法：每个 gene 都跑一次 bcftools consensus → 非常慢
（8000个基因对 × 60个样本 = 48万次查询）
优化策略：每个雌雄样本 → 1 条完整 genome consensus （后续直接切片，不再重复调用 bcftools）
优势：VCF只解析一次且避免重复的变异应用，I/O效率高更适合并行化处理
1.全基因组 Consensus：先为每个样本（如 Male 1, Female 2）生成一整套完整的染色体序列。这样做的好处是，后续无论我们要提取多少不同的基因对（8000多对基因对），都只需要从这几个已经生成好的“Consensus Sequence”里截取坐标即可，不需要反复去读取 VCF 文件，大大提高了效率。
2.避免过滤困扰：bcftools consensus 会自动应用 VCF 中已有的过滤结果。如果 VCF 中某个位点是缺失的（./.），它默认会保留参考基因组的碱基，除非额外指定了掩码参数
def generate_genome_consensus(sample_list, sex_tag):
consensus_files = {}
for sample in sample_list:
out_fa = OUT_DIR / f"{sample}_{sex_tag}.fa"
if out_fa.exists():
print(f"[INFO] {out_fa} already exists, skip")
consensus_files[sample] = out_fa
continue
cmd = f"{BCFTOOLS} consensus -s {sample} -f {FASTA_FILE} -o {out_fa} {VCF_FILE}" # bcftools consensus直接生成全基因组序列
subprocess.run(
cmd,
shell=True,
check=True,
stdout=subprocess.DEVNULL,
stderr=subprocess.DEVNULL #关闭 bcftools Warning
)
consensus_files[sample] = out_fa
return consensus_files
male_consensus = generate_genome_consensus(male_samples, "male")
female_consensus = generate_genome_consensus(female_samples, "female")2.4 第四步：基于 gene pair 提取序列片段
使用 pyfaidx 高效切片（Fasta 文件支持随机访问，不需要加载整个文件到内存）
gene_info = (
pd.read_csv(GPINFO_FILE)
.dropna(subset=["gene1", "gene2"])
)GPINFO_FILE：*.csv文件中每一行的一对需要比较的旁系同源基因对（gene1 与 gene2）
至少包含两列：gene1和gene2
在主循环中逐行读取：
for idx, row in gene_info.iterrows():
gene1, gene2 = row["gene1"], row["gene2"]
pair_id = f"{gene1}_{gene2}"2.4.1加载并缓存每个样本的全基因组 Fasta
fasta_cache = {}
for sample, path in male_consensus.items():
fasta_cache[f"male_{sample}"] = Fasta(str(path))
for sample, path in female_consensus.items():
fasta_cache[f"female_{sample}"] = Fasta(str(path))通过 pyfaidx 的 Fasta 接口，将每个样本的全基因组 consensus 序列加载为可随机访问对象，从而避免重复读取和高内存占用
Fasta() 是 pyfaidx 提供的随机访问接口
每个样本的 genome consensus 只加载一次
后续所有基因切片都直接复用
2.4.2根据 GFF 坐标提取基因序列
def extract_sequences_optimized(sample_name, gene, gene_coords, tag, out_handle):
chrom, start, end = gene_coords[gene] #直接从缓存的字典中取
fa = fasta_cache[sample_name]
seq = fa[chrom][start:end].seq
out_handle.write(f">{gene}_{tag}\n{seq}\n")gene_coords[gene]：从 GFF 中获取该基因在基因组上的坐标
fa[chrom][start:end]：直接从 fasta 中随机访问该区域
.seq：返回纯序列字符串
>{gene}_{tag}：为 MAFFT 提供清晰的样本标识
2.4.3对 gene pair 中的两个基因逐一提取（主循环）
with open(final_fa, "w") as out:
for sample in male_samples:
extract_sequences_optimized(f"male_{sample}", gene1, gene_coords, f"male_{sample}", out)
extract_sequences_optimized(f"male_{sample}", gene2, gene_coords, f"male_{sample}", out)
for sample in female_samples:
extract_sequences_optimized(f"female_{sample}", gene1, gene_coords, f"female_{sample}", out)
extract_sequences_optimized(f"female_{sample}", gene2, gene_coords, f"female_{sample}", out)对每个 gene pair，分别从雄性和雌性样本的 genome consensus 中提取对应的基因序列，并合并为一个 FASTA 文件用于后续多序列比对alignment
2.5 第五步：对每个 gene pair 进行多序列比对alignment
#为什么要 alignment？
序列长度不同、插入缺失（indel）、后续分析需要（计算FST（Fixation Index），π核苷酸多样性（nucleotide diversity）等等）
#调用MAFFT对每个gene pair进行比对，自动对齐
def run_mafft(in_fa, out_aln):
cmd = ( #自动检测序列反向互补，对齐时统一方向
f"mafft --auto "
f"--adjustdirectionaccurately "
f"--thread {MAFFT_THREADS} "
f"{in_fa} > {out_aln}"
)
subprocess.run(
cmd,
shell=True,
check=True,
stderr=subprocess.DEVNULL #关闭MAFFT Warning
)--auto：自动选择算法
--adjustdirectionaccurately：自动反向互补
--thread：加速（设置多线程）
3.异常处理与可更新续写
#如何写一个可恢复可更新后不从头开始跑的 pipeline？（脚本优化策略）
3.1 缺失基因的处理策略
由于*.csv文件中，有部分fragment gene（碎片化基因：一个完整的基因序列被切断，分布在了不同的 Scaffolds（支架）或 Contigs（连续序列）上）不在GFF 中，程序需要跳过并记录到一个missing_file.csv里。
# 记录缺失基因
missing_genes = []
#输出缺失基因信息到 CSV
if missing_genes:
missing_file = OUT_DIR / "missing_genes.csv"
with open(missing_file, "w", newline="") as f:
writer = csv.DictWriter(f, fieldnames=["pair", "missing_gene"])
writer.writeheader()
for row in missing_genes:
writer.writerow(row)
print(f"Missing genes recorded in {missing_file}")3.2断点续跑设计
原因：
（1）可以后续更新脚本加入如MK Test、Fst计算等新的模块后，之前已经跑完的步骤可以被跳过，不会从头开始重新跑一遍。
（2）可以避免跑脚本如果因为意外而中断，要重新从头开始跑（不过一般挂后台的话基本不会发生），因为数据量比较大，跑一次要花上好几天时间。
可以使已生成存在的 consensus 不再重复生成
已完成的 alignment 不再重复跑
1. FASTA 非空检查（防止“空文件假完成”）
#检查Fasta文件是否有序列
def is_fasta_nonempty(fa_file):
with open(fa_file) as f:
for line in f:
if line.startswith(">"):
return True
return False主循环中的使用：
if aln_fa.exists() and is_fasta_nonempty(aln_fa): #如果最终比对已存在，整个gene pair直接跳过
print(f"[INFO] Alignment exists, skip pair {pair_id}")
continue这一层断点续跑解决的问题是：
MAFFT 可能：被 kill、被 Ctrl+C、机器重启
结果文件存在，但内容是空的或损坏的
2. 三层断点续跑
（1）第一层：全基因组 consensus 的断点续跑
def generate_genome_consensus(sample_list, sex_tag):
consensus_files = {}
for sample in sample_list:
out_fa = OUT_DIR / f"{sample}_{sex_tag}.fa"
if out_fa.exists():
print(f"[INFO] {out_fa} already exists, skip")
consensus_files[sample] = out_fa
continue每个样本只生成一次 genome consensus
脚本中断后重新运行：已完成的样本直接跳过，不会重复跑 bcftools
（2）第二层：gene pair 的 alignment 断点续跑
if aln_fa.exists() and is_fasta_nonempty(aln_fa):
print(f"[INFO] Alignment exists, skip pair {pair_id}")
continue以 gene pair 为单位：已成功对齐的 pair 不会重新计算，支持从任意 pair 继续跑
（3）第三层：MAFFT 内部
def run_mafft(in_fa, out_aln):
if out_aln.exists():
print(f"[INFO] MAFFT result exists, skip: {out_aln.name}")
return4.输出结果的结构
示例目录结构：
result/
├── geneA_geneB/
│ ├── geneA_geneB.fasta
│ └── geneA_geneB.aln.fa
├── male.samples
├── female.samples
└── missing_genes.csv部分结果示例：
geneA_geneB.aln.fa为我们最终需要的文件，
里面包含：
（1）已经alignment的雄性geneA变异后序列和（2）其旁系同源基因geneB变异后序列，
（3）已经alignment的雌性geneA变异后序列和（4）其旁系同源基因geneB变异后序列。
部分结果示例：（30个样本雄性gene1+gene2序列→60条序列，30个样本雌性gene1+gene2序列→60条序列）