比较基因组——还是看我的教程吧！

一、运行orthofinder

首先 orthofinder使用的版本为2.5.* 不要使用2.2的，2.2默认比对是blast，速度非常慢，结果文件呈现形式也不让人满意。2.5默认用的diamond 速度非常快

第一步代码：
nohup orthofinder -t 40 -f data/ & # data中存放的是蛋白序列文件，建议测试时候最少使用4个物种也就是4份蛋白文件，因为我遇到了建立进化树过程中的报错，说小于4个物种不可以。

data文件内容：

运行结果示例：

第二步代码：
orthofinder -b WorkingDirectory # orthofinder运行结束后，再去运行这条命令，WorkingDirectory是orthofinder生成的，这条命令的目的是什么我也不知道，希望有大佬解读一下，参考的以下：
> https://www.jianshu.com/p/52c2b99615f6?ivk_sa=1024320u

运行完毕后：

##############################################伟大的分割线#####################################################

二、多序列比对并构建系统发育树：

cd Single_Copy_Orthologue_Sequences # Orthofinder的 Single_Copy_Orthologue_Sequences下存放着单拷贝同源基因的序列，我们可以利用这些序列构建系统发育树

$vi bash.sh

for i in *.fa
do muscle -in $i -out $i.1 # 多序列比对推荐使用muscle 记得测试下你的muscle有没有安装，我的安装了orthofinder后自动有这个软件了
done

$sh bash.sh

提取保守序列：
$vi bash2.sh
for i in *.1
do Gblocks $i -b4=5 -b5=h -t=p -e=.2 # Gblocks这个软件也是安装了orthofinder后自动有这个软件了
done
$sh bash2.sh

#####################################################################
以下不用运行，是解释参数，复制的别人的：
-t= default:p设置序列的类型，可选的值是 p,d,c 分别代表 protein， DNA， Codons 。
-b1= default:（ 序列条数的 50% + 1 ）,设定保守性位点必须有 >= 该值的序列数。该参数后接一个 integer 数，默认下比序列条数的 50% 大 1.
-b2= default: 序列条数的 85%,确定保守位点的侧翼位点时，其位点必须有 >= 该值的序列数。该值必须要比 -b1 的值要大。
-b3= default: 8,最大连续非保守位点的长度。
-b4= default: 10,保守位点区块的最小长度。该值必须 >=2 。
-b5= default: n, 设置允许含有 Gap 位点。可选的值有 n,h,a 分别代表 None, With Half, All 。 当为 h 时，表示
-e= default: .2, 设置输出结果的后缀。
作者：周小钊
链接：https://www.jianshu.com/p/52c2b99615f6


序列合并：
$vi bash3.sh
for i in *.2
do seqkit sort $i >$i.3
seqkit seq $i.3 -w 0 > $i.3.4 # seqkit这个软件也是安装了orthofinder后自动有这个软件了，没有的话自行安装，conda就可以
done
$sh bash3.sh

$mkdir new
$mv *.4 new/
$cd new
$paste -d " " *.4 > all.fa  # 这条命令可能会遇到以下报错：paste: OG0008914.fa.1.2.3.4: Too many open files。见图。
# 因此我将其改为R代码

# R代码：
# 获取当前目录下以".4"结尾的文件列表
file_list <- list.files(pattern = "\\.4$")

# 创建一个空的矩阵，用于存储合并后的数据
merged_data <- matrix(nrow = 8, ncol = length(file_list)) #由于我使用了4个物种，所以每个OG开头的*.4结尾的文件里有8行,这里行数要根据自己需要更改

# 逐个读取文件的内容并存储到矩阵中
for (i in 1:length(file_list)) {
  file_path <- file_list[i]
  file_contents <- readLines(file_path)
  
  # 将文件内容按行存储到矩阵的对应列中
  merged_data[, i] <- file_contents
}

# 将矩阵转换为数据框，并设置列名
merged_data <- as.data.frame(merged_data)
colnames(merged_data) <- file_list

# 将数据框写入到文件"all.fa"中，使用空格作为分隔符
write.table(merged_data, file = "all.fa", quote = FALSE, row.names = FALSE, sep = " ")

生成的文件长这样：

可以看到第一行是一堆文件名字，删除就可以。
第二、四、六、八行改成每个物种的名字，我的是之间按照物种蛋白质文件的顺序排列的，别人的应该也一样，如果不确定，使用以下命令检查：

cat ASM1339834v1.faa | grep '随便检索一个>名字' # 如果不确定就多检索几个。

以下是我操作截图：

修改后的文件长这样：

我的蛋白文件长这样：
OrthoFinder/ 文件夹是运行OrthoFinder后生成的。

接下来运行以下命令去除空格：

sed -i "s/ //g" all.fa # 去除all.fa中的空格

接下来使用Python将all.fa转为phy格式

import re
with open('all.fa', 'r') as fin:
    sequences = [(m.group(1), ''.join(m.group(2).split()))
    for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]
with open('all.phy', 'w') as fout:
    fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))
    for item in sequences:
        fout.write('%-20s %s\n' % item)
# 这里我是将以上复制到了一个phy.py的文件中，然后使用python phy.py运行,看以下截图 python版本为3.12.2

接下来使用iqtree2构建系统发育树：

iqtree -s all.phy -m MFP -B 1000 --bnni -T AUTO #参考链接：https://cloud.tencent.com/developer/article/1991339?areaSource=&traceId=        
# 感觉这个运行了好久~

生成了以下截图中的文件，将all.phy.treefile 这个树文件传到ITOL进行美化。


进化树制作到此结束。

##############################################伟大的分割线#####################################################

三、基因家族扩张与收缩分析

首先利用OrthoFinder的Orthogroups.GeneCount.tsv文件生成符合要求的输入文件:

cp Results_May02/Orthogroups/Orthogroups.GeneCount.tsv CAFE/
cd CAFE/
awk 'OFS="\t" {$NF=""; print}' Orthogroups.GeneCount.tsv > tmp && awk '{print "(null)""\t"$0}' tmp > cafe.input.tsv && sed -i '1s/(null)/Desc/g' cafe.input.tsv && rm tmp

awk 'NR==1 || $3<100 && $4<100 && $5<100 && $6<100 {print $0}' cafe.input.tsv >gene_family_filter.txt # 这里有几个物种就写到几加2，比方说我的4个物种，就写到$6<100

运行以下代码：

cafe5 -i gene_family_filter.txt -t SpeciesTree_rooted_node_labels.txt -o out -c 1 -k 2 -p # 这个树文件是不是用这个我不确定
# cafe用conda安装就行
# 结果文件在out中

结果文件展示：

结果文件解读：
Gamma_asr.tre ## 每个基因家族的树文件
Gamma_branch_probabilities.tab ## 每个分支计算的概率
Gamma_category_likelihoods.txt
Gamma_change.tab ## 每一个基因家族在每个节点的收缩与扩增数目
Gamma_clade_results.txt ##每个节点基因家族的扩增/收缩数目
Gamma_count.txt ## 每一个基因家族在每个节点的数目
Gamma_family_likelihoods.txt
Gamma_family_results.txt ## 基因家族变化的p值和是否显著的结果
Gamma_results.txt ## 模型，最终似然值，最终Lambda值等参数信息。

接下来画图：

使用CAFE_fig 软件画图，软件地址：https://github.com/LKremer/CAFE_fig

命令：

python ../CAFE_fig-master/CAFE_fig.py Gamma_report.cafe -pb 0.05 -pf 0.05 --dump test/ -g pdf --count_all_expansions  
# test/ 意思是生成一个test文件结果放在test下

解决CAFE_fig报错：

运行前要先运行以下：

pip install ete3  # 这个github上是pip3 install 'ete3==3.0.0b35' 建议还是使用我的命令，因为使用官网的我遇到了from ete3 import TreeStyle报错的问题
pip install six
pip install numpy

完事打开Python测试一下TreeStyle是否能导入

python3

from ete3 import TreeStyle
# 如果不报错那么也仅仅是成功了一半

如果继续报错说TreeStyle不能导入的问题，那么使用conda新创建一个名为pyqt的环境进行以下命令：

conda creat --name pyqt
conda install pyqt #如果慢可以使用mamba,自行搜索怎么用

这时候重新运行：

pip install ete3  # 这个github上是pip3 install 'ete3==3.0.0b35' 建议还是使用我的命令，因为使用官网的我遇到了from ete3 import TreeStyle报错的问题
pip install six
pip install numpy

完事打开Python测试一下TreeStyle是否能导入

python3

from ete3 import TreeStyle
# 还是那句话，如果不报错那么也仅仅是成功了一半，这时候运行这个你大概不会报错了

接下来你运行时候可能会遇到以下报错：

这个我也没解决，可能是没显卡问题？请大佬指点。因此我使用mac电脑来下载的CAFE_fig
之后还是刚才的操作：
创建pyqt环境并安装pyqt，进入环境运行Python，测试TreeStyle是否能导入，全部成功。
然后运行了CAFE_fig 没有报错：

python ../CAFE_fig-master/CAFE_fig.py Gamma_report.cafe -pb 0.05 -pf 0.05 --dump test/ -g pdf --count_all_expansions  
# test/ 意思是生成一个test文件结果放在test下

test下生成的文件展示：

生成的图很丑，可能是我物种数太少的原因，和CAFE_fig官网一点都不一样，就不展示了，自己调吧，或者有大佬指点一下给个代码？
################################################
找某个节点显著扩张/收缩的基因家族/基因信息

cat Gamma_change.tab |cut -f1,15|grep "+[1-9]" >sample.expanded #提取Gamma_change.tab第15列代表物种sample的扩张的orthogroupsID;  这里为什么要[1-9]呢有没有大佬解读一下
# 这里我的数据里没有加号，我用的cat Gamma_change.tab |cut -f1,15|grep "[1-9]" >sample.expanded 
cat Gamma_change.tab |cut -f1,15|grep "-" >sample.contracted  #提取Gamma_change.tab第15列代表物种sample的收缩的orthogroupsID
grep "sample<13>\*" Gamma_asr.tre > sample_significant_trees.tre # 根据sample ID和编号提取sample分支的基因家族显著扩张或收缩的基因家族树（Gamma_asr.tre文件中默认以p<0.05为标准判断变化是否显著）
grep -E -o "OG[0-9]+" sample_significant_trees.tre > sample_significant.ogs # 提取sample分支显著变化的OG IDs （默认以p<0.05为标准）
awk '$2 <0.01 {print $1}' Gamma_family_results.txt >p0.01_significant.ogs # 以p<0.01为标准提取所有显著扩张或收缩的orthogroupsID（根据情况调整，常用p<0.05或p<0.01）
grep -f sample_significant.ogs p0.01_significant.ogs > sample_p0.01_significant.ogs # 提取以p<0.01为标准判断显著性的sample分支基因家族显著变化的OG IDs。
grep -f sample_p0.01_significant.ogs sample.expanded |cut -f1 >s ample.expanded.significant #提取显著扩张的sample物种的orthogroupsID
grep -f sample_p0.01_significant.ogs sample.contracted |cut -f1 > sample.contracted.significant #提取显著收缩的sample物种的orthogroupsID
grep -f sample.expanded.significant ./OrthoFinder/Results_Oct14/Orthogroups/Orthogroups.txt |sed "s/ /\n/g"|grep "bv" |sort -k 1.3n |uniq >sample.expanded.significant.genes #提取显著扩张的基因列表，假设基因ID是bv的前缀。
# 问：这里要是有多个物种的基因前缀是'bv'呢 求大佬解读
grep -f sample.contracted.significant ./OrthoFinder/Results_Oct14/Orthogroups/Orthogroups.txt | sed "s/ /\n/g"|grep "bv" |sort -k 1.3n |uniq >sample.contracted.significant.genes #提取显著收缩的基因列表，假设基因ID是bv的前缀。
seqkit grep -f sample.expanded.significant.genes sample.pep.fa >sample.expanded.significant.pep.fa #提取显著扩张的基因序列
seqkit grep -f sample.contracted.significant.genes sample.pep.fa >sample.contracted.significant.pep.fa #提取显著收缩的基因序列
# 这段来源：https://yanzhongsino.github.io/2021/10/29/bioinfo_gene.family_CAFE5/

结果示例：

这一节结束！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！

#########################################伟大的分割线###############################################

四、选择压力分析

来啦来了，接下来选择压力分析！

首先准备你的基因组和注释文件，见图：

我只用了其中的
ASM1339834v1_genomic.fna
ASM1339834v1.gff
ASM1339956v1_genomic.fna
ASM1339956v1.gff
这四个文件
第一步，提取蛋白序列和cds序列：

gffread ASM1339834v1.gff -g ASM1339834v1_genomic.fna -x ASM1339834v1_cds.fa -y ASM1339834v1_pep.fa
gffread ASM1339956v1.gff -g ASM1339956v1_genomic.fna -x ASM1339956v1_cds.fa -y ASM1339956v1_pep.fa

生成以下文件：

第二步：diamond 比对 并筛选直系同源基因对

wgd dmd --eval 1e-10 -o dmd_out ASM1339834v1_cds.fa ASM1339956v1_cds.fa

这里可以输入多个物种的cds文件，但是建议两个两个的做。省的后续造成困扰。生成的文件如下：

把生成的文件复制到上一个文件夹内并改名：

cp dmd_out/*.rbh.tsv ./homo_pairs.txt # 为啥命令这么写呢，其实我也不知道，原理应当就是把生成的文件复制到上一个文件夹内并改名，
# 参考的以下：
# https://mp.weixin.qq.com/s/miGmtqR6Uk5vSf502x_ESw

第三步：合并两个物种的cds和蛋白文件：

cat ASM1339834v1_cds.fa ASM1339956v1_cds.fa > ASM1339834v1_ASM1339956v1_input_cds.fasta
cat ASM1339834v1_pep.fa ASM1339956v1_pep.fa > ASM1339834v1_ASM1339956v1_input_pep.fasta

第四步：蛋白比对转cds：

ParaAT.pl -h homo_pairs.txt -a ASM1339834v1_ASM1339956v1_input_pep.fasta -n ASM1339834v1_ASM1339956v1_input_cds.fasta \
-p proc.txt \
-o align_out -m muscle -f axt

生成以下文件夹：align_out
ParaAT.pl下载地址：
https://ngdc.cncb.ac.cn/tools/paraat
我用的2.0那个版本
这一步记得看一下homo_pairs.txt，ASM1339834v1_ASM1339956v1_input_pep.fasta和ASM1339834v1_ASM1339956v1_input_cds.fasta这三个文件内>后面的字符，我的有|这个竖线，它就把竖线作为分隔符了。
然后我用以下命令把竖线删除了：

sed -i 's/|//g' homo_pairs.txt
sed -i 's/|//g' ASM1339834v1_ASM1339956v1_input_pep.fasta
sed -i 's/|//g' ASM1339834v1_ASM1339956v1_input_cds.fasta

第五步：合并生成的文件：

cat align_out/*.axt > merge_align.axt

第六步：计算 kaks

/home/zhang/soft/KaKs_Calculator/KaKs_Calculator3.0/src/KaKs -m YN -i ./merge_align.axt -o result.txt

生成的文件如下：

这一节要感谢 Bioinfor 生信云 给的教程以及Bioinfor 生信云群里朋友给的帮助，感谢！
公众号： Bioinfor 生信云
微信号｜xhh99901
这一节到此结束，我再去肝物种分歧时间！
有缘定见！