WGDI的简单使用及预处理扩展

WGDI (Whole-Genome Duplication Integrated analysis) is a Python-based command-line tool designed to simplify the analysis of whole-genome duplications (WGD) and cross-species genome alignments.

WGDI（全基因组重复整合分析）是一款基于 Python 的命令行工具，旨在简化全基因组重复（WGD）和跨物种基因组比对分析。

https://github.com/SunPengChuan/wgdi

物种的进化总是有延续性的——这一点在基因组与染色体组上，就体现为多物种之间的同源区段及其动态变化。它们有时会产生多拷贝，有时会随着时间而打散甚至消失，还有一部分又重新整合，形成一个全新的基因簇。这些动态往往反映了生物的进化历史，是确定生物分化时间、亲缘关系等的直接证据之一。

WGDI就是在全基因组范围内，分析物种内与物种间重复（多拷贝）与整合（重组成基因簇）的工具。它的功能比较丰富，例如绘制 dot-plot、获取共线性参数、计算 Ks 值，等等。不同功能需要不同的配置文件，但是前置要求是一样的。你需要的文件包括：

• blast结果（采用format 6。可以用 NCBI-BLAST+，也可以用针对大规模数据优化的diamond）；
• GFF与LEN文件（需要特定脚本实现，有格式要求）。

其他参数将在后面介绍。

安装方法

你可以用pypi（pip）或者conda安装它：

pip3 install wgdi              # pypi
conda install bioconda::wgdi   # conda

如果使用conda但还没用过bioconda，请先在终端里执行这些指令：

conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

如果经历这种错误：

envs/wgdi/lib/python3.11/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.                       
......

卸载 WGDI 后重新执行安装：

conda create -n wgdi python=3.12 biopython wgdi

问题即可解决。

预处理过程

这里的 GFF 和 LEN 文件都不是标准格式文件。WGDI 对这两个文件的格式要求如下：

• 对于 GFF：

  列数	列名	解释
  1	染色体	该基因所在染色体代号
  2	ID	该基因的标识符
  3	起点	基因头部在该染色体上的碱基位置
  4	终点	基因尾部在该染色体上的碱基位置
  5	方向	基因所在链的方向（+/-）
  6	顺序	该基因在该染色体上的顺序。从 1 开始
  7	原始信息	原始ID或附加信息，可无。不用于分析

• 对于LEN：

  列数	列名	解释
  1	染色体	该染色体代号
  2	长度	该染色体的碱基长度
  3	基因数	该染色体上的基因数量

怎么重构？我们来看GFF原始文件。它包含主条目与各子条目，例如，这是取自 Ensembl Rapid 上的其中的“两个”条目：

###
1	ensembl	ncRNA_gene	749	1456	.	+	.	ID=gene:ENSFBZG00000004380;biotype=lncRNA;gene_id=ENSFBZG00000004380;version=1
1	ensembl	lnc_RNA	749	1456	.	+	.	ID=transcript:ENSFBZT00000006931;Parent=gene:ENSFBZG00000004380;biotype=lncRNA;tag=Ensembl_canonical;transcript_id=ENSFBZT00000006931;version=1
1	ensembl	exon	749	868	.	+	.	Parent=transcript:ENSFBZT00000006931;constitutive=0;exon_id=ENSFBZE00000031713;rank=1;version=1
1	ensembl	exon	947	1456	.	+	.	Parent=transcript:ENSFBZT00000006931;constitutive=0;exon_id=ENSFBZE00000031716;rank=2;version=1
###
1	ensembl	gene	6948	8703	.	+	.	ID=gene:ENSFBZG00000000663;biotype=protein_coding;gene_id=ENSFBZG00000000663;version=1
1	ensembl	mRNA	6948	8703	.	+	.	ID=transcript:ENSFBZT00000001038;Parent=gene:ENSFBZG00000000663;biotype=protein_coding;tag=Ensembl_canonical;transcript_id=ENSFBZT00000001038;version=1
1	ensembl	five_prime_UTR	6948	7244	.	+	.	Parent=transcript:ENSFBZT00000001038
1	ensembl	exon	6948	7566	.	+	.	Parent=transcript:ENSFBZT00000001038;constitutive=0;exon_id=ENSFBZE00000005188;rank=1;version=1
1	ensembl	CDS	7245	7566	.	+	0	ID=CDS:ENSFBZP00000001038;Parent=transcript:ENSFBZT00000001038;protein_id=ENSFBZP00000001038;version=1
1	ensembl	CDS	7841	8439	.	+	2	ID=CDS:ENSFBZP00000001038;Parent=transcript:ENSFBZT00000001038;protein_id=ENSFBZP00000001038;version=1
1	ensembl	exon	7841	8703	.	+	.	Parent=transcript:ENSFBZT00000001038;constitutive=0;exon_id=ENSFBZE00000005190;rank=2;version=1
1	ensembl	three_prime_UTR	8440	8703	.	+	.	Parent=transcript:ENSFBZT00000001038
###

是的，这并非是“12个”条目，而是两个—— ensembl 把不同的主条目用 ### 行隔开。在这里，每个条目的第一行就是我们要找的“主条目”了。我们只需要gene就好，另一个是ncRNA_gene，它不编码蛋白质，所以在这里并没有用。这可以靠检索 GFF 的第三列来确定。

再来看每一列，它们代表的应该是

染色体	来源	类型	起点	终点	得分	方向	CDS相位	其他属性

而 WGDI 需要的格式是

染色体	基因代号	起点	终点	方向	在该染色体上的顺序

第1、3、4、5可通过 GFF 直接获得（原第1、4、5、7），基因代号可通过原 GFF 的其他属性（第8列）获得，而最后一列通过染色体列间接获得。

先解决前五列：

    # GFF3 解析
    gff_data = pd.read_csv(args.ingff3, sep="\t", header=None, comment="#", dtype={0: str, 3: int, 4: int, 8: str})
    gff_data.columns = ['chr', 'source', 'feature', 'start', 'end', 'score', 'strand', 'phase', 'attribute']

    # 重新提取数据：'feature' = type 参数的 dataframe
    gff_data = gff_data[gff_data['feature'] == args.type].loc[:, ['chr', 'attribute', 'start', 'end', 'strand']]

    # 取各 feature 的实际 args.type
    gff_data['attribute'] = gff_data['attribute'].str.extract(r'ID=([^;]+)')

    # 重命名
    gff_data.columns = ['chr', args.type, 'start', 'end', 'strand']

最后一列可以通过重做索引的方法实现：

    # 重做索引
    gff_data = gff_data.sort_values(by=['chr', 'start']).reset_index(drop=True)
    gff_data['order'] = gff_data.groupby('chr').cumcount() + 1

我们还可以借此机会生成LEN文件，它需要所有染色体长度与基因数信息。染色体长度可通过该染色体最后一个基因的终点间接确定（虽然不太准确，但在没有DNA测序的情况下可作为近似）；基因数通过第一列里有该染色体的行数来确定。

    # 统计各染色体核苷酸长度
    chr, lens = [], []
    for rec in SeqIO.parse(args.infna, 'fasta'):
        chr.append(rec.id)
        lens.append(len(rec.seq))
    chr2length = pd.DataFrame(data={'chr': chr, 'lens': lens})

    # 聚合方法，统计 'attribute' 条目数（count）
    chr2count = gff_data.groupby('chr').agg({args.type: 'count'}).reset_index()

    # 合并
    lens = pd.merge(chr2length, chr2count, on='chr')

综上，生成GFF和LEN的任务可以合并到一个脚本来完成（1_get_gff_len.py）。

---

我们需要的3个文件，现在已经搞定了2个。还差个blastp的结果——可是蛋白质组数据的问题，还没有解决！这样如何获得blast结果嘛！

有人会说，获取蛋白质数据，完全可以直接用 gffread 来实现嘛！这确实可行，但实际操作的时候，还是会差点事。

例如，每一个蛋白质条目变成了转录本，基因与蛋白质的对应关系消失了（特别常见于 ensembl 版的 GFF）：

> transcript:ENSWHXT00000003709

另一种，基因已知，但分成了若干转录本（一般课题组自己获得的就是这种）：

> g11.t1
> g11.t2

我们显然不希望发生这两种情况。好在这俩其实从根本上是同一个问题——如果一个基因条目关联到多个转录本，如何选择？我们更加倾向于选择更长的转录本。例如这种：

###
1	ensembl	gene	123687	127517	.	+	.	ID=gene:ENSWHXG00000001281...
1	ensembl	mRNA	123687	127124	.	+	.	ID=transcript:ENSWHXT00000002004...
1	ensembl	five_prime_UTR	123687	123766	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	exon	123687	124332	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	CDS	123767	124332	.	+	0	ID=CDS:ENSWHXP00000002004...
1	ensembl	exon	124433	124766	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	CDS	124433	124766	.	+	1	ID=CDS:ENSWHXP00000002004...
1	ensembl	exon	124974	125247	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	CDS	124974	125247	.	+	0	ID=CDS:ENSWHXP00000002004...
1	ensembl	CDS	125767	126125	.	+	2	ID=CDS:ENSWHXP00000002004...
1	ensembl	exon	125767	127124	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	three_prime_UTR	126126	127124	.	+	.	Parent=transcript:ENSWHXT00000002004...
1	ensembl	mRNA	123687	127517	.	+	.	ID=transcript:ENSWHXT00000002017...
1	ensembl	five_prime_UTR	123687	123766	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	exon	123687	124332	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	CDS	123767	124332	.	+	0	ID=CDS:ENSWHXP00000002017...
1	ensembl	exon	124433	124766	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	CDS	124433	124766	.	+	1	ID=CDS:ENSWHXP00000002017...
1	ensembl	exon	124974	125247	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	CDS	124974	125247	.	+	0	ID=CDS:ENSWHXP00000002017...
1	ensembl	exon	125767	125948	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	CDS	125767	125948	.	+	2	ID=CDS:ENSWHXP00000002017...
1	ensembl	CDS	126507	126518	.	+	0	ID=CDS:ENSWHXP00000002017...
1	ensembl	exon	126507	127517	.	+	.	Parent=transcript:ENSWHXT00000002017...
1	ensembl	three_prime_UTR	126519	127517	.	+	.	Parent=transcript:ENSWHXT00000002017...
###

好复杂！但稍微精简一下，结构就一目了然：

###
1	gene	123687	127517	.	+	gene:ENSWHXG00000001281
1	mRNA	123687	127124	.	+	transcript:ENSWHXT00000002004
1	CDS	123767	124332	.	+
1	CDS	124433	124766	.	+
1	CDS	124974	125247	.	+
1	CDS	125767	126125	.	+
1	mRNA	123687	127517	.	+	transcript:ENSWHXT00000002017
1	CDS	123767	124332	.	+
1	CDS	124433	124766	.	+
1	CDS	124974	125247	.	+
1	CDS	125767	125948	.	+
1	CDS	126507	126518	.	+
###

很明显，一个gene带了两个transcript！怎么做呢？将它们做差再求和，然后取CDS最长的那个版本。

transcript:ENSWHXT00000002004
124332-123767+1 = 566
124766-124433+1 = 334
125247-124974+1 = 274
126125-125767+1 = 359
sum = 1533

transcript:ENSWHXT00000002017
124332-123767+1 = 566
124766-124433+1 = 334
125247-124974+1 = 274
125948-125767+1 = 182
126518-126507+1 = 12
sum = 1368

此时ENSWHXG00000001281对应的“唯一”转录本为ENSWHXT00000002004。我们首先要记住这个对应关系，然后将这个转录本“抽”出来，最后还要把转录本的名字换成基因的名字。——这可能吗？

过程显然要比上一个脚本复杂一些，但不一定无法完成。整理一下思路，我们首先要找到每个基因当中 CDS 最长的那个转录本，然后用这个转录本从基因组中翻译出蛋白质序列，最后用基因的名字命名序列。

同样的，我们还是要请 Pandas 来帮我们组织一下数据。我们的流程是——

• 利用 mRNA 做第一次”遍历”, 建立基因 -> 转录本的关系（虽然看上去是“遍历”，但其实是每一行同时进行处理）。

      # 建立基因到转录本(mRNA)的关系
      print(f"\033[94mAnalysing\033[0m make dictionary for gene to transcript...")
      g2t_df = pd.DataFrame(columns=['gene', 'transcript', 'chromosome'])
      g2t_df['gene'] = mrna_data['attribute'].str.extract(r"Parent=([^;]+)")
      g2t_df['transcript'] = mrna_data['attribute'].str.extract(r"ID=([^;]+)")
      g2t_df['chromosome'] = mrna_data['chr']
    
  

• 利用 CDS 做第二次”遍历”, 获得转录本 CDS 长度。

      # 计算 CDS 区长度，以便选出 CDS 最长的转录本
      print(f"\033[94mAnalysing\033[0m calculating transcript CDS length...")
      cds_data = cds_data.loc[:, ['attribute', 'start', 'end']].copy()
      cds_data['transcript'] = cds_data['attribute'].str.extract(r"Parent=([^;]+)")
      cds_data['start'] = cds_data['start'].astype(int)
      cds_data['end'] = cds_data['end'].astype(int)
      cds_data['length'] = cds_data['end'] - cds_data['start'] + 1
      cds_length_df = cds_data.groupby('transcript')['length'].sum().reset_index()
    
  

• 通过公共项“转录本名称”合并1和2的结果，然后，通过选择 CDS 最长的转录本，形成转录本 -> 基因一一对应的关系（字典）。

      # 合并两个dataframe，然后找最长的转录本
      print(f"\033[94mAnalysing\033[0m search the longest CDS...")
      transcript_table = pd.merge(g2t_df, cds_length_df, on='transcript')
      transcript_table = transcript_table.loc[:, ['gene', 'transcript', 'length', 'chromosome']]
      transcript_table_filtered = transcript_table.loc[transcript_table.groupby('gene')['length'].idxmax()]
      transcript_table_filtered['gene_order'] = transcript_table_filtered.groupby('chromosome').cumcount() + 1
    
      # 收集结果（转录本到基因、染色体编号等的关系）
      print(f"\033[94mAnalysing\033[0m collate result...")
      print(transcript_table_filtered)
      trans2gene = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene']))
      trans2chr = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['chromosome']))
      trans2order = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene_order']))
    
  

• 最后，在gffread的结果里：

  1. 删掉字典里没有出现的转录本；
  2. 将转录本名字换成基因名字。

    # 对生成的序列进行筛选
    # 若序列label与transcript库匹配，则取用
    filtered_records_cds = [
        record for record in SeqIO.parse(args.outfna, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]
    filtered_records_prot = [
        record for record in SeqIO.parse(args.outfaa, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]

    # 将transcript名称改为基因名称
    for idx, record in enumerate(filtered_records_cds):
        record.description = f"index={idx+1},gene={trans2chr[record.id]}g{trans2order[record.id]}"
        record.id = trans2gene[record.id]
    SeqIO.write(filtered_records_cds, args.outfna, "fasta")

    for idx, record in enumerate(filtered_records_prot):
        record.description = f"index={idx+1},gene={trans2chr[record.id]}g{trans2order[record.id]}"
        record.id = trans2gene[record.id]
        record.seq = record.seq.replace(old='.', new="X")
    SeqIO.write(filtered_records_prot, args.outfaa, "fasta")
    

具体流程可以参考第二个文件 2_get_prot.py。

至此，大部分命令所需的全部文件已经集齐。

正式的数据处理过程

WGDI功能非常丰富，它可以通过不同类型的配置文件，实现不同的功能：

-d：dot-plot，做同源基因点图，可作为文献插图使用。 -icl：collinearity，获得共线性结果，本质是改进版的 ColinearScan。进一步分析时可以利用。 -ks：计算同源基因对的 Ka/Ks。 -bi：blockinfo，将共线性与Ks结果合并到一起。

这里我们只演示 -d 和 -icl 的用法。

# dotplot.conf
[dotplot]
blast = blast file
blast_reverse = false
gff1 =  gff1 file
gff2 =  gff2 file
lens1 = lens1 file
lens2 = lens2 file
genome1_name =  Genome1 name
genome2_name =  Genome2 name
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 10 
position = order 
blast_reverse = false 
ancestor_left = none 
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = savefile(.svg,.png,.pdf)

之后再 wgdi -d dotplot.conf。

等待出图。

然后，执行wgdi -icl coll.conf。

# coll.conf
[collinearity]
gff1 = gff1 file
gff2 = gff2 file
lens1 = lens1 file
lens2 = lens2 file
blast = blast file
blast_reverse = false
multiple = 1 # 用红点显示的同源基因的最佳数量
process = 8 # 使用进程数量
evalue = 1e-5
score = 100 # 和点图保持一致
grading = 50,40,25 # 红，蓝，灰的不同分值,按照分值会优先连红色的点，蓝色点次之
mg = 40,40 # 两个基因对被认为能连起来的最大距离（罚分）
pvalue = 0.2 # 评估共线性的指标，如果想全部提取可以选1。共线块的compactness and uniqueness，范围为0-1，较好的共线范围为0-0.2
repeat_number = 10 # 和点图保持一致 
positon = order # 唯一选项
savefile = collinearity file # 输出文件

结果大概是这个样子（MCScanX风格）：

# Alignment 1: score=401 pvalue=0.1706 N=15 ped_chr001&Chr01 plus
g1362 1242 Pfi01G025700 2570 1
g1363 1243 Pfi01G025830 2583 -1
g1385 1262 Pfi01G025960 2596 1
g1401 1276 Pfi01G026040 2604 -1
g1434 1306 Pfi01G026060 2606 1
g1435 1307 Pfi01G026080 2608 1
g1451 1322 Pfi01G026310 2631 1
g1452 1323 Pfi01G026320 2632 1
g1457 1328 Pfi01G026350 2635 1
g1460 1330 Pfi01G026380 2638 1
g1465 1334 Pfi01G026400 2640 1
g1466 1335 Pfi01G026440 2644 1
g1470 1339 Pfi01G026460 2646 1
g1485 1354 Pfi01G026700 2670 1
g1492 1360 Pfi01G026710 2671 1
# Alignment 2: score=383 pvalue=0.0002 N=8 ped_chr001&Chr01 plus
g634 586 Pfi01G036910 3691 1
g635 587 Pfi01G036920 3692 1
g636 588 Pfi01G036940 3694 1
g637 589 Pfi01G036960 3696 1
g638 590 Pfi01G036970 3697 1
g639 591 Pfi01G036980 3698 -1
g640 592 Pfi01G036990 3699 -1
g641 593 Pfi01G037010 3701 -1
......

一点点的Q&A

如果 Chr 顺序不太对，或者名称太长互相叠叠乐？

想修改染色体及其顺序，可以直接动 LEN 文件。

改名可以通过另一个脚本来实现。

"""
map文件的参考格式：
旧染色体名1\t新染色体名1
旧染色体名2\t新染色体名2
旧染色体名3\t新染色体名3
......

获得方法：执行 cat 'FNA文件路径' | grep '>' > map.txt，然后补上 Tab，删掉多余内容即可。
"""

import argparse
import re

from Bio import SeqIO

def switch_dict(map_path: str) -> dict:
    dictionary = {}
    with open(map_path, "r") as map:
        for line in map:
            line_ls = line.strip().split("\t")
            dictionary[line_ls[0]] = line_ls[1]
    return dictionary

def default_output_name(filename):
    name, ext = filename.rsplit('.', 1)
    return f"{name}_changed.{ext}"

def main():
    print()
    parser = argparse.ArgumentParser(description="Renaming the chromosome name in fna/gff3 file.")

    parser.add_argument('--infna', required=True, help='Input FNA file name.')
    parser.add_argument('--ingff', required=True, help='Input GFF3 file name.')
    parser.add_argument('--inmap', required=True, help='Old and new Chr name, in TSV format.')
    parser.add_argument('--outfna', required=False, help='Output FNA file name. Optional.')
    parser.add_argument('--outgff', required=False, help='Output GFF3 file name. Optional.')

    args = parser.parse_args()
    outfna = args.outfna if args.outfna else default_output_name(args.infna)
    outgff = args.outgff if args.outgff else default_output_name(args.ingff)

    d = switch_dict(args.inmap)
    fna_record = []
    gff_record = ""
    
    # 更改 FNA 里的条目
    for idx, rec in enumerate(SeqIO.parse(args.infna, "fasta")):
        # 改名；若无对应则保持原状
        try: rec.id = d[rec.id]
        except KeyError: rec.id = rec.id
        # 添加记录
        fna_record.append(rec)

    # 写入
    SeqIO.write(fna_record, outfna, "fasta")

    # 更改 GFF 里的条目
    pattern = re.compile("|".join(re.escape(k) for k in d.keys()))
    with open(args.ingff) as gff_original, open(outgff, "w") as gff_changed:
        for idx, rec in enumerate(gff_original):
            gff_new = pattern.sub(lambda m: d[m.group(0)], rec)
            gff_record = gff_record + gff_new
        gff_changed.write(gff_record)
    
    print("DONE!")

if __name__ == "__main__":
    main()

右侧和下方有很多 Chr 叠在一起，很不好看？

可以直接在 LEN 文件里，删掉基因数量非常少的，没有组装到染色体上的 Chr，然后再重新作图。

接下来该如何做？例如美化共线性结果？或者针对某个基因座进行更精细的共线性分析？

这部分推荐交给 jcvi，或者 TBtools-II，然后使用 -icl 的结果。不过——如果想作图，建议先改 LEN 文件，把散在的 Chr 删掉，再跑 WGDI，效果会更好。

改 e-value 或者得分阈值，对结果影响大吗？

e-value的影响不是很大。得分阈100可行，200差异不大，但500会破坏共线性关系。因此直接用默认参数就好。

能否只使用一个配置文件？

可以。

运行 wgdi -xxxx ? > total.conf（-xxxx = -d/-icl/…），即可生成单功能的配置模板，然后就可以根据自己的需求，补充缺失的参数。继续生成时，把 > 改成 >> 即可。

但这就有了一个问题：同一个文件下，不同 section（一个 [] 即一个 section）下可能使用同一个参数，例如 gff len 之类。能否只动一处，就能改变全局的同名参数？

关键在于 .ini 文件写法。可以先设一个 [universal] section，把公共参数写进去，其他 section 调用时，则使用 ${universal:argument}：

[universal]
blast = BLAST_RES/Pedel2Pvulg.f6
gff1 =  GFF/Pedel_Simp.gff
gff2 =  GFF/Pvulg_Simp.gff
lens1 = LEN/Pedel.len
lens2 = LEN/Pvulg.len
......

[dotplot]
blast = ${universal:blast}
gff1 =  ${universal:gff1}
gff2 =  ${universal:gff2}
lens1 = ${universal:lens1}
lens2 = ${universal:lens2}

但代价是要在源码上做改动。在 WGDI 源码中，需要对 base.py 的第 29 行进行修改：

# envs/py312/lib/python3.12/site-packages/wgdi/base.py
- conf = configparser.ConfigParser()
+ conf = configparser.ConfigParser(interpolation=configparser.ExtendedInterpolation())

这能让 ConfigParser 同时兼容这种简并写法。

附录

处理脚本源码

附上这两个脚本的源代码，各位可以随意取用，或者进一步改进。这些代码是我在做本科毕设的共线性分析时写的。

"""
WGDI Pre-Process Step-1

以下预处理代码改进自原 WGDI 项目。
使用前请确认是否安装了 Pandas 和 BioPython。

GitHub主页: https://github.com/SunPengChuan/wgdi?tab=readme-ov-file
相关论文: https://doi.org/10.1016/j.molp.2022.10.018
脚本原始来源：https://lxz9.com/2021/09/06/wgdi/#%E8%84%9A%E6%9C%AC
"""

import argparse
import pandas as pd
from Bio import SeqIO

def main():
    print()
    print("WGDI Pre-Process \033[94mStep-1\033[0m")

    # 用法说明与参数解析
    parser = argparse.ArgumentParser(description="WGDI Pre-Process Step-1: Get GFF and LEN Files.")
    parser.add_argument('--infna', required=True, help='Input FNA file name.')
    parser.add_argument('--ingff3', required=True, help='Input GFF3 file name.')
    parser.add_argument('--outgff', required=True, help='Output simplified GFF file name.')
    parser.add_argument('--outlen', required=True, help='Output simplified LEN file name.')
    parser.add_argument('--type', required=True, help='Entries type.')
    args = parser.parse_args()

    print(f"\033[94mProcessing\033[0m [{args.ingff3}] >>> [{args.outgff}, {args.outlen}]")

    # --- GFF ---
    # GFF3 解析
    gff_data = pd.read_csv(args.ingff3, sep="\t", header=None, comment="#", dtype={0: str, 3: int, 4: int, 8: str})
    gff_data.columns = ['chr', 'source', 'feature', 'start', 'end', 'score', 'strand', 'phase', 'attribute']

    # 重新提取数据：'feature' = type 参数的 dataframe
    gff_data = gff_data[gff_data['feature'] == args.type].loc[:, ['chr', 'attribute', 'start', 'end', 'strand']]

    # 取各 feature 的实际 args.type
    gff_data['attribute'] = gff_data['attribute'].str.extract(r'ID=([^;]+)')

    # 重命名
    gff_data.columns = ['chr', args.type, 'start', 'end', 'strand']

    # 重做索引
    gff_data = gff_data.sort_values(by=['chr', 'start']).reset_index(drop=True)
    gff_data['order'] = gff_data.groupby('chr').cumcount() + 1

    # 输出 GFF
    print(f"\033[94mWriting\033[0m to {args.outgff}...")
    print(gff_data)
    gff_data.to_csv(args.outgff, sep="\t", header=False, index=False)

    # --- LEN ---
    print(f"\033[94mCalculating\033[0m chromosome lengths and gene counts...")
    # 基于 FNA 文件，统计各染色体核苷酸长度
    chr, lens = [], []
    for rec in SeqIO.parse(args.infna, 'fasta'):
        chr.append(rec.id)
        lens.append(len(rec.seq))
    chr2length = pd.DataFrame(data={'chr': chr, 'lens': lens})

    # 对原 gff_data 采用聚合方法，统计 'attribute' 条目数（count）
    chr2count = gff_data.groupby('chr').agg({args.type: 'count'}).reset_index()

    # 合并两个 dataframe
    lens = pd.merge(chr2length, chr2count, on='chr')

    # 输出 LENS
    print(f"\033[94mWriting\033[0m to {args.outlen}...")
    print(lens)
    lens.to_csv(args.outlen, sep="\t", header=False, index=False)

    print("\033[92mPre-Process Step-1 Complete!\033[0m")
    print()

if __name__ == "__main__":
    main()

"""
WGDI Pre-process Step-2

以下预处理代码改进自原 WGDI 项目。
使用前请确认是否安装了 gffread，
是否在 Python 环境中安装了 Pandas 与 BioPython。

GitHub主页: https://github.com/SunPengChuan/wgdi?tab=readme-ov-file
相关论文: https://doi.org/10.1016/j.molp.2022.10.018
脚本原始来源：https://lxz9.com/2021/09/06/wgdi/#%E8%84%9A%E6%9C%AC
"""

from Bio import SeqIO

import pandas as pd
import subprocess as sp
import argparse

def main():
    print()
    print("WGDI Pre-Process \033[94mStep-2\033[0m")

    # 用法说明与参数解析
    parser = argparse.ArgumentParser(description="WGDI Pre-Process Step-2: Get Protein FASTA Files.")
    parser.add_argument('--ingff3', required=True, help='Input GFF3 file name.')
    parser.add_argument('--infna', required=True, help='Input FASTA(Genome) file name.')
    parser.add_argument('--outfna', required=True, help='Output FASTA(CDS) file name.')
    parser.add_argument('--outfaa', required=True, help='Output FASTA(Protein) file name.')
    args = parser.parse_args()

    # GFF解析
    gff_data = pd.read_csv(args.ingff3, sep="\t", header=None, comment="#", dtype=str)
    gff_data.columns = ['chr', 'source', 'feature', 'start', 'end', 'score', 'strand', 'phase', 'attribute']

    # 提取包含mRNA和CDS的数据
    mrna_data = gff_data[gff_data['feature'] == "mRNA"]
    cds_data = gff_data[gff_data['feature'] == "CDS"]

    # 建立基因到转录本(mRNA)的关系
    print(f"\033[94mAnalysing\033[0m make dictionary for gene to transcript...")
    g2t_df = pd.DataFrame(columns=['gene', 'transcript', 'chromosome'])
    g2t_df['gene'] = mrna_data['attribute'].str.extract(r"Parent=([^;]+)")
    g2t_df['transcript'] = mrna_data['attribute'].str.extract(r"ID=([^;]+)")
    g2t_df['chromosome'] = mrna_data['chr']

    # 计算 CDS 区长度，以便选出 CDS 最长的转录本
    print(f"\033[94mAnalysing\033[0m calculating transcript CDS length...")
    cds_data = cds_data.loc[:, ['attribute', 'start', 'end']].copy()
    cds_data['transcript'] = cds_data['attribute'].str.extract(r"Parent=([^;]+)")
    cds_data['start'] = cds_data['start'].astype(int)
    cds_data['end'] = cds_data['end'].astype(int)
    cds_data['length'] = cds_data['end'] - cds_data['start'] + 1
    cds_length_df = cds_data.groupby('transcript')['length'].sum().reset_index()

    # 合并两个dataframe，然后找最长的转录本
    print(f"\033[94mAnalysing\033[0m search the longest CDS...")
    transcript_table = pd.merge(g2t_df, cds_length_df, on='transcript')
    transcript_table = transcript_table.loc[:, ['gene', 'transcript', 'length', 'chromosome']]
    transcript_table_filtered = transcript_table.loc[transcript_table.groupby('gene')['length'].idxmax()]
    transcript_table_filtered['gene_order'] = transcript_table_filtered.groupby('chromosome').cumcount() + 1

    # 收集结果（转录本到基因、染色体编号等的关系）
    print(f"\033[94mAnalysing\033[0m collate result...")
    print(transcript_table_filtered)
    trans2gene = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene']))
    trans2chr = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['chromosome']))
    trans2order = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene_order']))

    # 开始生成CDS/Prot序列
    print(f"\033[94mProcessing\033[0m [{args.ingff3}, {args.infna}] >>> [{args.outfna}, {args.outfaa}]")
    for opt, outfile in [("-x", args.outfna), ("-y", args.outfaa)]:
        sp.run(["gffread", args.ingff3, "-g", args.infna, opt, outfile], check=True)

    # 对生成的序列进行筛选
    # 若序列label与transcript库匹配，则取用
    filtered_records_cds = [
        record for record in SeqIO.parse(args.outfna, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]
    filtered_records_prot = [
        record for record in SeqIO.parse(args.outfaa, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]

    # 将transcript名称改为基因名称
    for idx, record in enumerate(filtered_records_cds):
        record.description = f"index={idx+1},gene={trans2chr[record.id]}g{trans2order[record.id]}"
        record.id = trans2gene[record.id]
    SeqIO.write(filtered_records_cds, args.outfna, "fasta")

    for idx, record in enumerate(filtered_records_prot):
        record.description = f"index={idx+1},gene={trans2chr[record.id]}g{trans2order[record.id]}"
        record.id = trans2gene[record.id]
        record.seq = record.seq.replace(old='.', new="X")
    SeqIO.write(filtered_records_prot, args.outfaa, "fasta")
    
    print("\033[92mPre-Process Step-2 Complete!\033[0m")
    print()

if __name__ == "__main__":
    main()

后记

回来之后，老师主动跟我说了 WGDI 的事情，以及做染色体核心演化的研究方向，于是重新捡起了这篇文章——但我已经忘得一干二净了！

万幸我那时工作留痕，写了这两篇文档，重新捡起了 WGDI，又对文章做了修改与大量补充。

回想起那时的经历——在三个月内恶补 Linux、BioPython 与 JCVI，从 WSL 到 Conda 再到写论文，坐电脑前昏天黑地没个头的日子，忙是真的忙，也非常地疯狂，但也真的有意思，越来越发觉我确实适合这个工作，这就是我应该坚持做的事情。

想到这，因为迷茫、生病而浑浑噩噩度过的本科生涯，也算是没白过了。

2025年8月29日。

---

之前写的代码多少过于凌乱，特别是对 DataFrame 的处理，甚至直接使用 Column 的索引号，而不是名称。两大坨“屎山”说是。

本着“经常洗洗更健康”的原则，我又对这两个脚本进行了优化，并在多个基因组数据中进行试验。目前已经做到能用，且可读性尚可了。

2025年11月10日。

参考文献

• Sun et al. (2022). WGDI: A user-friendly toolkit for evolutionary analyses of whole-genome duplications and ancestral karyotypes. Mol. Plant. doi: https://doi.org/10.1016/j.molp.2022.10.018.

  WGDI的原始文献。

• https://wgdi.readthedocs.io/en/latest/

  WGDI的在线帮助文档。

• https://lxz9.com

  李星泽的个人博客，2023年来已经不再更新，但原有内容仍有重读价值。它也是我附录那两个脚本的借鉴来源。