[Python知识库]WGDI 分析全基因组复制事件完整流程

简介

WGDI（全基因组重复识别），一种基于 Python 的命令行工具，可让研究人员深入了解递归多倍化并进行跨物种基因组比对分析。

官方文档
下一
篇会选择一个物种的分析结果做示例。

安装

## 1.使用conda安装
conda install -c bioconda  wgdi

## 2.使用pip安装
pip install wgdi

## 3.本地安装
git clone https://github.com/SunPengChuan/wgdi.git
cd wgdi
python setup.py install

依赖软件

PAML | MAFFT | MUSCLE | PAL2NAL | IQTREE

使用pip下载完成 需要配置文件目录
wgdi -conf /？ >conf.ini
里面是默认的文件路径，如果不对应打开修改即可
再次输入 wgdi -conf conf.ini
就将配置环境导入到程序中了

分析

数据预处理

数据处理是很有必要的，如果格式不正确，后面的分析很有可能会报错，大家可以自行处理数据得到gff文件以及基因组len文件

下面提供wgdi作者以及写的处理脚本,具体脚本内容将放到本文末尾

01.getgff.py
02.gff_lens.py
03.seq_newname.py
generate_conf.py # 脚本由徐洲更制作。下载链接

序列比对

通过blastp或者diamond进行蛋白之间比对。

makeblastdb -in ath.pep.fa -dbtype prot
blastp -num_threads 50 -db ath.pep.fa -query ath.pep.fa -outfmt 6 -evalue 1e-5 -num_alignments 20  -out ath.blastp.txt & 

点图

使用命令进入文件夹取出参数文件，主要是通过前面的[] 去识别程序所以可以将下列所有步骤放到一个文本文件中。

wgdi -d help >> total.conf

修改配置文件

[dotplot]
blast = blast file # blast
blast_reverse = false
gff1 =  gff1 file # gff和len使用python程序可以提取，格式有要求
gff2 =  gff2 file
lens1 = lens1 file
lens2 = lens2 file
genome1_name =  Genome1 name
genome2_name =  Genome2 name
multiple  = 1 # 最好同源的个数(用红色表示)，不清楚就写1
score = 100
evalue = 1e-5
repeat_number = 10 # 基因相对另一个基因组同源基因的数量限制
position = order # 与end相比，end的点图相对比较集中。
blast_reverse = false # 比对中基因对顺序是否翻转,
ancestor_left = none # 点图左侧物种的祖先染色体区域
ancestor_top = none # 点图顶部物种的祖先染色体区域
markersize = 0.5 # 调整图中点的大小
figsize = 10,10 # 保存图片大小比例
savefig = savefile(.svg,.png,.pdf) # 类型： 默认： *.PNG保存图片支持png、pdf、svg格式

运行命令 wgdi -d total.conf

任意修改lens1 和lens2的染色体的数量和顺序，可以得到任意染色体间的点图。
如果想修改图片中染色体的位置，则可以改动len文件进行操作。
也可以单独选择几条染色体做点图观察，都是通过改动len文件进行的

点图解读：
点图规则：点图是根据左边基因组的基因，最好同源的点为红色，次好的四个基因为蓝色，剩余部分（同源基因有限制个数）为灰色，是根据打分结果进行排序的。
（1）点图需要横向看和纵向看：这个点图是物种自身比对，所以只需横向看。片段深度应该为2，但最好同源个数为1，即红色点没有集中分布在某个片段上。常常可以认两个片段很相似，再加上自身。所以，认定为最近发生的加倍为三倍化。除此之外，蓝色或灰色的片段很少有，表明更古老的加倍很不明显。
（2）对角线上，本不应该出现片段。自身比自身显然是最好匹配，这些点已经去掉了。而剩余的这些由于串联重复造成的。即该基因临近的位置同源性很高，打在了对角线附近。可以计算这些串联重复的ks值，估算重复片段的爆发时间。
（3）这个点图是后续跑共线性的基石。所以，score， evalue， repeat_number判定的同源点的数量是共线性打分矩阵的来源。

可以看到点图的结果中，共线性片段，按照blast进行排序，红色，蓝色 ，灰色是最后剩下的，左边相对上面的义工只有十个同源基因对应
看到点图中共线性片段大多是由红色和蓝色构成的。

共线性

使用命令进入文件夹取出参数文件。wgdi -icl ? >> total.conf

[collinearity]
gff1 = gff1 file
gff2 = gff2 file
lens1 = lens1 file
lens2 = lens2 file
blast = blast file
blast_reverse = false
multiple  = 1 # 用红点显示的同源基因的最佳数量
process = 8 # 根据点图中的颜色分配不同的分数，红色默认为 50，蓝色为 40，灰色为 25。
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 # 输出文件

运行 wgdi -icl total.conf

Ks值计算

wgdi -ks ? >> total.conf

[ks]
cds_file =      cds file # 两个物种或物种自身的cds序列，如果要计算多个物种，可以将它们文件cat到一起
pep_file =      pep file # 两个物种或物种自身的pep序列
align_software = muscle # 选择mafft or muscle
pairs_file = gene pairs file #代表其中一种: colinearscan输出结果 / mcscanx输出结果 / collinearity输出结果 / tab分割的基因对
ks_file = ks result # 输出文件，格式有要求

运行 wgdi -ks total.conf

输出的ks结果，基因对算不出来的ks值，默认为-1。

Blockinfo

因为共线性和ks值结果不容易展示，而且，常常需要对共线性结果筛选，因此，将这些信息汇总到一个表中，方便删选。
这一步主要是将前面跑的共线性和ks值统一到一个文件当中。

首先将参数保存到文件当中
wgdi -bi ? >> total.conf

[blockinfo]
blast = blast file
gff1 =  gff1 file
gff2 =  gff2 file
lens1 = lens1 file
lens2 = lens2 file
collinearity = collinearity file
score = 100 # 和之前设置的保持一致
evalue = 1e-5 # 和点图、共线性计算的配置文件保持一致
repeat_number = 20 # 和点图、共线性计算的配置文件保持一致
position = order # 仅有order选项
ks = ks file
ks_col = ks_NG86 # 根据之前ks得出的列表选择，软件会生成ks_NG86或ks_YN00，也可以用其他方式计算的结果单独一列
savefile = block_information.csv # 将共线性和ks值汇总到一起

运行 wgdi -bi total.conf

结果生成一个表block_information.csv

第1列: id 共线性的结果的唯一标识
第2.4.5列: chr1,start1,end1 参考基因组（点图的左边）的共线性范围
第3.6.7列: chr2,start2,end2 参考基因组（点图的上边）的共线性范围
第8列: pvalue 评估共线性的结果，常常认为小于0.01的更合理些
第9列: length 共线性片段长度
ks_median 共线性片段上所有基因对ks的中位数（主要用来评判ks分布的）
ks_average 共线性片段上所有基因对ks的平均值
homo1 homo2 homo3 homo4 homo5 与multiple参数有关，即homo+multiple
主要规则是：基因对如果在点图中为红色，赋值为1，蓝色赋值为0，灰色赋值为-1。共线性片段上所有基因对赋值后求平均值，这样就赋予该共线性一个-1,1的值。如果共线性的点大部分为红色，那么该值接近于1。可以作为共线性片段的筛选。
block1,block2 分别为共线性片段上基因order的位置。
ks 共线性片段的上基因对的ks值
density1, density2 共线性片段的基因分布密集程度。值越小表示稀疏;如何计算的呢？就是length 除以对应的srart1 end1 差值

Correspondence

wgdi -c ? >> total.conf

[correspondence]
blockinfo =  block_information.csv #上一步blockinfo的输出csv文件作为输入文件
lens1 = lens1 file
lens2 = lens2 file
tandem = (true/false) #前一步展示过了不需要分析改为false
tandem_length = 200 #看对角线有没有加倍造成的片段，如果对角线上的ks值趋于0，那就不属于。
pvalue = 0.05 #要提取好的共线性 作者设置为了0.2
block_length = 5
multiple  = 1
homo = 0,1
savefile = block_information.new.csv # 名称自定义(*.csv即可)

wgd -c total.conf

筛选完的结果，是将对角线上的片段去除，剩余一些需要的片段。
得出筛选完之后的结果会发现新的csv文件会比之前的文件要小，主要就是将对角线上的片段给去除了。

blockks

wgdi -bk ? > blockks.conf

[blockks]
lens1 = lens1 file
lens2 = lens2 file
genome1_name =  Genome1 name
genome2_name =  Genome2 name
blockinfo = block_information.new.csv
pvalue = 0.05
tandem = true
tandem_length = 200
markersize = 1 #点的大小
area = 0,2 #显示范围
block_length =  minimum length
figsize = 8,8 #点图大小的比例值
savefig = save image #保存图片，默认： *.PNG保存图片支持png、pdf、svg格式

wgdi -bk total.conf

最终得到的结果会发现对角线上的点被去除，但是还是会有一些点会没有去除。但是也可以继续进行后续操作。

KsPeaks

计算Ks峰值

wgdi -kp ? >> total.conf

[kspeaks]

 blockinfo = block_information.new.csv
 pvalue = 0.05 #设置可自行调整
 tandem = true
 block_ length = int number
 ks_area = 0,10
 multiple = 1
 homo = 0,1 #结合最开始的点图，红色点代表1，查看共线性片段如果红色居多，整个共线性片段的值就接近于1，可以写成0.3/0.5,1
 fontsize = 9
 area = 0,3 
 figsize = 10,6.18
 savefig = saving image
 savefile = ks_median.distri.csv

wgdi -kp total.conf

结果图片如果发现前面趋于0 的值比较多，查看前面得出的点图，原因是由于Correspondence这一步设置的tandom值长度偏小，
可以继续拟合或者直接在ks_median.distri.csv文件中进行筛选，将ks_median和ks_average进行从小到大排序，将数值过低的去除，它们还有一个特性就是
chr1以及chr2 都是自己对自己。
将修改结果重新输入[kspeaks] blockinfo = 修改后的(*.csv)

PeaksFit

峰拟合，需要结合KsPeaks和blockks的结果看峰为一个还是多个，mode = median建议用每个block的中位数的ks值去做拟合。

wgdi -pf ? >> total.conf

[peaksfit] 
blockinfo = block information  #经过kspeaks处理的单个峰的输出文件(*.csv)
mode = median  #分为，对应Ks峰图的三块median，average，total
bins_number = 200  #柱状图的bins
ks_area = 0,10  #ks值的范围，可以去掉一些奇异的ks值
fontsize = 9 
area = 0,3  #横轴的范围
figsize = 10,6.18 
savefig = saving image

wgdi -pf total.conf

跑完得到 拟合优度以及对应的曲线
R-square: 0.9618516780863838
The gaussian fitting curve parameters are :
1.305477901404146 | 2.1185082504069626 | 0.41222580475422466
其中：
中间的第二个值就是峰值的大小，数据保存到单独文件中，就可以绘制多个物种种内和种间的ks峰值了。

关于一个物种发生多次加倍事件如何获取Ks峰值的策略:
每次加倍保留的共线性片段是一定的，是某一次产生的。可以通过Ks值进行区分。
应该对其分开处理，对共线性片段筛选（去除对角线的Ks峰或者有多次加倍的颜色发生变化的需要去除之后在进行拟合）

脚本化运行

# 启动环境
conda activate wgdi 
# 进入目录并解压
cd Arabidopsis_thaliana
gunzip *.gz

# 将generate_conf.py配置到环境变量
generate_conf.py -p ath Athaliana_447_TAIR10.fa Athaliana_447_Araport11.gene.gff3

head -n 3 ath.gff
Chr1	AT1G01010.Araport11.447	3631	5899	+	1	AT1G01010.1.Araport11.447
Chr1	AT1G01020.Araport11.447	6788	9130	-	2	AT1G01020.1.Araport11.447
Chr1	AT1G01030.Araport11.447	11649	13714	-	3	AT1G01030.1.Araport11.447

# 用 sed 删除多余信息
sed -i -e 's/.Araport11.447//' -e 's/.Araport11.447//' ath.gff

# 另外，pep.fa, cds.fa里面包含所有的基因，且命名特别的长， 同时我们也不需要基因中的 ".1", ".2" 这部分信息

head -n 5 Athaliana_447_Araport11.cds_primaryTranscriptOnly.fa
>ATCG00500.1 pacid=37375748 polypeptide=ATCG00500.1 locus=ATCG00500 ID=ATCG00500.1.Araport11.447 annot-version=Araport11
ATGGAAAAATCGTGGTTCAATTTTATGTTTTCTAAGGGAGAATTGGAATACAGAGGTGAGCTAAGTAAAGCAATGGATAG
TTTTGCTCCTGGTGAAAAGACTACTATAAGTCAAGACCGTTTTATATATGATATGGATAAAAACTTTTATGGTTGGGATG
AGCGTTCTAGTTATTCTTCTAGTTATTCCAATAATGTTGATCTTTTAGTTAGCTCCAAGGACATTCGCAATTTCATATCG
GATGACACCTTTTTTGTTAGGGATAGTAATAAGAATAGTTATTCTATATTTTTTGATAAAAAAAAAAAAATTTTTGAGAT

# 借助于seqkit, 对原来的数据进行了筛选和重命名
seqkit grep -f <(cut -f 7 ath.gff ) Athaliana_447_Araport11.cds_primaryTranscriptOnly.fa  | seqkit seq --id-regexp "^(.*?)\\.\\d" -i > ath.cds.fa

seqkit grep -f <(cut -f 7 ath.gff ) Athaliana_447_Araport11.protein_primaryTranscriptOnly.fa  | seqkit seq --id-regexp "^(.*?)\\.\\d" -i > ath.pep.fa

# 通过NCBI的BLASTP 或者DIAMOND 进行蛋白之间相互比对，输出格式为 -outfmt 6
makeblastdb -in ath.pep.fa -dbtype prot
blastp -num_threads 50 -db ath.pep.fa -query ath.pep.fa -outfmt 6 -evalue 1e-5 -num_alignments 20  -out ath.blastp.txt & 

# 绘制点阵图

脚本

getgff.py

import sys
import pandas as pd

data = pd.read_csv(sys.argv[1], sep="\t", header=None,skiprows=0)
data = data[data[2] == 'mRNA']
data = data.loc[:, [0, 8, 3, 4, 6]]
data[8] = data[8].str.split(':|=|;|\|',expand=True)[1]
# data.drop_duplicates(subset=[8], keep='first', inplace=True)
data[0] = data[0].str.replace('At','',regex=True)
data.to_csv(sys.argv[2], sep="\t", header=None, index=False)

gff_lens.py

#!usr/bin/env python
# conding: utf-8
import sys

import pandas as pd

data = pd.read_csv(sys.argv[1], sep="\t", header=None)

new = data[1].str.split('.').str
data['id'] = new[0].values
data['cha'] = data[3]-data[2]

for name, group in data.groupby(['id']):
    if len(group) == 1:
        continue
    ind = group.sort_values(by='cha', ascending=False).index[1:].values
    data.drop(index=ind, inplace=True)

# data[2] =data[2].astype(int)
# data[3] =data[3].astype(int)
# for name, group in data.groupby(0):
#     if len(group) == 1:
#         continue
#     start=0
#     # print(group.head())
#     group = group.sort_values(by=[2,3],ascending=[True,False])
#     for index, row in group.iterrows():
#         # print(row)
#         if row[3]>start or row[2]>start:
#             start=max([row[3],row[2]])
#             continue
#         data.drop(index=index, inplace=True)
    # ind = group.sort_values(by='cha', ascending=False).index[1:].values
    #print(name)
    # print(group.sort_values(by='cha',ascending=False))

    

# data = data[data[1].str.contains('\.mRNA1$')]
data['order'] = ''
data['newname'] = ''
print(data.head())
for name, group in data.groupby([0]):
    number = len(group)
    group = group.sort_values(by=[2])
    data.loc[group.index, 'order'] = list(range(1, len(group)+1))
    data.loc[group.index, 'newname'] = list(
        ['ath2s'+str(name)+'g'+str(i).zfill(5) for i in range(1, len(group)+1)])
data['order'] = data['order'].astype('int')
data = data[[0, 'newname', 2, 3, 4, 'order', 1]]
print(data.head())
data = data.sort_values(by=[0, 'order'])
data.to_csv(sys.argv[2], sep="\t", index=False, header=None)
lens = data.groupby(0).max()[[3, 'order']]
lens.to_csv(sys.argv[3], sep="\t", header=None)

seq_newname.py

import sys
import pandas as pd
from Bio import SeqIO

data = pd.read_csv(sys.argv[1], sep="\t", header=None, index_col=6)
data.index = data.index.str[:-3]
id_dict = data[1].to_dict()
print(data.head())
seqs = []
n = 0
for seq_record in SeqIO.parse(sys.argv[2], "fasta"):
    if seq_record.id in id_dict:
        seq_record.id = id_dict[seq_record.id]
        n += 1
    else:
        continue
    seqs.append(seq_record)
SeqIO.write(seqs, sys.argv[3], "fasta")
print(n)

generate_conf.py

#!/usr/bin/env python3

# GFF must have CDS feature
# GFF must have ID and Parent in column 9

import re
import argparse
from collections import defaultdict
from collections import OrderedDict

def get_parser():
    """Get options"""
    parser = argparse.ArgumentParser()
    parser.add_argument('fasta',
                        help="fasta file name")
    parser.add_argument('gff',
                        help="gff file name")
    parser.add_argument('-p','--prefix', type=str, default="output",
                        help="prefix for ouput ")

    return parser



# get the fasta  len
def get_fasta_len(fasta):
    fasta_dict = OrderedDict()
    handle = open(fasta, "r")
    active_sequence_name = ""
    for line in handle:
        line = line.strip()
        if not line:
            continue
        if line.startswith(">"): 
            active_sequence_name = line[1:]
            active_sequence_name = active_sequence_name.split(" ")[0]
        if active_sequence_name not in fasta_dict:
            fasta_dict[active_sequence_name] = 0
            continue
        sequence = line
        fasta_dict[active_sequence_name] += len(sequence)
    handle.close()
    return fasta_dict

# parse the gff 
def parse_gff(gff):

    gene_dict = OrderedDict()
    tx_pos_dict = defaultdict(list)
    CDS_dict = defaultdict(list)

    handle = open(gff, "r")

    for line in handle:
        if line.startswith("#"):
            continue
        content = line.split("\t")
        if len(content) <= 8:
            continue
        #print(content)
        if content[2] == "transcript" or content[2] == "mRNA":

            # use regual expression to extract the gene ID
            # match the pattern ID=xxxxx; or ID=xxxxx

            tx_id = re.search(r'ID=(.*?)[;\n]',content[8]).group(1)
            tx_parent = re.search(r'Parent=(.*?)[;\n]',content[8]).group(1)
            tx_parent = tx_parent.strip() # remove the 'r' and '\n'

            # if the parent of transcript is not in the gene_dict, create it rather than append
            if tx_parent in gene_dict:
                gene_dict[tx_parent].append(tx_id)
            else:
                gene_dict[tx_parent] = [tx_id]
            tx_pos_dict[tx_id] = [content[0],content[3], content[4], content[6] ]
        # GFF must have CDS feature
        if content[2] == 'CDS':
            width = int(content[4]) - int(content[3])
            CDS_parent = re.search(r'Parent=(.*?)[;\n]',content[8]).group(1)
            CDS_parent = CDS_parent.strip() # strip the '\r' and '\n'
            CDS_dict[CDS_parent].append(width)
    handle.close()
    return [gene_dict, tx_pos_dict, CDS_dict]

if __name__ == "__main__":
    parser = get_parser()
    args = parser.parse_args()
    fa_dict = get_fasta_len( args.fasta)
    gene_dict, tx_pos_dict, CDS_dict= parse_gff( args.gff )
    gene_count = {}

    # outfile
    len_file = open(args.prefix + ".len", "w")
    gff_file = open(args.prefix + ".gff", "w")

    for gene, txs in gene_dict.items():
        tmp = 0
        for tx in txs:
            tx_len = sum(CDS_dict[tx])
            if tx_len > tmp:
                lst_tx = tx
                tmp = tx_len
        tx_chrom = tx_pos_dict[lst_tx][0]
        if tx_chrom not in gene_count:
            gene_count[tx_chrom] = 1
        else:
            gene_count[tx_chrom] += 1
        tx_start = tx_pos_dict[lst_tx][1]
        tx_end   = tx_pos_dict[lst_tx][2]
        tx_strand = tx_pos_dict[lst_tx][3]

        print("{chrom}\t{gene}\t{start}\t{end}\t{strand}\t{order}\t{tx}".format(\
                chrom=tx_chrom,gene=gene,start=tx_start,end=tx_end,strand=tx_strand,order=gene_count[tx_chrom],tx=lst_tx), file=gff_file )

    for chrom,lens in fa_dict.items():
        print("{chrom}\t{lens}\t{count}".format(\
                chrom=chrom,lens=lens,count=gene_count[chrom]), file=len_file)
    len_file.close()
    gff_file.close()

参考文献

Sun P, Jiao B, Yang Y, et al. WGDI: A user-friendly toolkit for evolutionary analyses of whole-genome duplications and ancestral karyotypes[J]. bioRxiv, 2021.