文章

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

简单介绍WGDI的使用方法，以及预处理过程。

发表于 更新于
作者 HYLi360
27 分钟阅读
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安装它：

1
2

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

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

1
2
3

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

如果经历这种错误：

1
2

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 文件都不是标准格式文件。官方对这两个文件的格式要求如下：

  • 对于 GFF：

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

  • 对于LEN：

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

###
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的第三列来确定。

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

1

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

而WGDI需要的格式是

1

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

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

先解决前五列：

1
2
3
4
5
6
7
8

    # 先读
    data = pd.read_csv(args.input, sep="\t", header=None, skiprows=0, comment='#',
                       dtype={0: str, 2: str, 3: int, 4: int, 6: str, 8: str})
    # 再筛
    data = data[data[2] == args.type]
    data = data.loc[:, [0, 8, 3, 4, 6]]
    # 最后改正
    data[8] = data[8].str.split('=|;|\|',expand=True)[1]

第8行本来想做正则来着，但感觉不如用split直观，而且我用的也不太熟练。用正则显然是最佳方案。

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

1
2

    data = data.sort_values(by=[0, 3]).reset_index(drop=True)
    data['order'] = data.groupby(0).cumcount() + 1

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

1
2
3

	# 生成 LENS
	data[0] = data[0].astype(str).str.strip()
	lens = data.groupby(0).agg({4: 'max', 8: 'count'}).reset_index()

关键在最后一行。最后一行的 data.groupby(0).agg({4: 'max', 8: 'count'}) 我认为是不太常见但很有用的方法，非常适合做数据透视表。其功能是：

  • 通过第0列（染色体）分组；
  • 统计各组在基因终点（染色体末端）的最大值，以及基因数量；
  • 最后重置索引，得到 lens 这个dataframe。

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

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

有人会说，获取蛋白质数据，完全可以直接用gffread来实现。这确实可行，但实际操作的时候，还是会差点事。例如，每一个蛋白质条目变成了转录本，基因与蛋白质的对应关系消失了（特别常见于ensembl版的gff）：

1

> transcript:ENSWHXT00000003709

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

1
2

> g11.t1
> g11.t2

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

###
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
2
3
4
5
6
7
8
9
10
11
12
13
14

###
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最长的那个版本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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

    1
2
3
4
5
6
7
8
9
10
11

        print(f"\033[94mAnalysing\033[0m make dictionary for gene to transcript...")
    attributes_list = list(mrna_data[8])
    gene_to_transcripts = defaultdict(list)
    pattern = r"ID=([^;]+);Parent=([^;]+)"
  
    for i in attributes_list:
        match = re.search(pattern, i)
        if match:
            gene_to_transcripts[match.group(2)].append(match.group(1))
      
    g2t_df = pd.DataFrame([(gene, transcript) for gene, transcripts in gene_to_transcripts.items() for transcript in transcripts], columns= ['gene', 'transcript'])

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

    1
2
3
4
5
6
7
8

        print(f"\033[94mAnalysing\033[0m calculating transcript CDS length...")
    cds_data = cds_data.loc[:, [8, 3, 4]]
  
    cds_data.columns = ['transcript', 'start', 'end']
  
    cds_data.loc[:, 'transcript'] = cds_data['transcript'].str.split('=|;|\|',expand=True)[3]
    cds_data.loc[:, 'length'] = cds_data['end'] - cds_data['start'] + 1
    cds_length_df = cds_data.groupby('transcript')['length'].sum().reset_index()

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

    1
2
3
4
5
6
7
8
9
10
11

        # 将结果1套用到结果2
    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']]
      
    # 选择CDS总长最大的组合
    transcript_table_filtered = transcript_table.loc[transcript_table.groupby('gene')['length'].idxmax()]
  
    # 生成干净的字典
    print(f"\033[94mAnalysing\033[0m collate result...")
    trans2gene = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene']))

  • 最后，在gffread的结果里：

    1. 删掉字典里没有出现的转录本；
    2. 将转录本名字换成基因名字。
1
2
3
4
5
6
7
8
9

    filtered_records = [
        record for record in SeqIO.parse(args.outfaa, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]

    for record in filtered_records:
        record.id = trans2gene[record.id]
        record.description = ""
    
    SeqIO.write(filtered_records, args.outfaa, "fasta")

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

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

正式的数据处理过程

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

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

因为我们只有蛋白质序列，这里我们跑一下 -d 和 -icl 就行了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

# 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。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 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风格）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# 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

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

A：可以直接动LEN文件，修改染色体及其顺序。改名 可以在第一步里完成（通过--chrrm参数）。

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

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

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

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

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

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

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

A：可以。你只需要运行wgdi -xxxx ? > total.conf（-xxxx = -d/-icl/…），即可生成完整的配置文件模板，然后就可以根据自己的需求，补充缺失的参数。当然，如果你想为每一个功能都准备一个文件，也是完全可以的。

附录

处理脚本源码

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

"""
WGDI Pre-Process Step-1

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

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

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

    # receive args/接受参数
    parser = argparse.ArgumentParser(description="WGDI Pre-Process Step-1: Get GFF and LEN Files.")

    parser.add_argument('--input', 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='Entrie type.')

    parser.add_argument('--chrrm', required=False, help='Optional. Remove the content on First (Chr) column.')
    parser.add_argument('--chradd', required=False, help='Optional. Add the content on head of First (Chr) column.')

    args = parser.parse_args()

    # main process logic/主处理逻辑
    print(f"\033[94mProcessing\033[0m [{args.input}] >>> [{args.outgff}, {args.outlen}]")

    # generate/生成 GFF
    data = pd.read_csv(args.input, sep="\t", header=None, skiprows=0, comment='#',
                       dtype={0: str, 2: str, 3: int, 4: int, 6: str, 8: str})
    data = data[data[2] == args.type]
    data = data.loc[:, [0, 8, 3, 4, 6]]
    data[8] = data[8].str.split(r'=|;|\|',expand=True)[1]
    data = data.sort_values(by=[0, 3]).reset_index(drop=True)
    data['order'] = data.groupby(0).cumcount() + 1

    if args.chrrm is not None:
        data[0] = data[0].str.replace(args.chrrm,'',regex=True)
    if args.chradd is not None:
        data[0] = args.chradd + data[0]

    print(f"\033[94mWriting\033[0m to {args.outgff}...")
    data.to_csv(args.outgff, sep="\t", header=None, index=False)

    # generate/生成 LEN
    print(f"\033[94mCalculating\033[0m chromosome lengths and gene counts...")
    data[0] = data[0].astype(str).str.strip()
    lens = data.groupby(0).agg({4: 'max', 8: 'count'}).reset_index()
    print(f"\033[94mWriting\033[0m to {args.outlen}...")
    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()


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

"""
WGDI Pre-process Step-2

以下预处理代码改进自原 WGDI 项目。
使用前请确认是否在 Python 环境中安装了 argparse、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 subprocess
import re
from collections import defaultdict
from Bio import SeqIO
import pandas as pd
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('--ingff', required=True, help='Input GFF3 file name.')
    parser.add_argument('--infna', required=True, help='Input FASTA(Genome) file name.')
    parser.add_argument('--out', required=True, help='Output FASTA(Protein/CDS) file name.')
    parser.add_argument('--mode', required=True, help="Extract Protein(p) or CDS(c).", choices=['p', 'c'])
    args = parser.parse_args()

    # GFF与CDS统计
    data = pd.read_csv(args.ingff, sep="\t", header=None, comment="#", dtype=str)
    mrna_data = data[data[2] == "mRNA"]
    cds_data = data[data[2] == "CDS"]

    # 建立基因到转录本的关系
    print(f"\033[94mAnalysing\033[0m make dictionary for gene to transcript...")
    attributes_list = list(mrna_data[8])
    gene_to_transcripts = defaultdict(list)
    pattern = r"ID=([^;]+);Parent=([^;]+)"
    for i in attributes_list:
        match = re.search(pattern, i)
        if match:
            gene_to_transcripts[match.group(2)].append(match.group(1))
    g2t_df = pd.DataFrame(
        [(gene, transcript) for gene, transcripts in gene_to_transcripts.items() for transcript in transcripts],
        columns=['gene', 'transcript']
    )

    # 计算 CDS 区长度，以便选出 CDS 最长的转录本
    print(f"\033[94mAnalysing\033[0m calculating transcript CDS length...")
    cds_data = cds_data.loc[:, [8, 3, 4]].copy()
    cds_data.columns = ['transcript', 'start', 'end']
    cds_data['transcript'] = cds_data['transcript'].str.split(r'=|;|\|', expand=True)[3]
    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()

    # 找最长的转录本
    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']]
    transcript_table_filtered = transcript_table.loc[transcript_table.groupby('gene')['length'].idxmax()]

    # 收集结果
    print(f"\033[94mAnalysing\033[0m collate result...")
    trans2gene = dict(zip(transcript_table_filtered['transcript'], transcript_table_filtered['gene']))

    # 开始提取序列
    print(f"\033[94mProcessing\033[0m [{args.ingff}, {args.infna}] >>> [{args.out}]")
    if args.mode == 'c':
        command = f"gffread {args.ingff} -g {args.infna} -x {args.out}"
    else:
        command = f"gffread {args.ingff} -g {args.infna} -y {args.out}"
    subprocess.run(command, shell=True)

    # 对生成的序列进行筛选
    filtered_records = [
        record for record in SeqIO.parse(args.out, "fasta") if record.id in list(transcript_table_filtered['transcript'])
    ]

    for record in filtered_records:
        record.id = trans2gene[record.id]
        record.description = ""
    
    SeqIO.write(filtered_records, args.out, "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日。

参考文献

生物信息学, 工具
wgdi bioinfo-tools genome
本文由作者按照 CC BY 4.0 进行授权