计算旁系同源基因的Ks

旁系同源基因简单来说就是同一物种内，通过基因复制产生的不同的两种基因。

计算直系同源的办法网上已经很详细，今天我来补充一下计算旁系同源基因ks的办法。

最简单方法1 用wgd 这个软件

conda create -n py36 python=3.6.7 blast mcl muscle mafft prank paml fasttree cmake libpng mpi=1.0=mpich
conda activate py36
git clone https://github.com/arzwa/wgd.git
cd wgd
pip install .
#可能要下载i-ADHoRe-3.0  http://bioinformatics.psb.ugent.be/webtools/i-adhore/licensing/
#wgd mcl 生成.mcl文件  鉴定基因组内的同源基因
#只需要一个物种的cds就可以了 计算方法是Yn00 
#以毛果杨4.1版本为例
nohup wgd mcl -n 20 --cds --mcl -s Ptrichocarpa_533_v4.1.cds.fa  -o Ptrichocarpa_533.out  &
nohup wgd ksd  Ptrichocarpa_533_v4.1.cds.fa.blast.tsv.mcl  Ptrichocarpa_533_v4.1.cds.fa  -n 10  &
#wgd mcl 生成的blast.tsv.mcl放在和cds.fa文件一个文件夹

这里建议取出ks值用R画图
用awk取出第一列和第九列

cat populus_simonii_CDS.fa.ks.tsv | awk '{print $1"\t"$9}' > p_siminii.txt
data01<-read.table("simks.txt",header = T,sep="\t")
names(data01)<-c("Species","Ks")
data01[,1] <- "P.sim"
p1<-ggplot(data = data01,aes(x =Ks,fill=Species))+geom_histogram(position = "identity",bins=250,alpha=0.15,aes(y = ..density..))+xlim(0,3)+ylim(0,2)+stat_density(geom = "line",position = "identity",aes(colour = Species))+theme_bw()+
  geom_vline(xintercept = Peak1[[1]],colour="red",linetype="dashed") +
  geom_text(aes(x=Peak1[[1]], y= Peak1[[2]]+0.3,label=paste("Ks =",round(Peak1[[1]],4))),color="black") +
  theme(axis.title = element_text(size=8),axis.text=element_text(size=8),legend.position =  "top")
p1
densFindPeak <- function(x){
  td <- density(x)
  maxDens <- which.max(td$y)
  list(x=td$x[maxDens],y=td$y[maxDens])
}
S1 <- data01[data01$Species=="P.sim",]
S1_limit <- S1$Ks[S1$Ks>=0.005 & S1$Ks<=3]
S1_limit = na.omit(S1_limit)
Peak1 = densFindPeak(S1_limit)

第二种方法 KaksCalculator的办法利用MCScanX提取同源基因对
这里还是用从
https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Ptrichocarpa_er 下载数据

#需要下载 clustalw KaKs_Calculator2.0 ParaAT2.0 blast MCScanX
#先处理蛋白质数据只留下序列名id
python chang.py  Ptrichocarpa_533_v4.1.protein.fa >ptr.pep
 makeblastdb -in ptr.pep -dbtype prot
nohup blastp -query  ptr.pep -db ptr.pep -evalue 1e-10  -num_threads 20 -out ptr_ptr_blastp_out.m6 -outfmt 6 &
mv ptr_ptr_blastp_out.m6 ptr.blast
#处理gff
grep -v "#" Ptrichocarpa_533_v4.1.gene.gff3 >ptr.gff3
awk '$3 == "mRNA" {print $1 "\t" $9 "\t" $4 "\t" $5}' ptr.gff3 | sed 's/.v4.1/\t/g' |sed 's/ID=//g' | awk '{print $1 "\t" $2 "\t" $4  "\t" $5}'> ptr. gff
MCScanX ptr
 echo "20" > proc
grep -v '#' ptr.collinearity | awk '{print $3 "\t" $4}' > ptr.homolog
ParaAT.pl -h ptr.homolog -n ptr.fa  -a ptr.pep -m clustalw2 -p proc -f axt -o ptr_ptr
cd ptr_ptr/
for i in `ls *.axt`;do KaKs_Calculator -i $i -o ${i}.kaks -m YN;done
for i in `ls *.kaks`;do awk 'NR>1{print $1 "\t" $4}' $i >>ptrkaks.txt;done

计算4DTv值用下面perl脚本从githup修改来的

#!/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);

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',
'CCU'=>'P', 'CCC'=>'P', 'CCA'=>'P', 'CCG'=>'P','CCT'=>'P',
'ACU'=>'T', 'ACC'=>'T', 'ACA'=>'T', 'ACG'=>'T','ACT'=>'T',
'GCT'=>'A', 'GCC'=>'A', 'GCA'=>'A', 'GCG'=>'A',
'CGT'=>'R', 'CGC'=>'R', 'CGA'=>'R', 'CGG'=>'R',
'GGU'=>'G', 'GGC'=>'G', 'GGA'=>'G', 'GGG'=>'G','GGT'=>'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;
}