Ka/Ks和4DTv计算可视化

作者: wo_monic | 来源:发表于2020-11-30 20:25 被阅读0次

参考地址
Ka/Ks和4DTv用于分析基因组的复制情况，可共线性分析。属于比较基因组学。
Ka:非同义替换
Ks:同义替换
Ka＞＞Ks，基因受正向选择positive selection，说明基因受强烈正向选择，是近期正在快速进化的基因。
Ka=Ks，基因中性进化neutral evolution
Ka＜＜Ks，基因受纯化选择purify selection
Ka/Ks: Ka/Ks > 1，则认为存在正选择效应（positive selection）；若Ka/Ks = 1，则认为存在中性选择效应；如若Ka/Ks < 1，则认为存在负选择效应，即纯化效应或净化选择（purifying selection）。
一般我们认为非同义突变会遭受自然选择作用。Ka比较大，说明非同义突变比较多，即该基因家族正在受到强烈的自然选择作用，意思是正在快速进化的基因。

4DTv：四倍兼并位点,简单说就是一个密码子，3个碱基中，有一个位置的碱基无论是A,T,C,G中的哪一个，都是编码同一个氨基酸，这个位点就是四倍兼并位点。例如："TCA/TCT/TCC/TCG"中的第三位就是一个4DTV，在这个位点的碱基，无论怎么突变最终都是丝氨酸(Ser)。

转载自https://www.jianshu.com/p/9d28de3d18e6

检测4DTv使用 MCscan

1.准备数据
Arabidopsis_thaliana.cds、Arabidopsis_thaliana.pep、Arabidopsis_thaliana.gff3；Citrus_grandis.cds、Citrus_grandis.pep、Citrus_grandis.gff3；Citrus_sinensis.cds、Citrus_sinensis.pep、Citrus_sinensis.gff3；Malus_domestica.cds、Malus_domestica.pep、Malus_domestica.gff3；

数据处理
awk '$3=="gene"' Arabidopsis_thaliana.gff3|wc -l 查看基因的数量
awk '$3=="mRNA"' Arabidopsis_thaliana.gff3|wc -l查看转录本的数量
无论是共线性分析或计算KaKs,使用的都是需要和基因数量相同的蛋白序列。如果你的蛋白序列的条数和转录本的数量一样，则说明你有冗余，仔细检测是否真的成功获取了最长的转录本
取file.pep中最长转录本序列，去冗余。一个基因可能有多个转录本，从中找出最长的转录本。
removeRedundantProteins.py使用时候需要注意你的ID，需要修改为分割字符是 .,基因名是第1各字段，则li = line2[0],如果基因名是第2个字段，则li = line2[1]
gn005G428000.1这种是默认的格式，直接使用脚本即可。
>gn.005G428000.1这种基因名是第2个字段，需要调整li = line2[1]
removeRedundantProteins.py的代码如下：

#!/usr/bin/env python3
import sys,getopt
def usage():
    print('usage:python3 removeRedundantProteins.py -i <in_fasta> -o <out_fasta> <-h>')
    return
def removeRedundant(in_file,out_file):
    gene_dic = {}
    flag = ''
    with open (in_file) as in_fasta:
        for line in in_fasta:
            if '>' in line:
                line1 = line.strip('>\n')
                line2 = line1.split('.')
                li = line2[0]
                flag = li
                try:
                    gene_dic[li]
                except KeyError:
                    gene_dic[li] = [line]
                else:
                    gene_dic[li].append(line)
            else:
                gene_dic[flag][-1] += line
    with open (out_file,'w') as out_fasta:
        for k,v in gene_dic.items():
            if len(v) == 1:
                out_fasta.write(gene_dic[k][0])
            else:
                trans_max = ''
                for trans in gene_dic[k]:
                    a = len(list(trans))
                    b = len(list(trans_max))
                    if a > b:
                        trans_max = trans
                out_fasta.write(trans_max)
def main(argv):
    try:
        opts, args = getopt.getopt(argv,'hi:o:')
    except getopt.GetoptError:
        usage()
        sys.exit()
    for opt, arg in opts:
        if opt == '-h':
            usage()
            sys.exit()
        elif opt == '-i':
            in_fasta_name = arg
        elif opt == '-o':
            outfile_name = arg
    try:
        removeRedundant(in_fasta_name,outfile_name)
    except UnboundLocalError:
        usage()
    return
if __name__ == '__main__':
    main(sys.argv[1:])

处理完成后，获取基因组蛋白质最长转录本。
3.获取共线性基因对
使用的是blast,本地使用conda安装。

conda install -c bioconda blast
# blast安装perl模块的方法
conda install perl-digest-md5

（1）对蛋白序列构建索引
makeblastdb -in Citrus_sinensis.pep -dbtype prot -parse_seqids -out Citrus_sinensis
（2）blastp比对
nohup blastp -query Citrus_sinensis.pep -out Citrus_sinensis.blast -db Citrus_sinensis -outfmt 6 -evalue 1e-10 -num_threads 8 -qcov_hsp_perc 50.0 -num_alignments 5 2> blastp.log &
（3）gff3文件处理
grep '\smRNA\s' Citrus_sinensis.gff3 | awk '{print $1"\t"$4"\t"$5"\t"$9}' | awk -F 'ID=' '{print $1$2}' | awk -F 'Parent=' '{print $1}' | awk '{print $1"\t"$4"\t"$2"\t"$3}' | awk -F ';' '{print $1"\t"$3}' | awk '{print $1"\t"$2"\t"$3"\t"$4}' > Citrus_sinensis.gff
（4）运行MCScanX(.blast & .gff文件 )

下载MCScanX,下载地址

使用方法 MCscan2 github

unzip MCScanX.zip
cd MCScanX
make

安装MCScanX,遇到3处错误和1个缺少Java的警告。解决方法前三次是，修改报错的文件的头部，增加一行#include <getopt.h>.最后一处错误是需要安装jdk.错误解决详情参考
JDK下载地址 https://www.oracle.com/technetwork/java/javase/downloads/index.html
下载linux版得到tar.gz，直接解压，进入目录，添加环境变量就可以使用。
运行MCScanX，获取共线性基因对

MCScanX Citrus_sinensis -g -3 -e 1e-10

得到 Citrus_sinensis.collinearity、Citrus_sinensis.tandem文件及Citrus_sinensis.html文件夹，其中我们需要的信息就在Citrus_sinensis.collinearity结果文件中
（5）提取共线性基因对
cat Citrus_sinensis.collinearity | grep "Cs" | awk '{print $3"\t"$4}' > Citrus_sinensis.homolog
4、Ka、Ks及4Dtv值计算
（1）准备软件和脚本
KaKs_Calculator2.0安装使用二进制文件，省心省力。
ParaAT官网
下载后解压缩，添加环境变量即可直接运行。
需要用到的脚本：calculate_4DTV_correction.pl；axt2one-line.py
https://github.com/JinfengChen/Scripts/blob/master/FFgenome/03.evolution/distance_kaks_4dtv/bin/calculate_4DTV_correction.pl
https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
calculate_4DTV_correction.pl内容如下

#!/usr/bin/perl 
use strict;

##author: sun ming'an, sunma@genomics.org.cn
##modifier: fanwei, fanw@genomics.org.cn
##correction: LiJun, junli@genomics.org.cn
##Date: 2008-9-24

##4dtv (transversion rate on 4-fold degenerated sites) are calculated with HKY substitution models 
##Reference: M. Hasegawa, H. Kishino, and T. Yano, J. Mol. Evol. 22 (2), 160 (1985)

die "perl $0 AXTfile > outfile\n" unless( @ARGV == 1);
# %codons was changed to DNA codon according to the poster below. Thanks to 小霖123.
# https://www.jianshu.com/p/6d704378c342
my %codons=(
'CTT'=>'L', 'CTC'=>'L', 'CTA'=>'L', 'CTG'=>'L',
'GTT'=>'V', 'GTC'=>'V', 'GTA'=>'V', 'GTG'=>'V',
'TCT'=>'S', 'TCC'=>'S', 'TCA'=>'S', 'TCG'=>'S',
'CCT'=>'P', 'CCC'=>'P', 'CCA'=>'P', 'CCG'=>'P',
'ACT'=>'T', 'ACC'=>'T', 'ACA'=>'T', 'ACG'=>'T',
'GCT'=>'A', 'GCC'=>'A', 'GCA'=>'A', 'GCG'=>'A',
'CGT'=>'R', 'CGC'=>'R', 'CGA'=>'R', 'CGG'=>'R',
'GGT'=>'G', 'GGC'=>'G', 'GGA'=>'G', 'GGG'=>'G'
);

my %transversion = (
        "A" => "TC",
        "C" => "AG",
        "G" => "TC",
        "T" => "AG",
);

my $axtFile = shift;

open(AXT,"$axtFile")||die"Cannot open $axtFile\n";
$/ = "\n\n";
my @seqs = <AXT>;
$/ ="\n";
close AXT;

print "tag\t4dtv_corrected\t4dtv_raw\tcondon_4d\tcodon_4dt\n";
foreach my $line ( @seqs ){
        chomp $line;
        if( $line =~ /^(\S+)\n(\S+)\n(\S+)$/ ){
                my $tag = $1;
                my $seq1 =$2;
                my $seq2 =$3;
                my ($corrected_4dtv, $raw_4dtv, $condon_4d, $codon_4dt) = &calculate_4dtv($seq1, $seq2);
                print "$tag\t$corrected_4dtv\t$raw_4dtv\t$condon_4d\t$codon_4dt\n";
        }
}



sub calculate_4dtv {
        my($str1, $str2) = @_;

        my ($condon_4d, $codon_4dt) = (0,0);
        my ($V,$a,$b,$d) = (0,0,0,0); 
        my %fre=();
        for( my $i = 0; $i < length($str1); $i += 3){
                my $codon1 = substr($str1, $i, 3);
                my $codon2 = substr($str2, $i, 3);
                my $base1= uc(substr($str1, $i+2, 1));
                my $base2= uc(substr($str2, $i+2, 1));
               
                if( exists $codons{$codon1} && exists $codons{$codon2} && $codons{$codon1} eq $codons{$codon2} ){
                    $fre{$base1}++;
                    $fre{$base2}++;
                    $condon_4d++;
                    $codon_4dt++ if(is_transversion($codon1,$codon2));
                }
        }
        
        if($condon_4d > 0){
            $V=$codon_4dt / $condon_4d; ##this is raw 4dtv value
            ##correction the raw 4dtv values by HKY substitution model
            $fre{"Y"}=$fre{"T"}+$fre{"C"};
            $fre{"R"}=$fre{"A"}+$fre{"G"};
            foreach (keys %fre){
                $fre{$_}=0.5*$fre{$_}/$condon_4d;
            }

            if($fre{Y}!=0 && $fre{R}!=0 && $fre{A}!=0 && $fre{C}!=0 && $fre{G}!=0 && $fre{T}!=0){
                $a=-1*log(1-$V*($fre{T}*$fre{C}*$fre{R}/$fre{Y}+$fre{A}*$fre{G}*$fre{Y}/$fre{R})/(2*($fre{T}*$fre{C}*$fre{R}+$fre{A}*$fre{G}*$fre{Y})));
                if (1-$V/(2*$fre{Y}*$fre{R}) > 0) {
                    $b=-1*log(1-$V/(2*$fre{Y}*$fre{R}));
                    $d=2*$a*($fre{T}*$fre{C}/$fre{Y}+$fre{A}*$fre{G}/$fre{R})-2*$b*($fre{T}*$fre{C}*$fre{R}/$fre{Y}+$fre{A}*$fre{G}*$fre{Y}/$fre{R}-$fre{Y}*$fre{R});
                }else{
                    $d = "NA";
                }
            }else{
                $d = "NA";
            }


        }else{
            $V="NA";
            $d="NA";
        }

        return ($d,$V,$condon_4d, $codon_4dt);

}


sub is_transversion{
        my ($codon1,$codon2) = @_;
        my $is_transversion = 0;
        my $base1 = substr($codon1,2,1);
        my $base2 = substr($codon2,2,1);
        $is_transversion = 1 if (exists $transversion{$base1} && $transversion{$base1} =~ /$base2/);
        return $is_transversion;
}

axt2one-line.py内容如下：

#!/usr/bin/python
#name: axt2one-line.py
import sys
import re
mydict = {}
args = sys.argv
reader = args[1]
writer = open(args[2],"w")
with open(reader) as fh:
    for line in fh:
        if re.search(r"\d+",line):
            header = line.strip()
            mydict[header] = []
        else:
            mydict[header].append(line.strip())
for key,value in mydict.items():
    seqs = "".join(value)
    length = len(seqs)
    mid = int(length/2)
    seq1 = seqs[0:mid]
    seq2 = seqs[mid:]
    writer.write(key + "\n" + "\n".join([seq1,seq2]))

（2）数据准备
Citrus_sinensis.homolog、Citrus_sinensis.cds、Citrus_sinensis.pep；Citrus_grandis.homolog、Citrus_grandis.cds、Citrus_grandis.pep；Arabidopsis_thaliana.homolog、Arabidopsis_thaliana.cds、Arabidopsis_thaliana.pep；Malus_domestica.homolog、Malus_domestica.cds、Malus_domestica.pep；
（3）使用ParaAT程序将蛋白序列比对结果转化为cds序列比对结果
新建proc文件, 设定12个线程

echo "12" >proc

-m参数指定多序列比对程序为clustalw2 (可以是其中任一个 clustalw2 | t_coffee | mafft | muscle)，-p参数指定多线程文件，-f参数指定输出文件格式为axt

ParaAT.pl -h Citrus_sinensis.homolog -n Citrus_sinensis.cds -a Citrus_sinensis.pep -m clustalw2 -p proc -f axt -o Citrus_sinensis_out 2> ParaAT.log

file.homolog、file.cds、file.pep三个文件的基因ID需保持一致，且file.cds及file.pep为fasta格式
当有两个物种时，homolog，cds,pep这三个文件，都是2个物种的。使用
cat Citrus_sinensis.homolog Arabidopsis_thaliana.homolog >file.homology合并成为一个。
当一个物种内部计算时，就使用这一个物种的即可。
(4）使用KaKs_Calculator计算ka、ks值, -m参数指定kaks值的计算方法为YN模型

for i in `ls *.axt`;do KaKs_Calculator -i $i -o ${i}.kaks -m YN;done
# 将多行axt文件转换成单行
for i in `ls *.axt`;do axt2one-line.py $i ${i}.one-line;done
#（5）使用calculate_4DTV_correction.pl脚本计算4dtv值
ls *.axt.one-line|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
#（6）合并所有同源基因对的4dtv
for i in `ls *.4dtv`;do awk 'NR>1{print $1"\t"$3}' $i >>all-4dtv.txt;done

#（7）合并所有同源基因对的kaks值
for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done

#（8）给结果文件进行排序和去冗余
sort all-4dtv.txt|uniq >all-4dtv.results
sort all-kaks.txt|uniq >all-kaks.results
#（9）将kaks结果和4Dtv结果合并
join -1 1 -2 1 all-kaks.results all-4dtv.results >all-results.txt
#（10）给结果文件添加标题
sed -i '1i\Seq\tKa\tKs\tKa/Ks\t4dtv_corrected' all-results.txt

四、可视化作图
1、数据准备
Arabidopsis_thaliana-4dtv.txt、Citrus_grandis-4dtv.txt、Citrus_sinensis-4dtv.txt、Malus_domestic-4dtv.txt
2、数据处理
cat citrus_sinensisl-4dtv.txt | awk '{print "C.sinensis""\t"$2}' > C.sinensis-4dtv.txt
合并A.thaliana-4dtv.txt、C.grandis-4dtv.txt、C.sinensis-4dtv.txt、M.domestic-4dtv.txt文件
cat A.thaliana-4dtv.txt C.grandis-4dtv.txt C.sinensis-4dtv.txt M.domestic-4dtv.txt > 4species-4dtv.txt

给4species-4dtv.txt 文件添加标题

sed -i '1i\Species\tValue' 4species-4dtv.txt

3、RStudio作图

setwd("C:/Users/***/Desktop/4Dtv/")
### 读入4species-4dtv.txt文件
dtv <- read.table("4species-4dtv.txt", header = T, check.names =F, sep = "\t")
### 载入R包
library(ggplot2)
library(hrbrthemes)
library(dplyr)
library(tidyr)
library(viridis)

### 绘图
p <- ggplot(data=dtv, aes(x=Value, group=Species, fill=Species)) + geom_density(adjust=1.5, alpha=.4) + theme_ipsum() + labs(title = "Distribution of 4Dtv distances", x = "4Dtv Value", y= "Density", size = 1.5)

上述脚本已经打包成流程 KK4D github。目前可视化模块已完成，共线性使用的是jcvi,后期可能会增加MCScanX的实现。其实自己把MCScanX的共线性结果格式化后，重命名为标准格式，也可以使用KK4D进行后续分析。

Ka/Ks和4DTv计算可视化

检测4DTv使用 MCscan

下载MCScanX,下载地址

给4species-4dtv.txt 文件添加标题

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

基因组

基因组学

注释和富集

基因组组装

orthofinder

基因组组装

系统进化