查看: 9058|回复: 15

[动植物重测序] Phylogenomics-DensiTree绘制详细教程

  [复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
30
注册时间
2017.4.3
在线时间
117 小时

发表于 2021.3.24 21:24:35 | 显示全部楼层 |阅读模式
本帖最后由 静默人声412 于 2021.3.24 21:41 编辑

   所谓DensiTree,其实就是将多颗进化树的拓扑结构进行的叠加,以可视化进化树间的拓扑冲突(或基因树异质性)。绘制DensiTree绘制可以使用DensiTree软件(现在已经整合到BEAST2安装包中),也可以使用R包phangorn进行。下面记录一下DensiTree的绘制过程。

数据:某属物种的全基因组数据,基于氨基酸序列,先使用Orthofinder进行了同源基因鉴定,然后提取单拷贝直系同源基因,使用MAFFT进行多序列比对,Trimal提取保守区,过滤掉低于55aa序列,然后使用IQtree进行单基因树构建,最后使用DensiTree软件进行DensiTree绘制。
一、数据准备-物种氨基酸序列提取
# 1. 探索物种基因组注释文件,这里使用的是鞍带石斑鱼grep -v '^#' fEpiLan.gff|awk '{print $3}' |sort|uniq # 序列类型
awk -v FS="\t|;" '{print $10}' fEpiLan.gff |awk '/Dbxref=GeneID*/{print $0}'|sort|uniq|wc -l # 基因ID数量是否与基因数量对应,假基因也存在基因ID,ID数目应该大于gene数目
grep -v '^#' fEpiLan.gff | awk '{if ($3=="gene")print $0}' |wc -l
#基因数目awk -v FS=";" '{print $6}' fEpiLan.gff |grep "gene=*" |sort |uniq|wc -l #基因name的名称数量,## 基因ID是唯一的,基因name有重叠
# 2. 使用gffread提取CDS和对应的蛋白序列
gffread fEpiLan.gff  -F \
-g GCF_005281545.1_ASM528154v1_genomic.fna \
-x fEpiLan_cds.fa -y fEpiLan_pep.fa # 使用的是鞍带石斑鱼基因组数据
# 3. 对提取的序列名称进行重命名,只保留转录本与基因名
grep ">" fEpiLan_pep.fa |awk -v FS=" |," '{print $1}' |sort|uniq|wc -l
grep ">" fEpiLan_pep.fa |awk -v FS=" |," '{print $2}' |sort|uniq|wc -l
grep ">" fEpiLan_pep.fa |awk -v FS=" |," '{print $3}' |sort|uniq|wc -l # 图1,序列名称第二列是基因名称与基因ID对应数量一致,第一列为转录本名称
grep ">" fEpiLan_pep.fa |awk -v FS=" |," '{print $2}' |sort|uniq -d|head  # 图2#
* uniq -d;只输出重复行
#work directory(WD): ~/Sparus_latus/phylogenomics/Aligniment_Genomemkdir homology
python change_ensemble_pep_name.py Epinephelus_lanceolatus/fEpiLan.NCBI.pep.fa homology/fEpiLan.NCBI.pep.fa # 序列名称只保留转录本和基因名称


图1


图2


图3
# 4. 提取最长转录本## 序列名称本身带有|,需要先修改为_
#sed -i "s/|/_/g" Perca_fluviatilis/fPerFlu.NCBI.pep.fa
python extract_longest_cds.py homology/pep/fEpiLan.NCBI.pep.fa  homology/long_pep/fEpiLan.NCBI.long.pep.fa
sed -i 's/>/>fEpiLan\|/g' homology/pep/fEpiLan.NCBI.pep.fa
sed -i 's/>/>fEpiLan\|/g' homology/long_pep/fEpiLan.NCBI.long.pep.fa # 图5
# * 选择序列名称分割符时,要确保基因和转录本名称中没有相同符号


图4


图5
# 5. 最长转录本文件探索,gffread将CDS翻译为氨基酸序列时,出现了.号,可能是终止密码子,但是后面又发现U(终止密码子),所以我当时也有点懵,但是把.号去掉才能进行后续分析
##接下来就挨个对序列进行处理,序列末尾的.删除,中间的.之后再考虑(或者直接将序列作为假基因整体删除),对所有序列进行处理,过滤掉size<50的序列
## 查看文件中有多少.号,图6
tr -cd '/.' <fEpiLan.NCBI.long.pep.fa|wc -c
grep -o  "\." fEpiLan.NCBI.long.pep.fa|wc -l
awk -v RS='\.' 'END {print --NR}' fEpiLan.NCBI.long.pep.fa
## 序列名称中是否存在.号,图7
grep ">" fEpiLan.NCBI.long.pep.fa|grep -o "\." |wc -l## 序列末尾存在多少.号
grep "\.$" fEpiLan.NCBI.long.pep.fa|wc -l
## 先将序列末尾的.号删除
sed 's/\.$//g' fEpiLan.NCBI.long.pep.fa > ../con_homology/fEpiLan.NCBI.long.pep_2.fa
## 序列中间存在.的基因作为假基因被提取出来,只保留不带.的序列,并去掉size<50的序列
python ~/phylogenomics/Scripts/extract_Pseudogenes.py fEpiLan.NCBI.long.pep_2.fa fEpiLan.NCBI.long.pep_gene.fa fEpiLan.NCBI.long.pep_pseudogenes.fa fEpiLan.NCBI.long.pep_small.fa
grep -v ">" fEpiLan.NCBI.long.pep_gene.fa|grep "\." # 验证序列中是否还存在.号,图8
# * 可以将上述命令打包成一个bash文件,运,行图9。


图6


图7


图8


图9
二、直系同源蛋白编码基因分析
# 1. 同源基因鉴定,orthofinder,OMA等都可以进行同源基因鉴定##  将进行同源基因鉴定的物种链接到同一文件夹中
mkdir DB && cd DB
ln -s ~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/con_homology/*fa .
## orthfinder 进行同源基因分析,氨基酸序列中出现.会报错
#conda activate python27
#使用conda 安装orthofinder
#conda install orthofinder
#nohup bash orthofinder.sh > orthofinder.log 2>&1 & 后台运行
wget -t 0 -c https://github.com/davidemms/OrtthoFinder.tar.gz
tar -zxvf OrthoFinder.tar.gz # 手动安装orthofinder -h # 运行命令,显示帮助信息
#orthofinder -f con_homology/ -a 40 -M msa -T iqtree # 参数含义查询官方帮助文档
nohup orthofinder -f DB/ -a 40 -S blast -M msa -T raxml-ng -o ./orthofinder_results >orthofinder_14.log 2>&1 & # 使用blast比对,会有些慢,可以换成diamond
* 注意目录中不能有隐藏文件否则也会参与同源基因的鉴定



图10 orthofinder同源基因鉴定输出结果目录。
三、构建单基因树
# 根据氨基酸同源序列对应的提取同源核酸序列,
##这一步其实是为了对核酸序列按照codon进行比对,以进行后续的趋同进化分析,
##今天只绘制DensiTree,所以后面部分就等有时间了再把流程单独发出来
for i in `ls Orthogroup_Sequences/*fa`;do
grep ">" $i | while read id;do
echo $i $id;done
done # 封装为output_cda.sh
脚本nohup bash output_cda.sh >output_cda.tab & # 同源序列文件名与序列名称对应表
#WD : ~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/con_homology
cat *cds* > ../14_cds.fa # 所有物种的最长转录本合并为一个文件
#WD:~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/OMA/orthofinder_results/Results_Jan01
cat output_cda.tab|while read id;do
arr=($id);
file=`basename ${arr[0]}`;
name=${file%.*}
grep -A1 ${arr[1]} ~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/14_cds.fa >>Orthogroup_cds_Sequences/$name_cds.fa;
done
echo "split cds done" # 封装为split_cds.sh脚本nohup bash split_cds.sh & # 将cds序列按照同源氨基酸序列分类


# Step 1 同源氨基酸多序列比对,使用orthofiner-Orthogroup_Sequences结果文件中序列
#WD:~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/OMA/orthofinder_results/Results_Jan01
mafft --auto Orthogroup_Sequences/OG0003073.fa  > mafft_out/OG0003073.mafft.fa # 可以写脚本,批量进行比对
#mafft --anysymbol --auto Orthogroup_Sequences/OG0003073.fa  > mafft_out/OG0003073.mafft.fa # 序列中有U存在,mafft会报错,可以加上--anysymbol
#nohup bash /B313/public_db/Sparus_latus/phylogenomics/Aligniment_Genome/homology/scripts/mafft.sh Orthogroup_Sequences  >mafft.log 2>&1 &
# Step 2 trimal提取同源基因保守序列,将U替换为X,trimal才能进行修剪,图11
for i in `cat stop_codon_mafft.txt`;do
sed -i '/^[^>]/ s/U/X/' mafft_out/$i;
done
# stop_codon_mafft.txt,包含U的mafft比对结果文件名
nohup bash scripts/trimal.sh > trimal.log 2>&1 &
for i in `cat stop_codon_mafft.txt`;do
name=${i%.*}
grep "X" trimal_out/${name}.trimal.fa
done # 查看修剪后序列文件是否还存在X,结果:部分存在,部分不存在


# Step 3 构建单基因树,这里是同源基因都构建了进化树,后面绘制DensiTree只使用单拷贝直系同源基因树
#WD: ~/Sparus_latus/phylogenomics/Aligniment_Genome/homology/OMA/orthofinder_results/Results_Jan01mkdir AA_IQtree_tree
for i in `ls trimal_out/*fa`;do
name=`basename $i .fa`
echo "$IQtree -s $i -m TESTONLY -B 1000 -c AUTO --prefix AA_IQtree_tree/${name}_AA --seed 1" >> AA.IQtree.shdone # AA.IQtree.sh脚本
nohup bash AA.IQtree.sh >AA.IQtree.log 2>&1 & # 也可以一开始只使用Single_Copy_Orthologue_Sequences中序列构建单拷贝直系同源基因树
## 只保留具有contrees的gene/windows tree
ll *contree | awk '{print $9}'|awk -v FS='.' '{print $1"."$2"."$3"."$4"."$5".treefile"}' > usetreename.txt
perl fa/filter_windows_tree.pl fa usetreename.txt # 所有需要的进化树保留在tree目录中
* windows tree,对每条染色体,按固定距离滑动窗口提取序列(如10kb,或100kb,根据自己数据选择窗口大小),然后针对每个窗口序列构建进化树。


图11
四、绘制DenisTree
# step1 根据拓扑关系过滤树,使用newick_utils软件
### 满足(A,B)
mkdir mono
for tre in `ls tree`;do
   nw_clade -m tree/${tre} A B > mono/${tre} &
done # 进化树中A,B不为姐妹群,则输出文件为0K
cd mono
find . -size 0k -exec rm {} \; # 删除0k文件
ll *tre | awk '{print $9}' > ../mono.txt
cd ..
perl treefilter.pl tree mono.txt rmtre.txt # 保留A,B为姐妹群的基因树
# Step2 合并单基因树到一个文件中
cd tree
find . -type f -name "*treefile" |xargs sed 'a\' >all_genes.tre
# Step3 对基因树进行reroot和删除不需要的物种,也可以不进行这一步
nw_reroot all_genes.tre Loc > all_genes_reoot.tre # Loc为定根的物种名
nw_prune all_genes_reoot.tre Ipu Ola Ame > DensiTree/all_genes_reoot.tre_DensiTree.tre # 删除Ipu、Ola和Ame
# Step4 使用Figtree 删除bootstrap值,进化树之间必须有换行符
## 绘制DensiTree,进化树必须只有拓扑和分枝长度数据
# Step5 使用DensiTree.v2.2.5.jar
绘制DensiTree,见图12-17## 完整使用细节,看软件使用说明书


图12 DensiTree软件


图13 载入数据




图14 更改展示形式和细节调整



图15 调节分支顺序


图16 调节分支顺序



图17 更改DensiTree线条颜色。
DensiTree软件还有很多功能,这里就不一一展示了,想要用的可以去查看软件的说明文档。还可以跟地理分布图一起展示。

这篇教程只是写绘制DensiTree的流程,写的时候,想着数据的准备也应该写出来,但是这样内容太多了。所以有些脚本的代码没有写出来,只是展示一个数据准备的流程,大家可以照着做。有什么疑问或者需要某个软件/脚本,可以在后台留言脚本和软件名称或者点击“沟通交流”加好友,我看到了会回复的。   



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
努力奋斗
回复

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
30
注册时间
2017.4.3
在线时间
117 小时

 楼主| 发表于 2021.3.24 21:40:20 | 显示全部楼层
我真的心太累了,直接从微信公众号上上传,最后代码全连在一起了,一个个手动分开。什么时候这个问题能解决啊!
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
69 小时

发表于 2021.3.24 22:10:41 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
72 小时

发表于 2021.4.6 08:32:18 | 显示全部楼层
加油,加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2021.1.28
在线时间
12 小时

发表于 2021.4.10 11:07:01 | 显示全部楼层
打卡!!!!!!!!!!!!
回复

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
30
注册时间
2017.4.3
在线时间
117 小时

 楼主| 发表于 2021.4.12 13:41:18 | 显示全部楼层
努力奋斗
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
2
注册时间
2018.3.15
在线时间
21 小时

发表于 2021.4.14 23:34:54 | 显示全部楼层
谢谢~~~
回复

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
30
注册时间
2017.4.3
在线时间
117 小时

 楼主| 发表于 2021.4.17 20:14:54 | 显示全部楼层
努力奋斗
回复

使用道具 举报

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
210
注册时间
2017.7.3
在线时间
641 小时

活跃会员荣誉管理


发表于 2021.4.26 19:30:46 | 显示全部楼层
静默人声412 发表于 2021.3.24 21:40
我真的心太累了,直接从微信公众号上上传,最后代码全连在一起了,一个个手动分开。什么时候这个问题能解决 ...

在编辑窗口有获取公众号文章按钮
新的一天加油!
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
2
注册时间
2018.3.15
在线时间
21 小时

发表于 2021.5.6 18:46:39 | 显示全部楼层
赞赞赞~
回复

使用道具 举报

帝王蝶

Rank: 4

主题
5
注册时间
2016.4.7
在线时间
53 小时

发表于 2021.6.3 15:19:55 | 显示全部楼层
强大
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
90 小时

发表于 2021.7.5 11:29:01 | 显示全部楼层
不明觉厉
差点忘记签到了啊
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2021.7.12
在线时间
1 小时

发表于 2021.7.27 11:29:46 | 显示全部楼层
.............
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2019.6.21
在线时间
64 小时

发表于 2021.7.28 09:02:09 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2020.3.24
在线时间
6 小时

发表于 2021.8.21 16:31:56 | 显示全部楼层
okkkk
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.9.7
在线时间
3 小时

发表于 2021.9.10 11:15:29 | 显示全部楼层
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表