（11）Orthofinder同源比对&建树

琉璃月月呀

已于 2023-09-22 20:51:14 修改

阅读量1.9k

点赞数 7

文章标签： linux

于 2023-06-19 17:40:14 首次发布

本文链接：https://blog.csdn.net/weixin_44562189/article/details/131126135

版权

（1）更改蛋白文件的名字，以物种名来命名序列

一、原理相关

OrthoFinder 是一个快速、准确和全面的比较基因组学平台。
它找到正交群(orthogroups)和直系同源(orthologs)，推断所有正交群的有根基因树，并识别这些基因树中的所有基因复制事件。
它还为被分析的物种推断出一个有根的物种树，并将基因复制事件从基因树映射到物种树的分支。
OrthoFinder 还为比较基因组分析提供全面的统计数据。
OrthoFinder 使用简单，运行它所需的只是一组 FASTA 格式的蛋白质序列文件（每个物种一个）。

#总的来说，它将要分析的物种的蛋白质组作为输入，并从这些蛋白质组中：
#推断目标物种的正交群
#推断出一组完整的有根基因树
#推断有根物种树
#使用基因树推断基因之间的所有直系同源关系
#推断基因复制事件并将它们交叉引用到基因和物种树上的相应节点
#为目标物种提供比较基因组学统计数据
#除了大规模分析外，它还可以用于在实验研究之前仔细检查各个直系同源关系。

二、help文档解读

#开始分析的选项:
-f <dir>：从 FASTA 文件目录开始分析
-b <dir>：从 OrthoFinder 目录中的 BLAST 结果开始分析
-b <dir1> -f <dir2>：从 OrthoFinder dir1 中的 BLAST 结果开始分析并添加 FASTA 文件from dir2
-fg <dir>：从 orthogroups OrthoFinder 目录开始分析
-ft <dir>：从 OrthoFinder 目录中的基因树开始分析

#停止分析的选项：
-op：在为全对全序列搜索准备输入文件后停止（例如 BLAST/DIAMOND）
-og：在推断正交群后停止
-os：在为正交群写入序列文件后停止（需要 '-M msa'）
-oa：停止推断多张序列比对后orthogroups（需要“-M MSA”）
-ot：停止推断基因树orthogroups

#控制工作流程的选项
-M <opt>：使用 MSA 或 DendroBLAST 基因树推断，opt=msa,dendroblast [默认=dendroblast]

#控制所有程序的选项
-S <opt>：序列搜索程序 opt=blast,diamond,mmseqs,... 用户可扩展 [默认 = 菱形]
-A <opt>：MSA 程序 opt=mafft,muscle,... 用户可扩展（需要'-M msa') [默认 = mafft]
-T <opt>：树推理程序 opt=fasttree,raxml,iqtree,... 用户可扩展（需要 '-M msa'）[默认 = fasttree]

#更多选项
-d：输入是 DNA 序列 -t <int>：用于序列搜索、MSA 和树推理的线程数 [默认为机器上的内核数]
-a <int>：用于内部、RAM 密集型任务的并行分析线程数[默认值 = 1]
-s <file>：用户指定的有根物种树
-I <int>：MCL 膨胀参数 [默认值 = 1.5]
-x <file>：以 OrthoXML 格式输出结果的信息
-p <dir>：将临时 pickle 文件写入 <dir>
-1：仅执行单向序列搜索
-X：不要将物种名称添加到输出文件中的序列 ID
-y：将 HOG 根部下方的旁系同源进化枝拆分为单独的 HOG
-z：不修剪 MSA（列数>=90% 间隙，最小对齐长度 500）
-n <txt>：添加到结果目录的名称
-o <txt>：非默认结果目录
-h：打印这个帮助文本

三、OrthoFinder 工作流程

1.使用 DIAMOND 软件对输入序列进行 all-vs-all 序列比对。
2.使用 MCL（Markov Cluster Algorithm）算法，根据比对结果进行聚类，得到直系同源组（orthogroup）。
其中每个 orthogroup 的蛋白及序列信息存放在 Orthogroup_Sequences 文件夹中，
单拷贝 orthogroup 的蛋白及序列信息存放在 Single_Copy_Orthologue_Sequences 文件夹中，
orthogroup 的统计信息存放在 Comparative_Genomics_Statistics、Orthogroups 文件夹中。
#如 OrthoFinder 自带案例（ExampleData）中总共包含 2733 个基因，MCL 将 2202 个基因划分为 604 个 orthogroups（gene_num > 2），剩余 531 个基因为离散点（每个基因独立成组）。
3.使用 FastMe 软件对每个 orthogroup（gene_num >= species_num）构建无根基因树（gene tree）。
#如自带案例中总共生成 324 个基因树文件。
4.使用 STAG（Species Tree Inference from All Genes）软件根据 orthogroups（包含所有物种）推断无根物种树（species tree）
#如自带案例推断出的 604 个 orthogroups 中只有 316 个 orthogroups 中的同源基因在所有物种中均有分布，那就对这个316个orthogroups做推断
通过参数 -M dendroblast 或 -M msa，OrthoFinder 可以调用 STAG 中两种构建物种树的方法：DendroBLAST（默认）和 CMSA（Concatenated Multiple Sequence Alignment，联合多序列比对）。
5.使用 STRIDE（Species Tree Root Inference from Gene Duplication Events）通过基因复制事件的不可逆性为无根物种树、无根基因树赋根，
得到有根物种树、有根基因树、基因间的直系同源关系、基因复制事件。
结果存放在文件夹 Species_Tree、Gene_Tree、Orthologues、Gene_Duplication_Events、Comparative_Genomics_Statistics 中。

四、运行代码

（1）更改蛋白文件的名字，以物种名来命名序列

#先要去把用funannotate预测的蛋白文件的名字给改了，不然都长一个样，谁知道谁是谁

$ vim orthofinder_changename.sh #编辑一个更改名字的脚本
#脚本内容如下
#!/bin/bash
# 设置物种名变量
species_name="Umbilicaria_deusta" #这个应该是换哪个名就填哪个哈哈哈，也挺方便不是
# 读取蛋白文件中每一行文本
while read -r line; do
  # 如果当前行包含 ">" 符号，说明是序列头部
  if [[ "$line" == ">"* ]]; then
    # 将序列 ID 和描述信息拆分出来，并在其前面加上物种名并连接
    new_header="${species_name}${line#>}"
    #echo ">${new_header}" # 调试输出新的序列头部
  else
    # 否则说明是序列数据，直接输出即可
   #echo "${line}"
  fi
done < "Umbilicaria_deusta.proteins.fa" > "output.fasta"
#使用该脚本，记得提前改一下species_name的名字，done后面输入和输出的名字

（2）运行orthofinder

$ nohup orthofinder -f 10proteins.orthofinder -S diamond -M msa -T raxml -t 10 -a 10 &
#就是把所有要做orthorfinder的蛋白文件.fasta格式，放到一个文件夹下（data_pep），然后在该文件夹目录下操作
#就是要在有10proteins.orthofinder那个大目录，本次是~/Qxy/knowngene/下

五、Orthofinder输出文件

1. Orthogroups 文件夹###直系同源组目录
Orthogroups.tsv、Orthogroups.txt：记录了 MCL 中成功聚类（直系同源组中基因数 >= 2）的每个 Orthogroup 所包含的基因。
Orthogroups.txt类似于Orthogroups.csv，只不过是OrhtoMCL的输出格式
Orthogroups_UnassignedGenes.tsv：记录了 MCL 中未成功聚类（直系同源组中基因数 >= 1）的离散基因。
Orthogroups.GeneCount.tsv：记录了每个 Orthogroup 中基因在物种间的分布情况，可以用于分析同源基因在物种间的收缩和扩张
Orthogroups_SingleCopyOrthologues.txt：记录了单拷贝直系同源组。

2. Orthogroup_Sequences 文件夹###同源组的序列信息
均为 FASTA 格式文件，记录了每个 Orthogroup 所包含的基因 / 蛋白的序列信息

3. Single_Copy_Orthologue_Sequences 文件夹###单拷贝直系同源序列
均为 FASTA 格式文件，记录了每个单拷贝 Orthogroup 所包含的基因 / 蛋白的序列信息

4. Phylogenetic_Hierarchical_Orthogroups 文件夹 HOG
由于复制本在进化之间存在突变速率的异质性，所以在研究同源基因时更希望所研究的同源基因来自相同的复制本。
Hierarchical Orthogroups（HOG）就是为这一目的而设立的概念，
HOG 指由最近共同祖先中某一基因进化而来的一组直系同源基因，进化过程中不涉及基因复制，所以 HOG 中不包含旁系同源。
#看树，应该就是聚到某一枝上的那些
##OrthoFinder 以物种树中的节点（LCA）为标准，寻找有根基因树内由 LCA【最近共同祖先（last common ancestor，LCA）】中基因进化来的 HOG，
对原先 MCL 算法得到 orthogroup 进行细分。输出文件 N0.tsv，N1.txt，N2.tsv，… 分别指以物种树 N0，N1，N2，… 节点为标准推断出的 HOG。

5. MultipleSequenceAlignments 文件夹
此文件夹仅在 -M msa 模式下输出，均为 FASTA 格式文件。
记录了每个 orthogroup 中序列间的多序列比对结果。
记录了程序通过 CMSA 算法过滤后的 orthogroup 中各序列串联后的多序列比对结果，同时比对结果中空位数 > 50% 的列已被删除【不需要Mafft，只需要Gblocks就行】

6. Species_Tree 文件夹###物种树
SpeciesTree_rooted.txt：STAG、STRIDE 算法计算出的有根物种树结构。
SpeciesTree_rooted_node_labels.txt：相比上树在节点处具有标签（N0，N1...Nm）让后续的分析中可以方便的指定物种树节点。
Orthogroups_for_concatenated_alignment.txt：仅在 -M msa 模式下输出，列出了所有串联起来用于推断物种树的 orthogroup ID。

7. Gene_Trees 文件夹###基因树
记录了每个 orthogroup（gene_num >= 4）的有根基因树结构。

8. Gene_Duplication_Events 文件夹###基因重复事件
注意！OrthoFinder 只统计记录支持值（Support） >= 50% 的的复制事件。
支持值是指复制后两个基因副本未被丢失的比例，Support >= 50% 表示复制后至少有一半基因在演化中保留了下来。
Duplications.tsv：记录了程序推测出的所有基因复制事件的信息。
其中 Species Tree Node 表示基因复制事件发生时所对应的物种树节点（即复制是在该物种内发生的）；
Gene tree node 表示基因复制事件发生时所对应的基因树节点与基因复制事件对应的节点；
Support 表示复制后两个基因副本未被丢失的比例；
Type 中 Terminal 表示重复发生在物种树的末端分支上，Non-Terminal 表示重复发生在物种树的内部分支上，被多个物种共享；
Genes 1、Genes 2 为基因列表，其中 Genes 1 表示来自复制后基因的一个副本；Genes 2 表示来自复制后基因的另一个副本。
SpeciesTree_Gene_Duplications_0.5_Support.txt：记录了物种树每个节点、分枝上包含的基因复制事件的总和，格式为节点或物种名 + 数字（基因复制事件数量）。

9. Orthologues 文件夹###直系同源组
以物种为单位，记录了每个物种与其他物种间的直系同源基因。

10. Comparative_Genomics_Statistics 文件夹###比较基因组
Statistics_Overall.tsv：记录了有关 orthogroup 的常规统计信息。
Statistics_PerSpecies.tsv：以物种为单位，记录了有关 orthogroup 的常规统计信息。
Orthogroups_SpeciesOverlaps.tsv：记录了每个物种对之间共享的 orthogroup 数。
Duplications_per_Species_Tree_Node.tsv：记录了物种树中每个节点、物种中发生基因重复事件的数量。
Duplications_per_Orthogroup.tsv：记录了每个 orthogroup 中推断出的基因重复事件数量。
OrthologuesStats _ *：记录了每对物种之间一对一、一对多和多对多关系的直向同源物数量。

11. WorkingDirectory 文件夹
OrthoFinder 运行所需的必须中间文件，如 DIAMOND 比对结果，STAG 输出的无根物种树等。

12.其他
Citation.txt###若是引用OrthoFinder时，记得添加引用
Log.txt###用的命令行和数据的名字等信息
Phylogenetically_Misplaced_Genes
Putative_Xenologs
Resolved_Gene_Trees

六、使用单拷贝直系同源基因构建物种树

#利用上面Orthofinder结果目录中Single_Copy_Orthologue_Sequences下的文件
#该文件夹下有多个序列文件，就是所有的单拷贝直系同源基因，每一个序列文件中是有多条序列的，这些序列就来自于之前做Orthofinder同源比对的那些物种
#下面需要做的是先把这些序列比对、裁剪，再按照物种串联起来，建树

（1）使用Muscle进行序列比对

$ mkdir SingleCopyOrthologues_speciestree #新建一个树的文件夹 
（本次在~/Qxy/knowngenome/10proteins.orthofinder目录下，和OrthoFinder同级）
$ cd SingleCopyOrthologues_speciestree #进入该文件夹进行后续建树操作

#使用Muscle进行比对(速度快，准确率也可以，hh其实是因为MAFFT教程没看懂)
$ muscle -align seq.fasta -output aln.fasta #单独一个文件可以用这个代码进行比对（这个文件里包含多个序列）
$ vim Muscle_align.sh #多个序列文件就要用脚本进行比对了 !!!每次用之前看一下脚本啊！有些路径什么的是需要改的

#脚本Muscle_align.sh中的内容

#!/bin/bash
# 创建一个名为aligments的文件夹
mkdir -p alignments
# 遍寻Single_Copy_Orthologue_Sequences文件夹中的所有.fa文件
for file in ~/Qxy/knowngenome/10proteins.orthofinder/OrthoFinder/Results_May29/Single_Copy_Orthologue_Sequences/*.fa; do
  # 获取没有扩展名的文件名
  name=$(basename "$file" .fa) #意思就行XX.fa 文件，只保留XX

  # 使用muscle比对序列并将比对结果输出到alignments文件夹中
  muscle -align "$file" -output "alignments/${name}_aligned.fasta" #比对后的文件名都叫XX_aligned.fasta
done

#运行完毕将会得到比对后的序列文件 #会在前台流过数据，不要怕

（2）使用Gblocks对序列进行修剪

#安装Gblocks
$ conda create -n Gblocks #使用conda创建一个环境叫Gblocks
$ conda activate Gblocks #激活这个环境
$ conda install -c bioconda gblocks #conda安装gblocks 好像安在/ifs1/User/wuqi/anaconda3/envs/Gblocks这里了

$ conda activate Gblocks #激活这个环境

$ vim Gblocks_xiujian.sh #创建Gblocks脚本 !!!每次用之前看一下脚本啊！有些路径什么的是需要改的

#!/bin/bash

# 进入需要进行序列修剪的目录
cd ~/Qxy/knowngenome/10proteins.orthofinder/SingleCopyOrthologues_speciestree/alignments

# 循环处理每一个序列文件
for file in *.fasta
do
    # 定义输出文件名
    output="../Gblocks/ ${file%.fasta}_Gblocks.fasta" #好像没有按照这个输出
    
    # 运行 Gblocks 进行序列修剪
    Gblocks "$file" -t "$output" -b4=5 -b5=h -t=p -e=.2 
    
    # 输出处理完成的序列文件名
    echo "Processed $file -> ../Gblocks/"
done

#运行完毕得到两个文件，一个是XX.fasta.2，一个是XX.fasta.2.htm 都在当前目录下
#会在前台流过数据，不要怕

(3)使用seqkit串联序列

$ vim seqkit.sh  #写一个串联的脚本  

for i in *.2;do seqkit sort $i >$i.3 #对于所有后缀是.2的文件，做排序，并将结果输出，生成一个后缀为.3的文件
seqkit seq $i.3 -w 0 > $i.3.4 # 从文件$i.3中提取序列，并将序列输出为单行（行末无空格），生成新的文件.3.4中
done  #循环结束

$ mkdir seqkit #建一个新的文件夹
$ mv *.4 seqkit/ #移动所有后缀为.4的文件到该目录下
$ cd seqkit  #进入该目录
$ paste -d " " *.4 > all.fa #将所有扩展名为.4的文件按列粘贴，以空格为分隔符 
#可能会显示文件数目超出限制，查看一下，再修改一下
$ lsof | wc -l #查看当前打开了多少文件
$ ulimit -n 65535 #使用此命令设置一下文件数量
$ sed -i 's/ //g' all.fa #使用sed命令，将输出文件all.fa中所有空格替换为空，删除空格

#也会有数据在前台流过，不要怕 #最终得到的seqkit_all.fa文件就是串联好的序列啦

(4)使用记事本修改一下序列的名称

$ vim all.fa #用vim编辑器进行修改
#因为一行名字很多，所以一行很长很长，不用怕，视图模式下光标定位到要删除的位置，按d和$，然后这一行就全删掉了
#然后就可以把中间.2 /.3 /.4 /.htm等文件都删掉了，只保留最初单个比对的和比对后串联的

(5)使用raxml建树

也可以将上游弄好的.fa文件导出到本地，去CIPRES在线网站上用raxml建树

$ raxmlHPC-PTHREADS-SSE3 -f a -x 123456 -p 123456 -s all.fa -m PROTGAMMALGX -N 1000 -n output -k -T 8

-f a 用于选择RAxML运算的算法，a表示执行快速Bootstrap分析并搜索最佳得分的ML树
-x 12345 指定一个int数作为随机种子，以启用快速Bootstrap算法
-p 12345 指定一个随机数作为parsimony inferences 的种子
-# 1000 或 -N 1000 自展值为1000
-m GTRGAMMA 指定核苷酸替代模型，PORT表示氨基酸替代模型，GTR表示碱基替代模型， GAMMA表示使用GAMMA模型 就DNA填GTRGAMMA，AA填PROTGAMMALGX
-s input.phy 指定输入文件，.phy格式的多序列比对文件。也可以用.fasta文件，软件包中包含一个程序来将fasta格式转换为.phy格式
-n 输出文件
-T 30 指定多线程运行的CPUs
-k 指定引导的树应该打印分支长度。default:OFF

#建树完成后用iTOL美化，AI自己修，或者在R中用ggtree做（教程https://mp.weixin.qq.com/s/Smgl5CAtZI1xzw6_iOyJBQ）

七、其他

（1）将fasta格式转为phy格式，听说后面PAML也要用到

$ cat xx.fa |tr '\n' '\t'|sed 's/>/\n/g' |sed 's/\t/      /'|sed 's/\t//g'| awk 'NF > 0' > xx.phy.tmp

#cat xx.fa |tr '\n' '\t' (将换行符替换为制表符) 
#sed 's/>/\n/g' (将每个序列名前面的>符号替换为换行符)
#sed 's/\t/      /' (将每行第一个的制表符替换为多个空格) 
#sed 's/\t//g' (删除剩余的制表符,使序列连成一条线) 
#wk 'NF > 0' (删除空行)
#xx.phy.tmp (临时保存)

$ awk '{print "  "NR"  "length($2)}'  xx.phy.tmp|tail -n 1 | cat -  xx.phy.tmp >  all.phy

#awk '{print "  "NR"  "length($2)}'  supergene.phy.tmp (计算序列长度)
#tail -n 1 (返回最后一行,包括了行数与序列长度)
#cat -  xx.phy.tmp (-为上一管道数据) >  xx.phy(最终文件)

（2）统计同源基因数目

$ cd ~/Qxy/knowngenome/10proteins.orthofinder/OrthoFinder/Results_May29/ #进入到Orthofinder结果目录下
$ cd Orthogroups #进入该目录
$ less Orthogroups.GeneCount.tsv  #该表格是同源基因的比对结果，下载到本地然后可以统计同源基因数目

（3）根据其他同源基因数目建树（剔除某一个、两个物种）

示例：

#无bin0，bin1和其他为1的单拷贝基因建树
#根据Orthogroups.GeneCount.tsv，用EXCEL进行筛选，得到相关名称，写入nobin0.txt并上传到服务器

$ mkdri nobin0_tree
$ cd nobin0_tree

$ vim Ortho_byname_getfa.sh

#!/bin/bash
while read i; do
  cp "$i.fa" ~/Qxy/knowngenome/10proteins.orthofinder/SingleCopyOrthologues_speciestree/nobin0_tree
done < ~/Qxy/knowngenome/10proteins.orthofinder/SingleCopyOrthologues_speciestree/nobin0_tree/nobin0.txt

#将nobin0.txt中包含的序列名字，从某目录下提取出来

$ sed -i 's/\r$//' .txt #如遇报错，执行该代码，去除 Windows 格式的文本文件中的回车符号（^M 或 \r），使其在linux系统下正常显示

八、普通建树流程

1）准备好相关序列fasta格式

2）上传至服务器上

3）使用Muscle进行序列比对

$ muscle -align all.fasta -output aln_all.fasta
#输入待比对的序列，以及提前确定一下比对输出序列名称

注：如果muscle不行，那咱们改用mafft比对

$ mafft --auto input.fasta > output.fasta 
#之前orthofinder建树时比对用的是muscle，速度快，但是只能比对短的序列（长度<500,序列数<500）
#这次muscle比对出来是空，所以选择mafft来比对

4）Gblocks去除空白序列

$ conda activate Gblocks
$ Gblocks aln_all.fasta -t output.fasta -b4=5 -b5=h -t=p -e=.2
#用Gblocks去除空白序列，输入待去除文件和输出文件名称，但是好像输出文件是输入文件后缀加.2的

5）sed去除空格

$ sed -i 's/ //g' aln_all.fa
#因为Gblocks出来后会有一些空格，所以这一步要再去掉

6)再将序列传输到本地，利用CIPRES在线网站进行建树

琉璃月月呀

关注

7
点赞
踩
7

收藏

觉得还不错? 一键收藏
0
评论
（11）Orthofinder同源比对&建树

如 OrthoFinder 自带案例（ExampleData）中总共包含 2733 个基因，MCL 将 2202 个基因划分为 604 个 orthogroups（gene_num > 2），剩余 531 个基因为离散点（每个基因独立成组）。SpeciesTree_Gene_Duplications_0.5_Support.txt：记录了物种树每个节点、分枝上包含的基因复制事件的总和，格式为节点或物种名 + 数字（基因复制事件数量）。输出文件 N0.tsv，N1.txt，N2.tsv，…
复制链接

扫一扫