CentOS7系统使用小经验

发表于2020 年 1 月 28 日由chenlianfu

1. 增大系统字体

若屏幕分辨率为4k，则系统画面字体很小，眼睛看着难受。可以：（1）在gnome桌面点击右上角，再点击扳手图标打开控制面板；（2）再在控制面板中点击Devices，在Display中选择200%的Scale，则可以将整体画面都放大一倍，则感觉正常。

若在windows系统中使用虚拟机进行操作时，重新全屏一次或分辨率重新调节一次后，画面又恢复到正常的小字体情况，较为不便。则采用字体放大方法：（1）依次点击Applications，Accessories，Tweaks，打开Tweaks软件；（2）在左侧栏点击Fonts，将Scaling Factor设置为2，即可将所有字体放大一倍，则感觉正常。

2. 常用系统快捷键

windows + 左箭头：将选中的窗口放大到左半屏（左半侧50%全屏）；
windows + 右箭头：将选中的窗口放到到右半屏（右半侧50%全屏）；
windows + PgUp：切换到上一个桌面；
windows + PgDn：切换到下一个桌面；
windows + Tab：切换到同桌面中的下一个窗口；
windows + `：切换到下一个相同软件的窗口；
PrtScr：全屏截屏，保存到家目录；
alt + PrtScr：对当前窗口截屏，保存到家目录；
ctrl + PrtScr：对当前窗口截屏，保存到内存以供粘贴；
shift + PrtScr：对指定区域截屏，保持到家目录。

3. 手动设置系统快捷键

打开控制面板，依次点击Devices，Keyboard，可以查看并修改各快捷键设置。我个人喜欢的设置：

windows + d：隐藏所有窗口；
ctrl + 1：切换到第一个窗口（使用左边的ctrl键和小键盘数字）；
ctrl + 2：切换到第二个窗口；
ctrl + 3：切换到第三个窗口；
ctrl + 0：切换到第四个窗口；
alt + d：打开gnome-terminal终端软件（和ctrl + d退出较类似）。

4. 自动更新系统时间

在windows系统下使用CentOS虚拟机系统。由于经常休眠windows系统导致CentOS虚拟机中时间经常不正确，需要更新CentOS系统时间。使用root用户执行如下操作，能自动更新系统时间。

echo -e "*/5\t*\t*\t*\t*\t/sbin/ntpdate cn.pool.ntp.org" > /root/.crontab
crontab /root/.crontab

5. 一次性打开多个排列整齐的终端，并同时连接服务器

我个人喜欢一次性打开3个终端：屏幕左半屏（左侧50%面积区域）上下各一个终端，屏幕右半屏（右侧50%面积区域）一个终端。手工生成一个gnome-terminal配置文件：

[GNOME Terminal Configuration]
Version=1
CompatVersion=1
Windows=Window0;Window1;Window2;

[Window0]
Geometry=106x27+2+60
Terminals=Terminals0
[Terminals0]
Command=dell

[Window1]
Geometry=105x56+1925+60
Terminals=Terminal1
[Terminal1]
Command=dell

[Window2]
Geometry=106x28+2+1055
Terminals=Terminal2
[Terminal2]
Command=dell

使用如下命令即可一次性打开3个终端，并自动使用dell命令连接到戴尔服务器：

gnome-terminal --load-config=/home/chenlianfu/.terminal_open3

CentOS7安装mplayer

发表于2020 年 1 月 25 日由chenlianfu

sudo su -
yum -y install epel-release
rpm -ivh http://li.nux.ro/download/nux/dextop/el7/x86_64/nux-dextop-release-0-5.el7.nux.noarch.rpm
yum -y install mplayer

配置DELL MD3200存储服务器连接到多台主机

发表于2019 年 12 月 11 日由chenlianfu

戴尔MD3200磁盘存储服务器仅支持DELL SAS硬盘。该存储服务器可以上12块3.5英寸SAS硬盘、2个电源和2个控制器。每个控制器含有一个阵列卡带2G缓存，支持4个SAS IN接口用于连接4台主机。每台MD3200机器上2个控制器，则可以同时连接8台主机。此外，每个控制器含有一个SAS OUT口，用于连接到MD1200磁盘柜的SAS IN接口，从而支持拓展多达96块SAS硬盘。

1. 硬件组装

我本次经验，是将12块4TB的DELL SAS盘插入到DELL MD3200存储服务器中。然后，分别在3台服务器主机的PCI-E插槽上分别插上一片SAS直通卡。每片SAS直通卡有两个SF8088的SAS接口。使用6根SAS线（SF8088-SF8088）将MD3200存储上的SAS IN口和3台服务器的3片SAS直通卡进行相连。

此外，需要注意每张SAS直通卡贴纸上的一个编号，推荐拍照留存，以备后续使用。

2. 对DELL MD3200进行磁盘阵列设置

在个人windows10系统的笔记本电脑上安装软件DELL_MDSS软件，并设置笔记本电脑的有线网卡IP地址为192.168.128.120，使之和MD3200控制器IP地址192.168.128.101的IP段一致。然后将笔记本电脑的有线网口和控制器的网卡用网线连接，再打开DELL的MDSM软件。若两者有线网卡在同一IP地址段，则能自动找到存储阵列。

推荐重命名自动找到的存储阵列，然后双击对其进行管理，启动阵列管理窗口。再先后进行两项设置（1）阵列设置；（2）主机映射。

2.1 阵列设置

首先将最后两块盘设置为热备盘，然后以RAID5方式将其余所有的10块盘组成阵列，并使用所有的空间创建一个分区。

2.2 主机映射

在保证DELL MD3200存储服务器和3台服务器主机都开机并连接好数据线后，重启所有的服务器主机，以利于MD3200存储能识别3片SAS直通卡。

在MDSM的整列卡管理窗口中，创建一个名为chenlianfu的主机组，然后在该主机组下依次添加6个主机名。每添加一个主机时，需要设置一个SAS直通卡的机器码（需要和SAS直通卡上的机器码对应）和别名。每片SAS直通卡有两个SAS口，有两个机器码，对应添加的2个主机名。

主机映射做好后，即允许6个指定的SAS口可以访问磁盘存储。然后重启3台服务器，此时在重启界面中，在检测到SAS直通卡的时候，应该会提示其SAS直通卡检测到虚拟磁盘。在进入系统后，则会增加 /dev/sdb 和 /dev/sdc 两个设备，分别对应SAS直通卡的两个SAS口。

3. 对服务器主机进行多路径设置

在/dev/sdb和/dev/sdc两个路径都同时指向了存储服务器。需要通过CentOS系统中的多路径方式对这两个路径进行设置，使数据能同时通过这两个路径进行传输，才能充分达到6Gbps的数据。。

yum install https://dl.fedoraproject.org/pub/epel/epel-release-latest-7.noarch.rpm
yum install *device-mapper-multipath*
yum install *kpartx*
yum install *scsi_dh_rdac*
yum install *scsi_dh*
yum install *scsi_id*
yum install *udev*
yum install *modutils*
yum install *iscsi-initiator-utils*

在三台服务器上使用root用户都进行如下操作：

multipath -v3

使用如上命令得到/dev/sdb的wwid编号，然后，生成/etc/multipath.conf文件，其内容为：

blacklist {
    devnode "^sda"
# Begin Dell MD Modification
	device {
		vendor  "*"
		product "Universal Xport"
	}
	device {
		vendor  "*"
		product "MD3000"
	}
	device {
		vendor  "*"
		product "MD3000i"
	}
	device {
		vendor  "*"
		product "Virtual Disk"
	}
# End Dell MD Modification
}
defaults {
    user_friendly_names yes
    path_grouping_policy multibus
    failback immediate
    no_path_retry fail
# Begin Dell MD Modification
	max_fds             8192
# End Dell MD Modification
}
multipaths {
    multipath {
        wwid 36d4ae52000ad3c1800000c755df2d836
        alias chenlianfu_md3200
        path_grouping_policy multibus
        path_checker tur
        path_selector "round-robin 0"
    }
}
# Begin Dell MD Modification
devices {
	device {
		vendor                "DELL"
		product               "MD32xxi"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD32xx"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD36xxi"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD36xxf"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD34xx"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD38xxf"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
	device {
		vendor                "DELL"
		product               "MD38xxi"
		path_grouping_policy  multibus
		prio                  rdac
		polling_interval      5
		path_checker          rdac
		path_selector         "round-robin 0"
		hardware_handler      "1 rdac"
		failback              immediate
		features              "2 pg_init_retries 50"
		no_path_retry         30
		rr_min_io             100
		prio_callout          "/sbin/mpath_prio_rdac /dev/%n"
	}
}
# End Dell MD Modification

再启动multipath服务：

systemctl restart multipathd.service 
systemctl enable multipathd.service

这时，就可以在 /dev/mapper/目录下看到相应的设备了，在进行分区挂载操作等：

mkdir /disks/md3200/
mount /dev/mapper/chenlianfu_md3200p1 /disks/md3200/

最后，可以在 /etc/fstab 中增加一行，进行开机自动挂载：

/dev/mapper/chenlianfu_md3200p1 /disks/md3200   xfs     defaults        0       0

使用nginx运行wordpress

发表于2019 年 11 月 28 日由chenlianfu

apache和nginx都能搭建网页服务。apache更老牌，功能模块更多，运行更稳定；而nginx运行性能更高，内存消耗更少。使用apache时，每一个连接请求都对应服务器上一个进程，而在nginx中，多个连接(多达上万个)仅对应服务器上一个进程。相比于apache支持不超过3000个连接请求，nginx能支持高并发连接，每秒最多的并发连接请求理论可以道道5万个。本问讲解在不影响apache情况下额外使用nginx运行wordpress网站。

1. 安装wordpress软件

目前wordpress官网不再对中国内陆提供软件的下载了。需要使用代理下载wordpress软件。此外，最新版本的wordpress需要较高版本的PHP软件，而CentOS6或CentOS7自带的PHP软件版本较低，推荐使用较低版本的wordpress软件，例如：5.1.3版本。

cd /home/chenlianfu
wget https://cn.wordpress.org/latest-zh_CN.tar.gz
tar zxf latest-zh_CN.tar.gz

注意不管是使用apache还是nginx运行wordpress，最终其实都是使用apache用户来对软件目录进行读取或写入。需要使用root用户修改权
相关权限：

usermod -aG chenlianfu apache
chmod 750 /home/chenlianfu
chown -R apache:chenlianfu /home/chenlianfu/wordpress/

2. 安装nginx和php-fpm软件

使用root权限安装nginx和php-frm软件，后者用于让nginx识别php文件。此外，在安装CentOS系统时使用最大化安装，默认安装上了相应的一些mysql和php软件等。

yum install nginx php-fpm

3. 启动php-fpm软件

/etc/init.d/php-fpm restart
chkconfig php-fpm on

需了解的时是，在php-fpm的配置文件/etc/php-fpm.d/www.conf中，默认指定的使用者依然是apache。

4. 修改nginx配置文件

在nginx的主配置文件/etc/nginx/nginx.conf中修改端口号为8088。在CentOS6中可能是修改配置文件/etc/nginx/conf.d/default.conf。

    listen 8088;
    #listen       80 default_server;
    #listen       [::]:80 default_server;

再增加nginx对wordpress的配置文件/etc/nginx/conf.d/wordpress_chenlianfu.conf，其内容为：

server {
    listen 8088;
    server_name www.chenlianfu.com;
    root /home/chenlianfu/wordpress;
    access_log  /var/log/nginx/host.access.log  main;
    location / { 
        try_files $uri $uri/ /index.php?$args;  
        index  index.php;
    }   
    
    location ~ \.php$ {
        fastcgi_pass 127.0.0.1:9000;
        fastcgi_index index.php;
        include fastcgi_params;
        fastcgi_param SCRIPT_FILENAME $document_root$fastcgi_script_name;
    } 
}

最后，启动php-fpm和nginx服务并设置开机自动启动：

/etc/init.d/nginx restart
chkconfig nginx on

此时，可以在网页上使用8088端口访问nginx运行的wordpress网站了。这需要在浏览器中输入www.chenlianfu.com:8088来访问网站（这使还需要购买该域名并指向服务器的IP地址），比较不方便。

5. 修改apache配置文件实现apache和nginx共用80端口

默认情况下nginx和apache都使用了80端口，形成了冲突。我一般优先使用apache构建网站，有利于对服务器中的文件夹在网页中以目录形式访问。nginx对目录访问的效果较差，访问时其网址最后面一定要加一个目录符号正斜线 / ，否则不能访问其目录。因此，我优先让apache占用80端口，而让nginx占用8088端口。从而两种建站方式同时支持。

为了让nginx的网站直接使用80端口访问，这时，只需要在apache配置文件中设置将对www.chenlianfu.com:80的访问自动转到www.chenlianfu.com:8088即可。在/etc/httpd/conf/httpd.conf最尾部添加信息：

<VirtualHost *:80>
    ProxyPreserveHost On
    ServerName www.chenlianfu.com
    ProxyPass / http://127.0.0.1:8088/
    ProxyPassReverse / htt://127.0.0.1:8088/
</VirtualHost>

以上配置信息表示，每当访问www.chenlianfu.com的80端口时，自动转到本地机器的8088端口，即使用nginx提供的网页服务了。从而实现apache和nginx共用80端口的效果了。然后重启apache服务：

/etc/init.d/httpd restart

若敢兴趣的话，其实也可以优先让nginx占用80端口，让apache软件占用8088端口，然后在nginx的配置文件中设置端口转发，实现两个软件对80端口的共用。修改配置文件/etc/nginx/nginx.conf实现端口转发：

location ~ /.*/ {
    proxy_pass http://127.0.0.1:8088;
}

将以上信息添加到server { }块中。location能以这则表达式的方法对网址进行解析。上述 ~ / / 表示正则匹配，和perl类似；.* 表示匹配任意内容。因此，上述内容表示将任意网址都转到8088端口上。

搭建邮箱服务并批量化做多个邮箱帐号用于KAAS注释

发表于2019 年 11 月 27 日由chenlianfu

在使用KAAS网页工具进行注释时，一个邮箱同一时刻仅允许计算一个任务。若需要同时对许多物种进行KAAS分析时，则需要多个邮箱帐号。

若在各公共平台上申请邮箱会比较麻烦。此外，每使用一个邮箱帐号去进行KAAS分析时，然后登录一个个邮箱确认提交任务和查看结果是很麻烦的。因此，我在我的Aliyun服务器上直接创建很多邮箱帐号，然后用这些帐号进行KAAS分析，KAAS返回的邮件统一转到我的个人帐号，再进行后续的提交和结果查询。

1. 在万网上购买域名并进行域名解析

我购买的域名是chenlianfu.com，并将该域名解析中的A指向我的Aliyun服务器的IP地址，将MX指向chenlianfu.com。从而能让邮件识别我的域名并发送到我的Aliyun服务器上。

2. 在Aliyun服务器上通过postfix启动邮箱服务

默认设置下，CentOS6系统自带postfix软件。首先，修改postfix的配置文件/etc/postfix/main.cf内容，需要修改的几行如下：

myhostname = chenlianfu.com  # 用于设置主机名，表示邮箱账户为linux系统用户后加 @chenlianfu.com 。
inet_interfaces = all  # 能接受所有来源的邮件
#inet_interfaces = localhost  # 关闭默认的设置，默认仅接收来源本机的邮件
mydestination = $myhostname, localhost.$mydomain, localhost, chenlianfu.com

message_size_limit = 100000000
mailbox_size_limit = 200000000 # 在centos8系统中要添加该参数值比message_size_limit值大，否则postfix后台运行状态不正常，导致用户没法查看邮件信息。

然后，防火墙开放25号端口：

iptables -A INPUT -p tcp -m state --state NEW -m tcp --dport 25 -j ACCEPT
service iptables save
/etc/init.d/iptables restart

最后启动postfix服务，则CentOS6系统中的用户可以接收来自其余人发送的邮件了。

/etc/init.d/postfix restart

若需要查看或寄送邮件，则输入命令 mail 即可：

mail

3. 在Aliyun服务器上创建多个系统用户

使用root权限创建50个账户：

perl -e 'foreach (1..50) { $_ = '0' . $_ if length($_) == 1; print "um$_\n"; } ' > users.txt
for i in `cat users.txt`
do
    useradd $i
    echo '123456' | passwd --stdin $i
done

为了简单演示，以上创建的用户密码都为123456，这里不太推荐使用如此简单的密码。

同时，需要设置当邮件发送到这些用户时，会自动转发到我的常用帐号所对应的邮箱中。我在Aliyun服务器上的账号名称是chenlianfu，对应的邮箱名称为chenlianfu@chenlianfu.com。通过使用root权限修改配置文件/etc/aliases文件来实现邮件转发。

for i in `cat users.txt`
do
    echo "$i:    chenlianfu,$i"
done >> /etc/aliases

为了使配置文件/etc/aliases生效，重启postfix即可：

/etc/init.d/postfix restart

4. 现在，即可使用这50个账号了。

封禁一些IPs的SSH登录，保护服务器安全

发表于2019 年 11 月 20 日由chenlianfu

登录服务器的root用户时，可能会收到如下提示：

There were 3225 failed login attempts since the last successful login.

这表示有很黑客在尝试用一些用户和密码登录服务器。为了保障服务器安全，则可以在日志文件/var/log/secure中查看攻击的IPs来源和其失败登录的次数。例如：

Invalid user neisius from 211.144.12.75 port 34041
Failed password for invalid user server from 219.94.99.133 port 15139 ssh2
pam_unix(sshd:auth): authentication failure; logname= uid=0 euid=0 tty=ssh ruser= rhost=106.12.217.39

若某个IP次数很高，则将其IP信息添加到/etc/hosts.deny文件中进行封禁。例如，修改/etc/hosts.deny文件内容：

sshd:101.109.83.140
sshd:103.27.238.202
sshd:104.236.94.202
sshd:106.12.200.13
sshd:106.12.217.39
sshd:106.54.19.67
sshd:111.230.73.133
sshd:117.50.45.254
sshd:118.24.38.53
sshd:122.155.223.119

为了能自动地封禁IPs，编写名为IP_denied_for_ssh_security.pl的perl程序：

#! /usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    $0 <num>

    运行程序后，程序读取/var/log/secure文件中信息。若存在某个IP登录失败(含有Invalid user|Failed password|authentication failure等信息)次数超过指定的次数，则是需要封禁的IP。此外，程序再读取/etc/hosts.deny文件中的IP信息，得到所有需要封禁的IP。然后程序将这些IP信息全部写入到/etc/hosts.deny文件中进行封禁。
    程序需要以root权限运行，因为只有root权限才能读取/var/log/secure文件，对/etc/hosts.deny进行写入。
    程序运行完毕后，会清空/var/log/secure文件内容。

USAGE
if(@ARGV==0){die $usage}

open IN, "/var/log/secure" or die "Can not open file /var/log/secure, $!";
my %num;
while (<IN>) {
    if (m/Invalid user/ or m/Failed password/ or m/authentication failure/) {
        $num{$1} ++ if m/(\d+\.\d+\.\d+\.\d+)/;
    }
}
close IN;

my %deny;
my $deny_num = 0;
foreach (keys %num) {
    if ($num{$_} >= $ARGV[0]) {
        $deny{$_} = 1;
        $deny_num ++;
    }
}
print STDERR "$deny_num IPs were detected for blocking from file /var/log/secure.\n";

open IN, "/etc/hosts.deny" or die "Can not open file /etc/hosts.deny, $!";
while (<IN>) {
    $deny{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
};
close IN;

open OUT, ">", "/etc/hosts.deny" or die "Can not create file /etc/hosts.deny, $!";
my $total_num = 0;
foreach (sort keys %deny) {
    print OUT "sshd:$_\n";
    $total_num ++;
}
close OUT;
print STDERR "Total $total_num IPs were banned by file /etc/hosts.deny for SSH security.\n";

open OUT, ">", "/var/log/secure" or die "Can not  create file /var/log/secure, $!";
close OUT;

使用root用户执行该perl程序，若登录失败次数>=3次，则封禁对应的IP。

IP_denied_for_ssh_security.pl 3

基因正选择分析原理

发表于2019 年 11 月 18 日由chenlianfu

一、正选择分析的目的：

两两基因的密码子序列进行比较，从而计算dN/dS，即omega（ω）值。若该值<1，则表示纯化选择；omega = 1，则中性进化；omega > 1，则正选择。若分析基因在两个物种中的序列，可以计算dN/dS的值，若omega > 1，即表明该基因在物种进化过程中，即由其祖先物种分化成这两个物种时，基因受到了正选择。对于两个物种/序列的正选择分析，比较简单。而实际情况中，要分析的物种数量很多，包含多个类群。这个时候的正选择分析相对复杂些。
对多个物种的基因序列进行正选择分析，若仍然按照两个物种时的要求，即分析该基因在物种进化中是否受到了正选择？这种结果可能不好说清楚。因为该基因可能在某一类群中序列很相似，其两两比较时，omega <= 1；而在另外一类群中两两比较时，很多时候omega > 1。最后软件可以从总体上给一个omega值，该值不可以拿来简单地评价该基因是否受到了正选择。所以，对多个物种进行正选择分析时，没法直接评价该基因是否受到了正选择。正选择只有在进行两两序列比较的时候，才能计算omega值，从而得到结果。
对基因在多个物种上的正选择分析，分析的目的则是：比较某个分枝上祖先节点和后裔节点（可以理解成，对无根树上某分枝两侧的两组物种进行比较，依然属于两两比较），从而计算该分枝的omega值。而在实际数据中，基因在不同的进化分枝上具有不同的omega值，同时在序列不同的位点也具有不同的omega值。目标分枝两侧的物种数量较多时，可以对序列上的每个位点进行omega值分析，从而鉴定正选择位点。所以，对基因在多个物种上的正选择分析，需要同时分析分析目标分枝的omega值和序列位点的omega值，从而判断基因是否受到正选择压。

二、使用PAM对基因进行正选择分析，有三种方法：

PAML site model: 主要用于检测基因中的正选择位点。该方法分析时，认为进化树中各分枝的omega值是一致的，并比较两种模型：(1)模型m1是null model，认为所有位点的omega值<1或=1; (2)模型m2是正选择模型，存在omega <1、=1或> 1的位点。比较两个模型的似然值（lnL）差异，利用卡方检验（自由度为2）算出p值。若p值 < 0.05，则否定null model，认为存在正选择位点。此外，推荐采用比较模型m7和m8，它们将omega值分成了10类，其p值结果比上一种比较方法更宽松，能检测到更多的正选择基因。使用PAML site model方法能在整体水平上检测基因的正选择位点，而不能表明基因在某个进化分枝上是否受到正选择压。
PAML branch-site model: 主要用于检测基因在某个进化枝上是否存在的正选择位点。该分析方法认为目标分化枝具有一个omega值，其它所有分枝具有一个相同的omega值，然后再检测正选择位点。同样对两种模型进行比较：（1）第一种模型为模型2，将omega值分成<1、=1、>1的三类，这和site model中的一样；（2）第二种模型和前者一致，只是将omega固定成1，作为null model。比较两种模型的似然差异，利用卡方检验（自由度为2）算p值（chi2命令算出的值除以2）。若p值< 0.05，则能通过Bayes Empirical Bayes (BEB)方法计算正选择位点的后验概率，若存在概率值 > 0.95正选择位点，则表示基因在目标分枝上受到正选择压。PAML软件在branch-site模式下，并不给出分枝上的omega值。这表示branch-site模式虽然考虑了目标分枝上具有不同的omega值，但仍然以分析位点上的omega为主。值得注意的是，在branch-site模式下可能检测到正选择位点，但在目标分枝上的omega值仍然可能低于1。可能软件作者基于这点考虑，就没有给出目标分枝上的omega值，以免影响一些人对正选择结果的判断。
PAML branch model: 主要用于检测在某个分枝上，其omega值是否显著高于背景分枝，即基因在目标分枝上进化速度加快。该方法认为基因序列上所有位点的omega值是一致的，对两种模型进行比较：（1）第一种模型为null model，所有分枝具有相同的omega值；（2）第二种模型认为目标分枝具有一个omega值，其它所有分枝具有一个相同的omega值。比较两种模型的似然差异，利用卡方检验（自由度为1）算p值。若p值 <= 0.05，且目标分枝上的omega值高于背景值，则认为该基因为快速进化基因。一般情况下，该方法计算得到的p值会低于第二种方法的结果。

三、其它注意事项

Branch-site model相比于site model的优点是考虑了不同的分枝具有不同的选择压，即具有不同的omega值。该方法让目标分枝具有一个不同的omega值，并没有让所有分枝的omega值独立进行计算（理论上这样是最好的）。这样算法很复杂，程序运行非常非常消耗时间。但其实也没必要这样做，因为正选择分析其实是两条序列比较后，分析dN/dS，再找正选择位点，其分析结果就应该是某个分枝上基因是否受到正选择，在序列那个位点上受到正选择。

若在目标分枝上，其omega值小于1，但是却能找到正选择位点。即该基因在该分枝上的dN/dS < 1，但是在某些位点上，dN/dS > 1。那么该基因是否属于正选择基因？我认为：属于。之所以为正选择基因，主要是因为基因的个别位点或多个位点存在正选择。当只有个别位点受到正选择压时，而其它多个位点存在纯化选择时，可能导致整体上的omega值小于1。此时，该基因也应该是属于正选择基因。

四、参考文献中的正选择分析方法描述

Science文章（https://science.sciencemag.org/content/364/6446/eaav6202）中的正选择基因分析方法：To estimate the lineage-specific evolutionary rate for each branch, the Codeml program in the PAML package (version 4.8) (134) with the free-ratio model (model = 1) was run for each ortholog. Positive selection signals on genes along specific lineages were detected using the optimized branch-site model following the author’s recommendation. A likelihood ratio test (LRT) was conducted to compare a model that allowed sites to be under positive selection on the foreground branch with the null model in which sites could evolve either neutrally and under purifying selection.[ The p values were computed based on Chi-square statistics, and genes with p value less than 0.05 were treated as candidates that underwent positive selection. We identified PSGs at the ancestral branch of Ruminantia (table S22), the ancestral branch of Pecora (table S23), each ancestral family branch of Ruminantia (tableS24), and each ancestral subfamily branch of Bovidae (table S24). We also compared the dN/dS values of Ruminantia families with outgroup mammals (fig. S52).

Science文章（https://science.sciencemag.org/content/364/6446/eaav6202）中快速进化基因分析方法：The branch model in PAML was used, with the null model (model=0) assuming that all branches have been evolving at the same rate and the alternative model (model=2), allowing the foreground branch to evolve under a different rate. An LRT with df =1 was used to discriminate between alternative models for each ortholog in the gene set. Genes with a p value less than 0.05 and a higher ω value for the foreground than the background branches were considered as evolving with a significantly faster rate in the foreground branch.

https://journals.plos.org/plosntds/article?id=10.1371%2Fjournal.pntd.0007463文献中对正选择基因的分析方法：

A calculation of mutational rate ratio ω between two gene sequences was the basis for the positive selection analysis. The ω was calculated as a ratio of nonsynonymous to synonymous mutational rates. The ratio indicates negative purifying selection (0 < ω < 1), neutral evolution (ω = 1), and positive selection (ω > 1) [54]. A set of selected genes from complete genomes was tested relative to positive selection using the maximum likelihood method using the CODEML of the PAML software package [55]. PAML version 4 [56] and its user interface PAMLX [57] were used in our study. For each analyzed gene, its maximum likelihood phylogenetic tree was used as an input tree. The CODEML offers several different codon evolutionary models, and the statistical likelihood ratio test (LRT) was used to compare the codon evolutionary model to the null model. The Bayes empirical Bayes method (BEB) [58] was then used to evaluate the posterior probability of sites considered to have been positively selected.

The CODEML models could produce different results (i.e., a list of sites under positive selection) since they calculate different parameter estimates. Site models allow ω to vary in each site (codon) within the gene. Statistical testing was required for sites with ω > 1. Two pairs of models were predominantly used since their LRTs have low false-positive rates. M1a (nearly neutral evolution) was compared to M2a (positive selection) [58,59] and M7 (beta) was compared to M8 (beta & ω) [60]. Our preliminary testing found that the two model pairs gave the same or very similar results. Therefore we chose to use the M7-M8 model pair. The M7 model is a null model that allows 10 classes of sites with a ω beta-distribution within the interval 0 ≤ ω ≤ 1. Sites with ω > 1 are not allowed. The alternative M8 model adds an eleventh class of sites with ω > 1. Each site was tested to determine the class to which it belongs. The LRT compares twice the log-likelihood difference 2Δl = 2(l1-l0) between the M7 model (log likelihood value l0) and the M8 model (log likelihood value l1) to the χ2 distribution [61]. If the twice log-likelihood difference is above a critical χ2 value, then the null model is rejected, and the positive selection is statistically significant.

A considerable disadvantage of the site models is that ω was calculated as an average over all codons of the site. Therefore, the site models are not suitable for the data where ω also varies between lineages. In contrast, the branch-site models search for positive selection in sites and pre-specifies lineages where different rates of ω may occur [62]. Sequences of lineages are a priori divided into a group of foreground lineages where positive selection may occur and group of background lineages where only purifying selection or neutral evolution occurs. We used branch-site model A, which allows four classes of sites and different setups of foreground lineages to be tested depending on the gene phylogeny. In branch-site model A, all lineages under purifying selection with a low value of ω0 belong to site class 0. Weak purifying selection and neutral evolution with ω1 near to value 1 are allowed in site class 1. In site class 2a, a proportion of class 0 sites in foreground lineages is under positive selection with ω2 > 1. Similarly, site class 2b is a proportion of class 1 sites under positive selection with ω2 > 1. The null model for LRT has ω2 = 1. Critical values of LRT (2Δl) are 2.71 at 5% and 5.41 at 1% [63]. The posterior probabilities of suggested sites under positive selection were calculated using the BEB method.

使用squid搭建代理服务器

发表于2019 年 11 月 16 日由chenlianfu

使用squid可以搭建代理服务器，可以用于翻墙访问国外网站。这里讲解通过在校内服务器上搭建squid代理服务，用于校外通过校内的代理服务器下载文献。

1. 在校内服务器上安装squid软件，并启动squid服务。

使用root权限在学校内网中安装有CentOS7系统的服务器上安装squid软件：

yum install squid

再修改squid配置文件，以让所有其它IP地址的用户可以访问校内服务器：

perl -p -i -e 's/http_access deny all/http_access allow all/' /etc/squid/squid.conf

启动squid服务：

systemctl start squid.service

2. 将代理服务器的3128端口映射到外网Aliyun服务器上的31280端口

squid代理服务器使用的默认端口是3128。若此时将内网服务器的防火墙开放3128端口，则在个人电脑上代理设置上填写内网服务器的固定IP和3128端口（见本文第4部分），在个人电脑上的浏览器访问网络时，就是通过代理服务器来访问了。

若是国外服务器上做的代理，则防火墙开放3128端口，可以在国内访问国外代理服务器进行翻墙。

firewall-cmd --add-port=3128/tcp --permanent 
systemctl restart firewalld.service

由于学校内网一般对外封闭了所有的端口，在校外依然是无法使用代理服务器的。因此，需要通过一个具有外网固定IP的服务器来进行中转。我这儿就使用了Aliyun上购买的一台具有固定IP的云服务器进行中转。中转时不需要校内的服务器开放防火墙3128端口。

中转的方法是使用ssh进行反向隧道连接，在内网服务器上输入反向隧道命令，将内网服务器的3128端口映射到Aliyun服务器的31280端口。这样当访问Aliyun服务器的31280端口时，就等同于在内网服务器上直接访问3128端口了。

ssh -f -N -R 31280:localhost:3128 chenlianfu@115.29.105.xxx

3. 在校外有固定IP的Aliyun服务器防火墙上开放31280端口的访问

使用root用户在Aliyun服务器上操作，开放防火墙31280端口，以利于个人电脑能访问Aliyun服务器的31280端口。

firewall-cmd --add-port=31280/tcp --permanent 
systemctl restart firewalld.service

4. 个人Windows10电脑上设置代理访问

按windows+x键，再按字母n，打开系统设置界面；点击网络和internet，再点击代理；在下半个页面的手动设置代理部分，打开使用代理服务器开关，然后在地址栏填写Aliyun服务器的IP地址115.29.105.xxx，在端口栏填写31280，最后点击保存。

然后打开浏览器，输入网址，即可通过校内的代理服务器访问网络。这时等同于在学校连接校园网后的网络访问。当用IE浏览器访问网址，若代理服务器没生效，IE浏览器会报告正在使用某个IP某个端口进行代理访问，代理不成功等信息。

成功后，示例访问文献https://doi.org/10.1016/j.molp.2020.05.005，若能查看文献全文和下载PDF全文，则表示成功。

使用PAML软件利用密码子序列进行正选择分析

发表于2019 年 10 月 25 日由chenlianfu

PAML软件中的codeml命令可以使用Maxmum Likelihood方法对密码子或氨基酸的多序列比对结果进行分析和进化参数计算。这两种数据共用一个命令，程序运行时的参数绝大部分是公用的，但仍需要分清两种数据各自的参数意义。PAML软件中进行正选择分析的示例数据位于examples/lysozyme/目录下。

1.输入文件（1/3）：含有拓扑结构的树文件

输入一个树文件input.trees。该文件可以是有根树，也可以是无根树，也可以带有枝长信息。codeml命令用于分析各分枝上的omega值，若是有根树的话，对于分析root节点所在的分枝，其结果可能是不对的。因此，作者推荐使用无根树进行正选择分析。无根树是3分叉的树，最少需要3个物种；而3个物种的无根数有3个分枝（2n – 3），对其进行分析，仅能得到某个物种所对应分枝的omega值，这种情况下没太大意义；而使用4个物种的无根树进行分析，可以得到两个物种共同祖先分枝的omega值。故作者推荐使用至少4~5个物种进行分析。一个树文件内容示例：

7  1
((Hsa_Human, Hla_gibbon) #1,((Cgu/Can_colobus, Pne_langur) #1, Mmu_rhesus), (Ssc_squirrelM, Cja_marmoset));

文件内容分两行：第一行表述树中有5个物种，共计1个树，两个数值之间用任意长度的空分割。此外，Newick格式的树尾部一定要有分号，没有的话程序可能不能正常运行。

2. 输入文件（2/3）：密码子序列比对结果文件 input.txt

7  390
Hsa_Human        AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Hla_gibbon       AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Cgu/Can_colobus  AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Pne_langur       AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Mmu_rhesus       AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Ssc_squirrelM    AAGGTCTTCGAAAGGTGTGAGTTGGCCAGAACT...
Cja_marmoset     AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...

PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。

3. 输入文件（3/3）：codeml命令的配置文件 codeml.ctl

软件难点在于配置文件的理解，常用示例如下：

*输入输出参数：
      seqfile = input.txt    *设置输入的多序列比对文件路径。
     treefile = input.trees  *设置输入的树文件路径。
      outfile = mlc          *设置输出文件路径。
        noisy = 9            *设置输出到屏幕上的信息量等级：0，1，2，3，9。
      verbose = 0            *设置输出到结果文件中的信息量等级：0，精简模式结果（推荐）；1，输出详细信息，包含碱基序列；2，输出更多信息。
        getSE = 0            *设置是否计算并获得各参数的标准误：0，不需要；1，需要。
 RateAncestor = 0            *设置是否计算序列中每个位点的替换率：0，不需要；1，需要。当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现。

*数据使用说明参数：
      runmode = 0      *设置程序获取进化树拓扑结构的方法：0，直接从输入文件中获取进化树拓扑结构（推荐）；1，程序以输入文件中的多分叉树作为起始树，并使用星状分解法搜索最佳树；2，程序直接以星状树作为起始树，并使用星状分解法搜索最佳树；3，使用逐步添加法搜索最佳树；4，使用简约法获取起始树，再进行邻近分枝交换寻找最优树；5，以输入文件中的树作为起始树，再进行邻近分枝交换寻找最优树；-2，对密码子序列进行两两比较并使用ML方法计算DnDs，或对蛋白序列进行两两比较计算ML距离，而不进行其它参数（枝长和omega等）的计算。
* fix_blength = 1      *设置程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行进行计算；1，以输入的枝长信息作为初始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
      seqtype = 1      *设置输入的多序列比对数据的类型：1，密码子数据；2，氨基酸数据；3，输入数据虽然为密码子序列，但先转换为氨基酸序列后再进行分析。
    CodonFreq = 2      *设置密码子频率的计算算法：0，表示除三种终止密码子频率为0，其余密码值频率全部为1/61；1，程序分别计算三个密码子位点的四种碱基频率，再算四种碱基频率在三个位点上的算术平均值，计算61种密码子频率时，使用该均值进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；2，程序分别计算三个密码子位点的四种碱基频率，计算61种密码子频率时，使用相应位点的碱基频率进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；3，直接使用观测到的各密码子的总的频数/所有密码值的总数，得到所有密码子的频率。一般选择第三种方法，设置CodonFre的值为2。第一种方法让所有密码子频率相等，不符合实际轻卡；第二种方法没有考虑到三种位点各碱基频率的差异，得到的密码子频率不准确；第三种方法能比较准确计算出各密码子的频率；第四种方法由输入数据得到真实的各密码子频率，但有时候有些非终止密码子的频率为0，可能不利于后续的计算。
    cleandata = 1      *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。
*       ndata = 1      *设置输入的多序列比对的数据个数。
        clock = 0      *设置进化树中各分支上的变异速率是否一致，服从分子钟理论：0，变异速率不一致，不服从分子钟理论（当输入数据中物种相差较远时，各分支变异速率不一致）；1，所有分枝具有相同的变异速率，全局上服从分子钟理论；2，进化树局部符和分子钟理论，程序认为除了指定的分支具有不同的进化速率，其它分枝都具有相同的进化速率，这要求输入的进化树信息中使用#来指定分支；3，对多基因数据进行联合分析。
        Mgene = 0      *设置是否有多个基因的多序列比对信息输入，以及多各基因之间的参数是否一致：0，输入的多序列比对文件中仅包含一个基因时，或多个基因具有相同的Kappa和Pi参数；1，输入文件包含多个基因，这些基因之间是相互独立的（这些基因之间具有不同的Kappa和Pi值，且其进化树的枝长也不相关）；2，输入文件包含多个基因，这些基因具有相同的Kappa值，不同的Pi值；3，输入文件包含多个基因，这些基因具有相同的Pi值，不同的Kappa值；4，输入文件包含多个基因，这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时，多个基因的进化树枝长虽然长度不一样，但是呈正比关系。这点和参数值等于1是不一样的。
        icode = 0      *设置遗传密码。其值1-10和NCBI的1-11遗传密码规则对应：0，表示通用的遗传密码。
   Small_Diff = .5e-6  *设置一个很小的值，一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果，需要该参数值对结果没什么影响。

*位点替换模型参数：
        model = 1            *若输入数据是密码子序列，该参数用于设置branch models，即进化树各分枝的omega值的分布：0，进化树上所有分枝的omega值一致；1，对每个分枝单独进行omega计算；2，设置多类omega值，根据树文件中对分枝的编号信息来确定类别，具有相同编号的分枝具有相同的omega值，没有编号的分枝具有相同的omega值，程序分别计算各编号和没有编号的omega值。
                             *若输入数据是蛋白序列，或数据是密码子序列且seqtype值是3时，该参数用于设置氨基酸替换模型：0，Poisson；1，氨基酸替换率和氨基酸的观测频率成正比；2，从aaRatefile参数指定的文件路径中读取氨基酸替换率信息，这些信息是根据经验获得，并记录到后缀为.dat的配置文件中。这些经验模型(Empirical Models)文件位于PAML软件安装目录中的dat目录下，例如，Dayhoff(dayhoff.dat)、WAG(wag.dat)、LG(lg.dat)、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等；
   aaRatefile = dat/wag.dat  *当对蛋白数据进行分析，且model = 2时，该参数生效，用于设置氨基酸替换模型。
       aaDist = 0            *设置氨基酸之间的距离。
      NSsites = 0            *输入数据时密码子序列时生效，用于设置site model，即序列各位点的omega值的分布：0，所有位点具有相同的omega值；1，各位点上的omega值小于1或等于1（服从中性进化neutral）；2，各位点上的omega值小于1、等于1或大于1（选择性进化selection）；3，discrete；4，freq；5:gamma；6，2gamma；7，beta；8，beta&w；9，beta&gamma；10，beta&gamma+1；11，beta&normal>1；12，0&2normal>1；13，3normal>0。
                             *可以一次输入多个模型进行计算并比较，其结果输出的rst文件中。
    fix_alpha = 1            *序列中不同的位点具有不同的碱基替换率，服从discrete-gamma分布，该模型通过alpha（shape参数）和ncatG参数控制。该参数设置是否给定一个alpha值：0，使用ML方法对alpha值进行计算；1，使用下一个参数设置一个固定的alpha值。
	                     *对于密码子序列，当NSsites参数值不为0或model不为0时，推荐设置fix_alpha = 1且alpha = 0，即不设置alpha值，认为位点间的变异速率一致，否则程序报错。若设置了alpha值，则程序认为不同密码子位点的变异速率不均匀，且同时所有位点的omega值一致，当然各分枝的omega值也会一致，这时要求NSsites和model参数值都设置为0（这一般不是我们需要的分析，它不能进行正选择分析了）。
        alpha = 0            *设置一个固定的alpha值或初始alpha值（推荐设置为0.5）。该值小于1，表示只有少数热点位置的替换率较高；该值越小，表示位点替换率在各位点上越不均匀；若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。
       Malpha = 0            *当输入的多序列比对结果中有多基因时，设置这些基因间的alpha值是否相等：0，分别对每个基因单独计算alpha值；1，所有基因的alpha值保持一致。
        ncatG = 5            *序列中不同位点的变异速率服从GAMMA分布的，ncatG是其一个参数，一般设置为5，4，8或10，且序列条数越多，该值设置越大。
                             *对于密码子序列，当NSites设置为3时，ncatG设置为3；当NSites设置为4时，ncatG设置为5；当NSites值设置>=5时，ncatG值设置为10。
    fix_kappa = 0            *设置是否给定一个Kappa值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的Kappa值。
        kappa = 2            *设置一个固定的Kappa值，或一个初始的Kappa值。
    fix_omega = 0            *设置是否给定一个omega值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的omega值。 
        omega = .4           *设置一个固定的omega值，或一个初始的omega值。
       method = 0            *设置评估枝长的ML迭代算法：0，使用PAML的老算法同时计算所有枝长，在clock = 0下有效；1，PAML新加入的算法，一次对一个枝长进行计算，该算法仅在clock参数值为1，2或3下工作。

4. 运行codeml进行分析

运行程序的命令：

codeml codeml.ctl

5. 正选择基因分析

对基因进行正选择分析，可以按如下步骤进行：

(1) 收集该基因在多个物种中的CDS和protein序列（可以使用orthoMCL分析结果）。
(2) 对该基因的蛋白序列进行多序列比对，再根据CDS序列转换成Codon序列比对结果。这个步骤需要自己编写一些程序来实现。
(3) 根据Codon序列的比对结果构建系统发育树。推荐使用RAxML软件使用ML算法对所有的单拷贝同源基因进行物种树构建。然后，输入物种树信息、Codon多序列比对信息，使用codeml进行omega计算和选择压模型检验。
(4) 快速过滤非正选择基因：设置runmode = -2运行codeml命令对两两序列比较，使用ML方法计算dNdS，若存在omega值大于1，则认为该基因可能属于正选择基因，进行后续分析。
(5) 正选择基因的初步鉴定方法：设置model = 0、NSsites = 1 2运行codeml命令（site models）分别对正选择模型(M2a)和近中性进化模型(M1a)进行检测，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。从程序结果中（例如：lnL(ntime: 11  np: 16):   -870.867266      +0.000000）可以找到M2a模型的似然值l₁和M1a模型的似然值l₀，再计算2Δl = 2(l₁-l₀)，通过命令“chi2 2 2Δl”根据卡方检验算出p值。此外，也可以设置NSsites = 7 8进行计算，则表示分别对M8(beta & omega)和M7(beta)模型进行检测，和前者类似，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。根据PAML作者所说，M2a和M1a的比较，比M7-M8的比较更严格。即若想得到更多的正选择基因，可以使用M7_VS_M8分析。此外，根据M2和M8模型的BEB（Bayes empirical Bayes）分析结果，还可以得到在多序列比对结果中第一条序列上的正选择位点。根据PAML作者所说，这种site models之间的比较，能用于检测是否在序列上不同的位点具有不同的omega值，而不是用于正选择检测（We suggest that The M0-M3 comparison should be used as a test of variable ω among sites rather than a test of positive selection）。
(6) 指定分化枝上的正选择基因鉴定方法：经过上一步初步鉴定后，设置model = 2、NSsites = 2、fix_omega = 0、omega = 2.0运行codeml命令（branch-site model A）；再设置model = 2、NSsites = 2、fix_omega = 1、omega = 1运行codeml命令（modified branch-site model A / null model）。进行这两种模型分析时，要求输入的树文件中对目标分化枝进行标注。对这两种模型进行LRT分析，计算2Δl = 2(l₁-l₀)，注意是前者的似然值减后者（null model）的似然值；再使用自由度为1的卡方检验，通过命令“chi2 1 2Δl”计算出的值再除以2，即得到p值。
(7) 目标分化枝上的快速进化基因（Rapidly evolving gene）鉴定方法：设置model = 0、NSsites = 0运行codeml命令，再设置model = 2、NSsites = 0运行codeml命令（这时要求输入的树文件中对目标分化枝进行标注）。然后比较两种运行模式下的似然值，通过自由度为1的卡方检验，可以鉴定该基因在目标分化枝的omega值和背景差异显著。

使用PAML软件利用核酸序列进行系统发育参数分析

发表于2019 年 10 月 23 日由chenlianfu

PAML软件中的baseml命令可以利用Maxmum Likelihood方法计算系统发育中的相关参数。basml命令能以分支树拓扑结构信息和多序列比对信息作为输入，在参数配置文件中设置一些控制程序运行的参数，来对一些进化过程中的参数进行分析。这些参数主要有：

1. 序列中各碱基频率。
2. 两两物种之间的距离。
3. 系统发育树的枝长信息。
4. 碱基替换矩阵。

PAML软件的baseml示例数据位于软件安装目录下。

1. 输入文件（1/3）：含有拓扑结构的树文件 input.trees

输入一个树文件input.trees。该文件可以是有根树，也可以是无根树，也可以带有枝长信息。例如：

5  1
(((Human:0.1, Chimpanzee:0.2):0.8, Gorilla:0.3):0.7, Orangutan:0.4, Gibbon:0.5);

2. 输入文件（2/3）：核酸序列比对结果 input.txt

5   895

Human       AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTT...
Chimpanzee  AAGCTTCACCGGCGCAATTATCCTCATAATCGCCCACGGACTT...
Gorilla     AAGCTTCACCGGCGCAGTTGTTCTTATAATTGCCCACGGACTT...
Orangutan   AAGCTTCACCGGCGCAACCACCCTCATGATTGCCCATGGACTC...
Gibbon      AAGCTTTACAGGTGCAACCGTCCTCATAATCGCCCACGGACTA...