WeFun blog
0%

Orthofinder查找共有和独有基因

发表于
本文字数： 6.6k 阅读时长 ≈ 6 分钟

为什么要查找共有和独有基因？
共有基因一般是同源基因，独有基因是指该物种因为环境而被迫演化出的的基因，这意味这在对应的环境下，必定有一些适应该环境的特殊基因，这可以帮助我们弄清楚物种的之间的演化过程，也可以在育种方面提供原材料。
共有基因往往对应的是一些基因家族，在漫长的演化过程中，基因家族中有一些基因减少，或者增多，这其实就对应着基因家族的扩张收缩，这也是研究物种演化的关键数据之一。
查找共有和独有基因其实就是做基因家族的聚类，分析基因家族的扩张收缩，常有的软件有orthofinder和orthomcl，orthofinder的集成度更好，还可以自动完成建树因此此处我们学习orthofinder的使用。

1.安装

下载和官方教程地址：https://github.com/davidemms/OrthoFinder
可以使用conda安装 

1
conda install orthofinder -c bioconda

也可以手动安装 

1
2
wget https://github.com/davidemms/OrthoFinder/releases/latest/download/OrthoFinder.tar.gz
tar xzvf OrthoFinder.tar.gz

2.准备数据

本次所使用的数据是Magnolia_biondii.chr.pep和Msin.chr.pep两物种的蛋白质数据，但这两个数据往往不可以直接获得，就算直接获得的，也往往不是最长转录本，因此我们先要进行一些处理
首先需要根据数据提取最长转录本
为什么要提取最长转录本？
由于可变剪切的存在，通常一个基因可以转录为多个转录本。但是如果将多个转录本同时进行分析，那么分析会因此受到影响。所以，目前的解决办法是，选取一个最具代表性的转录本（最长转录本）来进行分析。
最长转录本其实可以从序列文件中获取，也可以从gff文件中获取，获取的方式也有很多，可以直接使用gettranstool，用

1
GetLongestTransFromGencode --file example.fa --outfile longest_trans_gencode.fa

官方的教程：https://github.com/junjunlab/GetTransTool，有需要自己可以查阅
在此介绍对于小白最友好的方法，即使用之前提过多次的软件，TBtools： 

sequence toolkit->fasta tools->fasta get representative

note部分不需要管，ID grouping pattern处为正则表达式，这是序列名称的编号规则，一般不需要更改，然后是输入蛋白质序列数据以及输出文件名的设置
建议直接使用TBtools最简单，只需要弄清楚ID号的正则表达式即可

##蛋白质序列的序列名应该尽量简单，否则可能会出现报错

ValueError: too many values to unpack

这种情况一般在10个以上物种进行聚类的时候出现，出现这种情况时，可以直接使用python脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#!/usr/bin/env python3
import os
import re
from Bio import SeqIO
import sys

def simplify_fasta_headers(input_file, output_file):
    """简化单个FASTA文件的序列头，保留第一个有效标识符"""
    with open(output_file, "w") as out_handle:
        for record in SeqIO.parse(input_file, "fasta"):
            new_id = record.id.split()[0]  # 取第一个空格前的部分
            record.id = new_id
            record.description = new_id
            SeqIO.write(record, out_handle, "fasta")

def process_directory(input_dir, output_dir):
    """批量处理目录下所有.fasta文件"""
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    for filename in os.listdir(input_dir):
        if filename.endswith(".fasta") or filename.endswith(".fa"):
            input_path = os.path.join(input_dir, filename)
            output_path = os.path.join(output_dir, filename)
            simplify_fasta_headers(input_path, output_path)
            print(f"处理完成: {filename} -> {output_path}")

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print(f"用法: {sys.argv[0]} <输入文件夹> <输出文件夹>")
        print(f"示例: {sys.argv[0]} ./input_fasta ./cleaned_fasta")
        sys.exit(1)

    input_dir = sys.argv[1]
    output_dir = sys.argv[2]
    process_directory(input_dir, output_dir)
    print(f"\n所有文件处理完成！结果保存在: {output_dir}")

该脚本是终端脚本，使用方式为

1
python simplify_headers.py ./原始文件夹 ./输出文件夹

删除所有非标准氨基酸字符

1
for i in data4/*;do sed -i '/^>/! s/[^A-Za-z]//g' "$i"; done

3.家族聚类及结果解读

1
2
3
4
5
6
orthofinder -f data -S diamond -M msa -T fasttree -t 20 
-f： 数据目录，先将所有蛋白放进文件夹data中
-S： 比对方法，blast, mmseqs, blast_gz, diamond
-M： 基因树推断方法，dendroblast，msa
-T：建树软件，iqtree, raxml-ng, fasttree, raxml
-t： 线程数，根据服务器配置选择，默认为56

聚类结果如下：

Results_XX #结果文件夹名
Citation.txt
Comparative_Genomics_Statistics #统计结果
Gene_Duplication_Events
Gene_Trees # 每个家族基因树
Log.txt
MultipleSequenceAlignments # 多序列比对结果 包含单拷贝基因supergene
Orthogroups # 主要结果目录
Orthogroup_Sequences # 每个家族序列文件
Orthologues
Phylogenetically_Misplaced_Genes
Phylogenetic_Hierarchical_Orthogroups
Putative_Xenologs
Resolved_Gene_Trees
Single_Copy_Orthologue_Sequences # 单拷贝基因家族序列
Species_Tree # 物种树构建结果
WorkingDirectory

Results_XX/Orthogroups#主要结果
Orthogroups.GeneCount.tsv # 基因数目
Orthogroups_SingleCopyOrthologues.txt # 单拷贝基因家族列表
Orthogroups.tsv # 聚类结果
Orthogroups.txt # 聚类结果
Orthogroups_UnassignedGenes.tsv # 未聚类基因

4.处理聚类结果并绘制韦恩图

1
2
3
4
awk '{ if(NR==1){ for(i=2;i<NF;i++ ){printf $i"\t"} }else{ \
for(i=2;i<NF;i++){ if($i>0){printf $1 } ; printf "\t" }}; \
printf "\n"}' Orthogroups.GeneCount.tsv | \
sed 's/\t$//' > Orthogroups.GeneCount.venn

脚本解释： 

NR==1: 当读取到文件的第一行时（通常是列标题行），执行以下操作：
    for(i=2;i<NF;i++): 从第2列开始（跳过第1列），循环处理每一列直到倒数第2列。
    printf $i"\t": 输出每列的内容，后面跟一个制表符（\t），这些输出结果是列标题（从第2列开始），将用于标记维恩图中不同组的名称。
    else{ for(i=2;i<NF;i++){ if($i>0){printf $1 }; printf "\t" }}
非第一行时
    for(i=2;i<NF;i++): 对每一行，从第2列开始直到倒数第2列进行循环。
    if($i>0){printf $1 }: 如果当前列的值大于0，则输出该行的第1列内容（通常是一个基因组的名称或ID）。
    printf "\t": 在输出内容之后，加入一个制表符。

每处理完一行后，脚本在行尾输出一个换行符（\n），以确保每行的数据换行。

1
2
3
4
sed 's/\t$//'
    sed命令用于对awk生成的输出进行后处理。s/\t$//是一个替换操作：
        \t$表示行末的制表符（\t）。
        //表示将匹配到的行末制表符删除。

或者直接使用dos2unix进行转换

处理好的.venn文件可以直接拿来做韦恩图，可以使用在线网站也可以用R，通常情况下我会使用jvenn网站绘制完，根据数据使用AI自己重新绘制一个
jvenn网址：https://jvenn.toulouse.inrae.fr/app/example.html
下面介绍使用R绘制
下面是脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# 加载必要的库
library(ggVennDiagram)
library(ggplot2)

# 读取文件数据
file_path <- "Orthogroups.GeneCount.venn"  # 请根据实际文件路径进行替换
data <- read.table(file_path, header = TRUE, sep = "\t", stringsAsFactors = FALSE)

# 提取非空值的集合
set_JNYL <- data$JNYL[data$JNYL != ""]
set_WCYL <- data$WCYL[data$WCYL != ""]

# 创建集合列表
sets <- list(
  JNYL = set_JNYL,
  WCYL = set_WCYL
)

# 使用ggVennDiagram绘制维恩图
ggVennDiagram(sets,
              category.names = c("JNYL", "WCYL"), # 集合名字
              set_color = c("blue", "red"),       # 集合颜色
              set_size = 6,                       # 集合名字大小
              label = "both",                     # 标签显示"count"和"percent"
              label_geom = "label",               # 标签样式
              label_alpha = 0,                    # 标签背景透明
              label_color = "firebrick",          # 标签颜色
              label_percent_digit = 2,            # 百分比保留2位小数
              edge_lty = "dashed",                # 边框线型
              edge_size = 1.2                     # 边框粗细
) +
  scale_fill_gradient(low = "grey90", high = "grey60") +  # 填充色
  scale_color_manual(values = c("grey10", "grey10"))      # 边框颜色

结果可以使用ggplot2进行优化，这也是为什么要调用ggplot2的原因个人感觉其实到这就可以了
#先将文件转为数据框
venn <- Venn(sets)
df <- process_data(venn)
df
#ggplot2绘图
ggplot()+
  geom_sf(data = venn_region(df), aes(fill=count))+
  geom_sf(data = venn_setedge(df),size=2,lty="dashed",color="grey")+
  geom_sf(data = venn_setlabel(df),aes(label=name))+
  geom_sf_label(data = venn_region(df),aes(label=id),fontface="bold")+
  scale_fill_distiller(palette = 5)+
  theme_void()

其实到这我对图还是不满意，我后面又用AI换了一下颜色

提取独有共有基因的OGid
使用jvenn，点击venn图中想要继续分析的区域的数字，即可获得该区域内的所有OGid

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#将OGid与基因名相匹配
cp ../../Results_May02/Orthogroups/Orthogroups.tsv ./	#将之前结果文件复制到当前目录下
grep -f Ms.txt Orthogroups.tsv | cut -f3 | sed "s/ /\n/g" | sed "s/\t/\n/g" | sed "s/,//g" | sort | uniq > Ms.specific	
#Ms.txt是得到的OGid，必须是unix格式，即使用cat -A打开后结尾是^I^M$，如果不是使用dos2unix转换一下
	-f表示从文件中读取关键词
	| cut -f3 cut 用于按列提取内容，-f3 表示提取每行的第 3 列，因为Ms在Orthogroups.tsv的第三列
"s/ /\n/g"：
    s：表示替换命令。
    / /：表示匹配空格字符。
    \n：表示将空格替换为换行符。
    g：表示全局替换，即替换行中所有的空格
| sed "s/\t/\n/g"：将制表符（\t）替换为换行符
| sed "s/,//g"：将逗号（,）移除
| sort：用于对文本内容进行排序
| uniq： 用于去除相邻的重复行

提取得到基因列表，可以按之前的教程继续后面的GO、KEGG分析