实验篇——鉴定不同基因的复制模式

实验篇——鉴定不同基因的复制模式

---

前言

这几天了解了一下，如何使用“DupGen_finder"工具进行基因复制模式的鉴定。其作者的相关代码发布在
https://github.com/qiao-xin/DupGen_finder

一、工具下载

在Linux中的操作

#将代码克隆过来
git clone https://github.com/qiao-xin/DupGen_finder.git


#环境配置
make
chmod 775 DupGen_finder.pl
chmod 775 DupGen_finder-unique.pl
chmod 775 set_PATH.sh
./set_PATH.sh


#查看
./DupGen_finder.pl
Usage: DupGen_finder.pl -i data_directory -t target_species -c outgroup_species -o output_directory
#####################
Optional:
-a 1 or 0(are segmental duplicates ancestral loci or not? default: 1, yes)
-d number_of_genes(maximum distance to call proximal, default: 10)
#####################
The following are optional MCScanX parameters:
-k match_score(cutoff score of collinear blocks for MCScanX, default: 50)
-g gap_penalty(gap penalty for MCScanX, default: -1)
-s match_size(number of genes required to call a collinear block for MCScanX, default: 5)
-e e_value(alignment significance for MCScanX, default: 1e-05)
-m max_gaps(maximum gaps allowed for MCScanX, default: 25)
-w overlap_window(maximum distance in terms of gene number, to collapse BLAST matches for MCScanX, default: 5)

二、数据文件准备

假设对乌药Lindera aggregata(LIN)中的基因复制模式进行划分，使用山胡椒Lindera glauca (GLA)作为外类群。

总共要准备四个文件：

1. gff文件

1. 目标基因组(LIN)的gff文件(LIN.gff)

2. 目标基因组(LIN)与外类群基因组(GLA)的合并gff文件(LIN_GLA.gff)

cat LIN.gff GLA.gff >LIN_GLA.gff

注意 ：
这里要求的gff文件中的数据格式如下：

LIN-chr1	LaggChr1G00000010	5815	15157
LIN-chr1	LaggChr1G00000010.1	5815	15157
LIN-chr1	LaggChr1G00000010.1	5815	5928
LIN-chr1	LaggChr1G00000010.1	9867	9955
LIN-chr1	LaggChr1G00000010.1	14803	15157
LIN-chr1	LaggChr1G00000020	55363	57323
LIN-chr1	LaggChr1G00000020.1	55363	57323
LIN-chr1	LaggChr1G00000020.1	55363	55513
LIN-chr1	LaggChr1G00000020.1	57214	57323
LIN-chr1	LaggChr1G00000030	79251	85730

但实际上一般的gff文件的数据格式为:

我们只需要"所属染色体片段的名称"、“基因ID”、“起点位置”、"终点位置"这四个信息，故要对原本的gff文件进行修改
代码示例如下：

import pandas as pd

def gff_xiugai(wenjian,new_wenjian,qianzhui):
    wenjian = pd.read_csv(wenjian,skiprows=1,sep="\t",header=None)
    wenjian["ID"] = wenjian[8].str.extract(r"ID=([^;\n]+)")
    wenjian["Parent"] = wenjian[8].str.extract(r"Parent=([^;\n]+)")
    #如果没有ID，用Parent的值代替(具体情况具体分析）

    wenjian["ID"] = wenjian["ID"].fillna(wenjian["Parent"])
    wenjian[0] = qianzhui + wenjian[0]
    new_df = pd.DataFrame(wenjian[[0,"ID",3,4]])
    new_df.to_csv(new_wenjian,sep="\t",header=False,index=False)

#示例
gff_xiugai(r"D:\yuceji\04.Lindera_glauca\Lindera_glauca.gff","D:\GLA2.gff","GLA-")

函数的内容可以随着具体的数据格式而改变，以适应具体的要求。

2. blast文件

3. 自身比对得出的blastp文件(LIN.blast)；
   将目标物种基因组(LIN)的pep文件当作参考数据库与自身比对

4. 目标物种(LIN)和外类群物种(GLA)的blastp文件(LIN_GLA.blast);
   将GLA的pep文件当作参考数据库与LIN的pep文件比对。

对于blast文件的获得，可以在linux服务器上使用blast软件获得

1. 自身比对文件：

#创建数据库
makeblastdb -in LIN.pep -dbtype prot -title LIN -parse_seqids -out LIN
#比对
blastp -query LIN.pep -db LIN -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out LIN.blast

2. 跨基因组比对文件：

#创建数据库
makeblastdb -in GLA.pep -dbtype prot -title GLA -parse_seqids -out GLA
#比对
blastp -query LIN.pep -db GLA -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out LIN_GLA.blast

3. blast文件格式

每一列代表：

1. Query acc.ver：查询序列的Accession号或标识符。

2. Subject acc.ver：参考序列的Accession号或标识符。

3. % Identity：查询序列与目标序列之间的匹配相似度百分比。(表示在对齐的位置上，两个序列中相同的核苷酸或氨基酸的比例。)

4. Alignment length：对齐的长度，即匹配上的核苷酸或氨基酸数量。

5. Mismatches：在对齐过程中，查询序列和目标序列之间的不匹配数量。这包括替换和插入/删除错误。

6. Gap opens：在对齐中开放的间隙数，表示在查询序列或目标序列中存在的缺失或插入发生的次数。

7. q. start：查询序列的开始位置，在对齐中的序列坐标上。

8. q. end：查询序列的结束位置，在对齐中的序列坐标上。

9. s. start：参考序列的开始位置，在对齐中的序列坐标上。

10. s. end：参考序列的结束位置，在对齐中的序列坐标上。

11. Evalue（期望值）：表示在随机情况下获得与查询序列具有相同或更好匹配的目标序列的期望数量。较小的E值表示匹配的统计意义较高。

12. Bit score：比特分数是根据匹配的统计意义计算得出的一个评分。它表示匹配的强度和可靠性，较高的比特分数表示更显著的匹配。

4. 额外补充：

在Linux中可以使用screen工具来创建会话,使程序在后台运行。
在远程连接或断开连接的情况下，也可以保持会话的运行。

#创建一个新会话(name为自定义的会话名)
screen -S name   

#查看创建的会话（会显示Attached或Detached状态)
screen -ls

#恢复会话窗口(继续查看后台运行的进程）
screen -r name

更多详细内容，参考
https://blog.csdn.net/weixin_41010198/article/details/121166007

三、执行

数据文件准备好后，将他们放于同一目录下。

1. DupGen_finder.pl

使用命令如下(在DupGen_finder工具的目录下）：

perl DupGen_finder.pl -i data -t LIN -c GLA -o output

命令解释：

1. perl：表示使用 Perl 解释器来执行 DupGen_finder.pl 脚本。

2. DupGen_finder.pl：这是一个 Perl 脚本，用于运行 DupGen_finder 工具进行基因复制模式的鉴定。

3. -i data：指定输入目录的路径，data表示准备好的数据文件所在目录的路径。

4. -t LIN：指定用于鉴定基因复制模式的待比较基因组的简称或标识符，如LIN

5. -c GLA：指定用于鉴定基因复制模式的参考基因组的简称或标识符,如GLA

6. -o output：指定输出文件的路径

2. DupGen_finder-unique.pl

为了减少不同模式中的冗余复制基因，可使用DupGen_finder-unique这个更严格的版本
命令如下：

perl DupGen_finder-unique.pl -i data -t LIN -c GLA -o output

四、结果文件

使用DupGen_finder-unique这个更严格的版本

# ls
LIN.collinearity            LIN.pairs.stats-unique     LIN.tandem.pairs-unique
LIN.dispersed.genes-unique  LIN.proximal.genes-unique  LIN.transposed.genes-unique
LIN.dispersed.pairs-unique  LIN.proximal.pairs-unique  LIN.transposed.pairs-unique
LIN.gff.sorted              LIN.singletons             LIN.wgd.genes-unique
LIN_GLA.collinearity        LIN.tandem.genes-unique    LIN.wgd.pairs-unique

1. Duplicated gene pairs

WGD复制（LIN.wgd.pairs-unique）
串联重复（LIN.tandem.pairs-unique）
邻近复制（LIN.proximal.pairs-unique）
转座复制（LIN.transposed.pairs-unique）
分散复制（LIN.dispersed.pairs-unique ）

示例文件内容如下：
Duplicate 1:第一个重复基因的标识符
Location: 第一个重复基因的位置
Duplicate 2: 第二个重复基因的标识符
Location: 第二个重复基因的位置
E-value: 相关性或相似性评估的E值

2. LIN.singletons

该文件包含在目标物中没有同源基因的基因

文件内容如下：

3. LIN.pairs.stats-unique

包含不同模式中得到的复制基因对的数量，相当于一个统计文件

文件内容如下：

Types: 基因对的类型
NO. of gene pairs: 对应类型的基因对数量

4. Collinearity files

共线性文件
LIN.collinearity (自身）
LIN_GLA.collinearity （与外类群的）

文件内容如下：

---

总结

这篇文章算是我完成 "鉴定不同基因的复制模式"这一任务的总结，从工具的下载到数据文件的准备再到最后的执行，成功得到了结果文件。（在这一过程中，可能会因为操作环境的不同而导致出现错误）。

仙人抚我顶，结发受长生

–2023-8-8 实验篇