实验篇——Ka/Ks分析

实验篇——Ka/Ks分析

---

前言

鉴定不同基因的复制模式
本文得到的共线性基因对文件 来自于上一篇文章中的LIN.collinearity共线性文件

参考文章：
https://www.jianshu.com/p/a3a39b2f341b

---

一、名词解释

Ka/Ks分析是一种用于评估基因进化速度的方法。较高的Ka/Ks比值可能表明基因在进化过程中经历了功能上的重要变化，而较低的Ka/Ks比值可能表明基因在进化过程中具有保守的功能。

更多相关的解释信息

二、实操

1. 安装软件

安装ParaAT、KaKs_Calculator2.0以及muscle（或者其它序列比对软件） 软件，并配置环境
或者在TBtools软件中操作

2. 准备文件

准备存储共线性基因对的文件 homolog、cds文件、pep文件（注意，cds及pep为fasta格式，“>”后面只接基因ID号）

以LIN.collinearity共线性文件为例

得到homolog文件

cat LIN.collinearity | grep "Lagg" | awk '{print $3"\t"$4}' >LIN.homolog

得到cds文件
原来的文件格式：

sed 's/^>\([^\t]*\)\t.*/>\1/' Lagg.gene.cds >Lagg.cds

要求的文件格式：

得到pep文件
原来的文件格式

sed 's/^\(>[^ ]*\).*/\1/' Lindera_aggregata.gene.pep >Lagg.pep
sed -i 's/\*//g' Lagg.pep

要求的文件格式：

或者也可以用别的编程语言实现，如python

import re
with open("D:\yuceji\shili.pep", 'r') as file:
    lines = file.readlines()
output_file = 'output.fasta'

with open(output_file, 'w') as file:
    for i in range(len(lines)):
        line = lines[i].strip()
        if re.match(r'^>', line):
            gene_id = line.split()[0]
            seq = ''
            j = i + 1
            while j < len(lines) and not re.match(r'^>', lines[j].strip()):
                seq += lines[j].strip()
                j += 1
            file.write(gene_id + '\n')
            file.write(seq + '\n')

3. 使用ParaAT2.0 + KaKs_Calculator2.0

ParaAT.pl -h test.homologs -n test.cds -a test.pep -p proc -m muscle -f axt -g -k -o output 2> ParaAT.log &

-h test.homologs：指定包含同源序列的文件

-n test.cds：指定包含核苷酸序列的文件

-a test.pep：指定包含蛋白质序列的文件

-p proc：指定并行处理的进程数, (proc 可以自行设置）

-m muscle：指定用于多序列比对的软件，这里使用的是muscle

-f axt：指定输出文件的格式，这里使用的是axt格式

-g：启用gap stripping功能，即移除比对序列中的缺失片段。

-k：启用Ka/Ks比值的计算

-o output：指定输出文件的名称

2> ParaAT.log：将标准错误输出重定向到名为"ParaAT.log"的日志文件中。

&：将命令放入后台运行。

4. 使用TBtools软件

基于TBtools软件，可以使用" Sinple Ka/Ks Calculator" 程序

同样依然要准备那三个文件

从左到右每一列的介绍：

1. Seq_1：表示序列对中第一个序列的标识符或名称。

2. Seq_2：表示序列对中第二个序列的标识符或名称。

3. Ka：表示非同义突变（氨基酸替换）的数量，也称为Ka值。

4. Ks：表示同义突变（氨基酸保守替换）的数量，也称为Ks值。

5. Ka/Ks：表示Ka值除以Ks值得到的比值。（Ka/Ks比值用于衡量非同义突变和同义突变的相对丰度，从而推断基因或序列的选择压力）

6. EffectiveLen：表示序列对的有效长度，即用于计算Ka和Ks的比对序列的长度。

7. AverageS-sites：表示平均同义突变位点的数量。

8. AverageN-sites：这一列表示平均非同义突变位点的数量。

9. cN：这一列表示非同义突变位点的校正计数。

10. cS：这一列表示同义突变位点的校正计数。

11. pN：表示非同义突变位点的概率。

12. pS：这一列表示同义突变位点的概率。

其实得到的文件中除了这12列外，还有一列”Note",来记录一些“high sequence divergence value (ps>=0.75)” 这类的信息。

三、额外

总之，最主要的是得到了Ka/Ks的计算结果，达到了我们的目标，其余的也没必要全都了解。
另外，对于结果文件的每一列，我是这样记忆的：
对于Ka与Ks 的区分，可以Ks中的“s"来记忆，因为这可以很快与"same"这个单词联想起来。即相同的意思，故Ks就是同义突变。至于另一个Ka就是非同义突变了。
AverageS-sites 与AverageN-sites 的记忆也是如此。
cN 与 cS ，也是如此，另外"c"可表示"count",即计数
pN与pS，也是如此，另外"p" 可表示"probablity",即概率的意思

---

总结

本文主要讲述的是Ka/Ks 值的计算，其实无论是使用ParaAT2.0 + KaKs_Calculator2.0来计算，还是使用TBtools软件中的" Sinple Ka/Ks Calculator" 程序。都是为了得到结果。我个人还是比较推荐使用TBtools软件的，因为足够简单（不需要在Linux中再下载一些什么软件后再配置环境了，只要输入三个文件，就能得出结果文件）。当然或许它也会有一些其它问题。这样我们可以具体情况再具体分析。

无丝竹之乱耳，无案牍之劳形。

–2023-8-21 实验篇