漳州网站建设公司推荐,建设很多网站能赚到钱,wordpress手机单页面模板,wordpress英文站群一、运行orthofinder
首先 orthofinder使用的版本为2.5.* 不要使用2.2的#xff0c;2.2默认比对是blast#xff0c;速度非常慢#xff0c;结果文件呈现形式也不让人满意。2.5默认用的diamond 速度非常快 第一步代码#xff1a;
nohup orthofinder -t 40 -f data/ # …一、运行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_sa1024320u运行完毕后
##############################################伟大的分割线#####################################################
二、多序列比对并构建系统发育树
cd Single_Copy_Orthologue_Sequences # Orthofinder的 Single_Copy_Orthologue_Sequences下存放着单拷贝同源基因的序列我们可以利用这些序列构建系统发育树$vi bash.shfor 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 -b45 -b5h -tp -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?areaSourcetraceId
# 感觉这个运行了好久~生成了以下截图中的文件将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 tmpawk NR1 || $3100 $4100 $5100 $6100 {print $0} cafe.input.tsv gene_family_filter.txt # 这里有几个物种就写到几加2比方说我的4个物种就写到$6100
运行以下代码
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 ete33.0.0b35 建议还是使用我的命令因为使用官网的我遇到了from ete3 import TreeStyle报错的问题
pip install six
pip install numpy完事打开Python测试一下TreeStyle是否能导入
python3from ete3 import TreeStyle
# 如果不报错那么也仅仅是成功了一半如果继续报错说TreeStyle不能导入的问题那么使用conda新创建一个名为pyqt的环境进行以下命令
conda creat --name pyqt
conda install pyqt #如果慢可以使用mamba,自行搜索怎么用这时候重新运行
pip install ete3 # 这个github上是pip3 install ete33.0.0b35 建议还是使用我的命令因为使用官网的我遇到了from ete3 import TreeStyle报错的问题
pip install six
pip install numpy完事打开Python测试一下TreeStyle是否能导入
python3from 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 sample13\* Gamma_asr.tre sample_significant_trees.tre # 根据sample ID和编号提取sample分支的基因家族显著扩张或收缩的基因家族树Gamma_asr.tre文件中默认以p0.05为标准判断变化是否显著
grep -E -o OG[0-9] sample_significant_trees.tre sample_significant.ogs # 提取sample分支显著变化的OG IDs 默认以p0.05为标准
awk $2 0.01 {print $1} Gamma_family_results.txt p0.01_significant.ogs # 以p0.01为标准提取所有显著扩张或收缩的orthogroupsID根据情况调整常用p0.05或p0.01
grep -f sample_significant.ogs p0.01_significant.ogs sample_p0.01_significant.ogs # 提取以p0.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/结果示例 这一节结束
##############################################伟大的分割线##################################################### 明天再做物种分歧时间那个有点看不懂怎么填那个基准时间
