一、 写在前面
2023年4月中旬自己开始做基因家族的分析,对于这块自己没有接触过,因此也是一个挑战,没事!!!(安慰自己),对于基因家族的分析网上的教程很多,跟着步骤走就可以。在这部分,我自己主要是做生信这块,实验验证是师姐在做,所以论文结构自己不用操心。此外,可视化的工具很多,也很方便,不需要自己特意去学。我们这里就60%使用TBtools软件进行可视化和分析。
此外,本次分析80%的内容都是基于TBtools。确实牛X!!自己开始接触TBtools是在2019年吧,也是通过一个师兄的推荐才知道的。2019年CJ还没将TBtools发表在MP上,那时还是预印版本吧。但是,引用已经有了很多,了不起哦。后面TBtools一直在开发新的的小“软件” or“程序”,将生信分析的门槛一降再降。点赞点赞!!!
--Du
如想获得本文文档,可看文末!!
注意:此教程有些话语可能会带有自己的方言,读不通时也不要在意!![泪目!]
一,在Pfam数据中获得基因家族
我们这里预测作物中某一个基因家族的基因,目前在此作物中未报道。因此,使用Pfam数据库中一致的基因进行同源搜索(其实,你也可以使用已知作物中的基因进行同源搜索,获得结果基本一致)。那么我们就根据文章中和报道的Pfam数据库中的基因作为基序,进行同源搜索。
-
在Pfam数据库中下载FBNs基因家族(Pfam 04755),Pfam网址:https://pfam-legacy.xfam.org/
-
打开网址:http://www.ebi.ac.uk/interpro/entry/pfam/?search=0477#table
- 点击进入
PF04775
,下载所有的Proteins序列
以上只是其中的一种方法,但为获得FBN基因家族的蛋白序列。下面使用Pfam数据中搜索
- 打开网页。https://pfam-legacy.xfam.org/
-
搜索
-
进入
- 搜索后获得PAP_fibrillin
下载Reviewed
的PF04755序列
二、同源序列检索预测
对于同源基因的搜索,很多基因家族的文章都使用HMMER进行检索,也有一些文章是使用BLAST。你任选其中一个即可,都能获得你想要的结果同源基因。在做分析的时候,我将使用Hmmer寻找同源基因的文章分享在公众号中,在评论区有一个大佬对HMMER和BLAST之间的差异给出回答。
这两个方法原理上区别,balstp是基于序列同源性进行打分的,有打分矩阵,hmm是基于隐马尔可夫模型,对序列结构域进行比对。来自“泼皮混混”的评论。
2.1 HMMER同源结构域搜索
2.1.1 Hmmer的安装
安装,主要是使用源码安装或是是使用conda
进行安装即可。
- conda安装
conda install -y hmmer
-
源码安装:
官网:http://www.hmmer.org/
任意下载一个版本即可,安装步骤不再做说明。
2.1.2 使用hmmbuild构建.hmm文件
在有些数据库中是有.hmm
文件,只需要下载即可。但是,这仅仅只限于有些大数据库。对于我们自己使用,不可能全部都有,这就需要我们自己构建,很多教程到这步就是让你收费了.......。
在本教程,讲述其中一种方法吧,希望对大家有所帮助。
hmmbuild构建时,需要使用
.sto
文件进行构建。因此,我们必须获得.sto
文件。
- 使用mafft软件进行间序列进行对齐
mafft --auto --clustalout ../Pfam_PF04755_reviewed.fasta > Hmmbuild_index/Pfam.FBNs.align.clustal
转换:
http://sequenceconversion.bugaco.com/converter/biology/sequences/fasta_to_phylip.php
- hmmbuild构建文件
hmmbuild Pfam.FBNs.hmm sample.stockholm
- hmmsearch
hmmsearch Hmmbuild_index/Pfam.FBNs.hmm Potato/DM_1-3_516_R44_potato.v6.1.working_models.pep.fa > ../02_Result/Potao.hmmer.out.txt
- 筛选出最佳的结果,E-value值小于
1e-5
,Score值大于“> 90” - 对于筛选结果,可以直接使用Hmmsearch获得结果;也可以如上所示根据自己需求进行筛选,自己做的话,如果搜索的目的基因太多,而自己不需要这么多的同源基因,自己会进行手动过滤一些同源性较弱的基因。
cat Potato.hmmer.out.txt |grep -v "#" | awk '{if($4 < 1e-5 && $5 > 90) print $9}' | sort | uniq | grep -v "+" > Potato.hmmer.best.out.txt
2.2 提取目的基因序列
日志:通过Hmmsearch获得同源基因的ID,那么后面对目的同源基因进行进化树、结构域、motif等的分析,这些分析都需使用目的同源基因的序列。
如何获得同源基因序列??
- 使用脚本获得
- 使用ggffead获得,需要获得同源基因的
.gtf
文件等信息。 - 生信工具获得、如TBtools等。
对于这步、我们就多做讲解,使用自己拿手的方式获得即可。
问:后面的分析使用核酸序列 or蛋白序列呢??
答:都可以。
FBN 家族的分析日志。使用Pfam、拟南芥(11)和水稻的FBN家族基因同源搜索马铃薯中的FBN同源基因
## 水稻中的FBN家族基因
cat all.pep | grep ">" | grep fibrillin |awk -F "|" '{print $1}' | awk -F " " '{print $1}' | sed 's/>//g' > O_sativa.FBN.id.txt
##拟南芥中FBN家族基因
可以在拟南芥网址中的同源搜索,也可以在拟南芥蛋白数据中搜索
cat Araport11_pep_20220914 | grep FBN | awk -F "|" '{print $1}' | sed 's/>//g' > Araport11_FBN.id
2.3 使用TBtools提取目的基因
说实话,TBtools确实是个很牛的生信工具,基本可以让你不写代码获得你想要的东西。以及,各种类型的小脚本软件都一直在开发。赞赞!!
2.3.1 TBtools软件的下载
- 网址:https://github.com/CJ-Chen/TBtools
- 安装。
- 动手运行
2.3.2 提取序列
- 准备作物所有的蛋白序列文件(or基因文件)
-
目的基因的ID
- 打开TBtools,
Fasta Extract or Filter (Qyick)
-
获得结果
2.4 目的同源基因motif分析
2.4.1 使用MEME进行motif预测
-
上传相关的fa文件,以及修改相关的参数,进行提交
-
输出结果
输出结果很快,有以下几个结果文件。
2.4.2 motif可视化
对于motifi分析可以参考一下文章:
- TBtools | 多图合一至强版教程!进化树 + Motifs + 结构域 + 启动子 + 基因结构 + ....,TBtools开发者本人的教程
- TBtools | 基因家族分析 (进化树、Motifs、结构域)
- 或是本篇教程
MEME网址结果可以给我们的seqlogo信息和motif信息。
-
Seqlogo
结果文件中就有seqlogo文件信息。
也可以自己的下载后绘制。
按以下操作即可下载序列。
也可以下载已有的seqlogo图片。
下载后所有的motif序列信息。
- 使用R语言对Seqlogo序列进行可视化
这里借用这篇教程,基因结构及motif分析。批量生产Seqlogo可视化。
我们可以根据自己的motifi数量进行命名,我自己只有10个motif信息。所以命名为motif1-10.txt
。
## 加载所需要的包
library(ggplot2)
#BiocManager::install("ggseqlogo")
library(ggseqlogo)
## 批量生产文件名
filelist = c(paste0('motif',1:10,'.txt'))
filelen <- length(filelist)
##批量读取
data.list <- list()
for (i in 1:filelen) {
data.list[[paste0('motif',i)]]=scan(filelist[i],what = '')
}
ggseqlogo(data.list,col_scheme="clustalx", ncol = 5)+
theme(axis.line = element_line(colour = 'black'),
axis.text.x = element_blank(),
legend.title = element_blank())
ggplot()+
geom_logo(data.list, col_scheme = "clustalx")+
theme_logo()+
facet_wrap(~seq_group,ncol = 5,scales = "free_x")+
theme(axis.line = element_line(colour = 'black'),
axis.text.x = element_blank())
对比一下MEME网站中的图形。
对于Seqlogo的绘制,美化,可以根据很多优秀的教程。在网上上一搜,都可以找到。
2.4.3 motif的分析
-
下载结果文件
MAST XML output
,使用TBtools软件进行可视化。
-
打开TBtools中的
Gene Structure View
,只需上传MEME中的XML文件即可,上传上去直接点击Start
。
--
操作:
结果:
保存!!
注意:我们这里保存的时候最好保存为PDF或SVG格式,输出为矢量图。
如果我们的教程只是到这里,那么就没有什么意义了。因为,类似非常优秀和详细的教程很多。绘制出图形是一方面、美化可是重头戏。
在MEME输出文件中,也提供了motif的图形,也可直接使用。
2.5 基因家族保守结构域分析
- 使用
Batch CD-Search
进行预测,网址:https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi
-
提交序列信息即可
- Batch CD-search只支持目的基因蛋白序列信息, 以及序列数量少于1000。
Warning: Batch CD-Search accepts only protein sequences. The maximal number of query sequences per request is 1000. A single query sequence can not exceed a length of 40,000 residues.
可以提供你的邮箱,等运行结束后,直接发送到你的邮箱。如果序列较多,建议提供邮箱。
-
下载文件
结果文件:
- 打开TBtools中的
Visualize NCBI CDD DOmainPattern
-
输入结果文件和fa文件
-
根据自己的需求进行调整即可。
-
输出文件。
2.6 进化树分析
进化树分析,在基因家族中是必须的,以及在很多图中都是需要的。进化树分析和绘制,也有很多教程,参考iqtree+ggtree绘制进化树教程、或是你也可以使用MEGA来做分析。
2.6.1 iqtree+ggtree绘制进化树教程
- iqtree获得树文件
所需软件
- mafft
- iqtree
mafft安装
我是使用服务器中运行的,安装可以使用conda
- iqtree
conda install mafft
iqtree官网
http://www.iqtree.org/
iqtree功能很强大,大家可以查看软件的官方文档。
安装
conda install iqtree
软件安装好后直接运行即可。
- 序列准备
进化树序列可以使用蛋白序列或核酸序列即可,格式按其准备即可。
>B2LU34
MTSIAFWNAFTVNPFPAAARRSPPPLTPFTSGALSPARKPRILEISHPRTLPSFRVQAIAEDEWESEKKALKGVVGSVAL
AEDETTGADLVVSDLKKKLIDQLFGTDRGLKATSETRAEVNELITQLEAKNPNPAPTEALSLLNGRWILAYTSFAGLFPL
LGAESLQQLLKVDEISQTIDSEGFTVQNSVRFVGPFSSTSVTTNAKFEVRSPKRVQIKFEEGIIGTPQLTDSIVIPDKFE
FFGQNIDLSPFKGVISSLQDTASSVAKTISSQPPIKFPISNSNAQSWLLTTYLDDELRISRADGGSVFVLIKEGSPLLT
>B4F6G1
MTSIAFCNAFTVNPFLAAARRSPPPLTPLTSVALSPARKPRILAIFHPRTFPSFRVQAIAEDEWESEKKTLKGVVGSVAL
AEDEKTGADLVVSDLKKKLIDQLFGTDRGLKATSETRAEVNELITQLEAKNPNPAPTEALSLLNGKWILAYTSFVGLFPL
LGAESLQQLLKVDEISQTIDSEGFTVQNSVRFVGPFSSTSVTTNAKFEVRSPKRVQIKFEEGIIGTPQLTDSIVIPDKVE
FFGQNIDLSPFKGVISSLQDTASSVAKTISSQPPIKFPISNSNAQSWLLTTYLDDELRISRADGGSVFVLILESSPLLT
>O49629
MATVQLSTQFSCQTRVSISPNSKSISKPPFLVPVTSIIHRPMISTGGIAVSPRRVFKVRATDTGEIGSALLAAEEAIEDV
EETERLKRSLVDSLYGTDRGLSASSETRAEIGDLITQLESKNPTPAPTEALFLLNGKWILAYTSFVNLFPLLSRGIVPLI
KVDEISQTIDSDNFTVQNSVRFAGPLGTNSISTNAKFEIRSPKRVQIKFEQGVIGTPQLTDSIEIPEYVEVLGQKIDLNP
IRGLLTSVQDTASSVARTISSQPPLKFSLPADNAQSWLLTTYLDKDIRISRGDGGSVFVLIKEGSPLLNP
- mafft比对
使用mafft将序列对齐。
mafft test.fa > test.aligend.fa
我们获得对齐后的数据格式。
- iqtree构建树
iqtree -s test.aligend.fa -m MFP -bnni -nt AUTO -cmax 15 -redo -bb 1000
关于iqtree的使用,可以看这篇教程IQ-TREE的使用 - 超快速用极大似然法构建进化树,讲的很详细。
必须参数:
-s 输入多序列比对文件
-nt 多线程,AUTO是自动多线程
-bb 1000 指定了要用快速BS法做1000次
最终,我们可以获得以下结果文件。
- ggtree绘制进化树
这里,我们使用基迪奥的教程,如何绘制添加分类色块的进化树?,这个教程也是讲解得很详细。
注意:我们这里使用iqtree输出文件test.aligend.fa.treefile
作为输入文件。
#载入相关的R包;
library(ggtree)
library(treeio)
library(ggplot2)
#读入newick格式的进化树文件;
tr = read.newick("test.aligend.fa.treefile")
ggtree(tr)
#为进化树添加叶标签;
p1 <- p0 + geom_tiplab(size=2,color="grey10")
p1
#为进化树添加圆形顶点;
p2 <- p0+ geom_tiplab(size=2,offset=0.03, color="grey10")+
geom_tippoint(color="#6bc72b",fill="#6bc72b",
alpha=0.4, size=3,shape=21)
p2
后面的教程参数调整,按着教程即可如何绘制添加分类色块的进化树?
2.6.2 MEGA制作进化树
此部分内容来自:TBtools | 基因家族分析 (进化树、Motifs、结构域)
输入数据为目标基因家族的蛋白质序列。
先进行多序列比对,用MUSCLE默认参数。
图片将比对好的结果保存为.meg格式。
重新打开比对后的文件,构建进化树,使用最大似然法,根据需要选择建树方法。再构建之前可以进行模型的预测,这里节省时间直接使用默认参数。
现在就构建好了一棵进化树,导出为.nwk格式。接下来最后一步就是再TBtools中展示所有结果。
2.7 使用Figtree绘制进化树
在2.6
和2.7
小节中,我们讲述了使用ggtree和MEAG绘制进化树,这些软件都是比较常用的。在这次作图过程中,自己的无意间也查询到使用Figtree
可视化工具绘制进化树。主要是看到这张图,平时自己看到的图都是矩阵类型或是圆形,类似这个半圆看着是比较好看。
Figtree网址:http://tree.bio.ed.ac.uk/software/figtree/
软件下载可以到GitHub中下载:https://github.com/rambaut/figtree/releases
下载后无需安装,即可使用(根据自己的版本调整)。
将
FigTree v1.4.4
快捷键发送到桌面即可对于Figtree软件的使用,全网依旧是很一定数量的教程,大家可以自行进行查找,或观看帮助文档。
2.7.1 Figtree绘制进化树基础图形
打开Figtree界面是比较简单,这个软件的获得的图形的类型也是相对比较少,只适合小众类型的进化树绘制。对于很复杂类型进化树还是不推荐使用Figtree绘制。
- 点击
File
-Open
,导入数据
-
获得进化树
-
调整。全部参数可以在左侧调整即可。包括,大小、间距、距离参数等。
以上参数,仅仅只是必要调整的参数,具体看自己的分析进行调整即可,无固定模式。
2.7.2 Figtree绘图的模式
我在前面说过Figtree绘制进化树的图类型很少,只有三种大类型
。具体如下所示。
-
一般的聚类类型
- 圆形
circular
2.7.3 Figtree绘制进化树美化图形
如何进行美化,是我们一直在追求的方向。在进化树中分支的上色是必须的,在Figtree中依旧可以做。注意:我们这里只是简单的说明如何上色,具体操作自己进行。
最终图形可以获得如下图所示。
2.7.4 Figtree导出图形
调整好图形参数,如何导出图形呢?操作如下所示。File
-Export JPEG/PNG/PDF.....
,导出适合的的图形格式即可,但是建议导出的矢量图。后期AI进行调整。(通过上面导出图形,我们可以看到图形的颜色长度是不同的,这个问题要如何解决,暂时没有找到好的方法。在ggtree绘制中自己也遇到这里的问题。如果在的图形软件中无法解决,只能通过后期解决。)
2.7.4 重新文章中图形
那么如何绘制类似的图形呢?根据前期的参数,只需要进一步优化即可。
- 主图
(1) 将图形性状选择圆形
(2) 调整Root Angle
和Angle Rangle
调整到适合的形状。
-
分类附图
在这个图中,我个人将其进化树分为进化树分类附图。这个图也是使用的Figtree进行绘制。具体操作如下所示。
-
选择分类图形
-
调整参数
-
树枝的宽度可以宽1-2个size
- 调整自己喜欢的
Trabsform Branches
-
继续调整
--
注意:
进化树的分支,主图和附图要一致。为了进一步确定明确两个图的一致性,建议直接在附图中,对分支进行填充颜色。操作与上述一致。
2.7.5 AI合并美化
- 打开AI
-
新建图形
-
导入进化树图形
-
Ctrl + R
打开AI中的标尺、拖出x轴或Y轴参考线
-
调整半圆进化树,做到“横平竖直”
-
Ctrl + A
全选,选择图形,Ctrl + C
进行复制,或直接进行拖拽到新建图形中。 - 调整适合的图形大小,调整时,一直按住
shuft
,避免图形横纵大小改变。
-
建议,在图形中如有新的图形产出,建议每个新的图形都新建立一个图层,利于后期的修改。
-
随后就进行进化的调整,我们在这里,需要对AI有一定的基础知识,才可以。比如,如何随意修改图形的形状,类似图例所示。这里操作很繁琐,具体操作自己进行。
-
导入进化树分支
-
如何线条太细,可以进行调整适合粗细。
-
分支添加颜色
-
新建图形
-
选择椭圆工具
-
绘制椭圆,调整适合的分支位置和的添加分支颜色
-
适当的调整颜色
-
依次绘绘制即可
-
字体调整(如果在图形中梯子较小,也可以在AI中调整)
使用选择工具,选择调整字体,直接进行修改即可。
-
调整图形大小
-
最终出图
-
也可以直接间监矩形进化树进行进行合并,相比育德圆形或半圆,调整颜色柱就很容易,直接拉成一样长度即可。
--
细节自己调整。
2.8 目的基因结构可视化
需要文件:
- 目的基因注释文件(GFF or GTF)
- 进化树文件(可选)
2.8.1 使用ID和基因组注释文件绘制
- 使用TBtools直接操作,依次点击:
Gene Structure View
结果如图所示:
2.8.2 提取目的基因的注释文件(推荐)
我们会发现,输入ID
处也是可以输入进化树文件信息。因此,我们推荐直接提取获得目的基因的注释文件信息,单独使用GTF文件信息或是GFF信息进行绘制。
- 获得GFF注释信息
使用已有的目的基因的ID与基因组注释文件进行匹配获得。
cat Araport11_GTF_genes_transposons.current.gtf | grep -wf TAR11.test.id > TAR11.test.gtf
$ cat Araport11_GTF_genes_transposons.current.gtf | grep -wf TAR11.test.id | head
Chr1 Araport11 mRNA 18935301 18937665 . + . transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18935380 18935673 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18935743 18935796 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18935908 18935982 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18936083 18936205 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18936278 18936469 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18936552 18936635 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18936723 18936815 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18936903 18936956 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
Chr1 Araport11 CDS 18937039 18937118 . + 0 transcript_id "AT1G51110.1"; gene_id "AT1G51110";
-
进化树获得
同上的方法获得 -
MEMExml or MAST.xml文件
同上 -
绘图
依次提交相关的文件即可
2.9 进化树、Motifs、结构域、基因结构合图绘制
以上的操作,都可以获得单张图形,那么如何多图绘制在一起呢?TBtools也提供了相关的教程,TBtools | 多图合一至强版教程!进化树 + Motifs + 结构域 + 启动子 + 基因结构 + ....,我们可以根据此教程进操作。具体如下:
获得结果(来自CJ教程):
2.10 图形美化
到这里,我们的整张图形就可以获得。但是,只是这样的话,我觉得自己的这个教程就没有意义。我前面说过,我的这个教程重点是图形美化。自己是更喜欢,TBtools单张出图的类型,然后进行AI或PS美化的。软件默认的颜色,我自己不是很喜欢,但是也可以自己调整,也是很方便的哦。
2.10.1 TBtools图形颜色的调整
我们这里只是随意进行调整,图形无任何意义。
-
步骤一、点击图形中的方块、右键
- 调整色块
3、选择先要的色块、点击Selecteed
4、更改成功
但是你会大发现,图中所有一样的颜色色块都会改变。
类似的功能、自己逐渐去摸索。
2.10.2 单张出图
如果上面的方式没有很好实现自己想要的效果。那么,我们就只能单张出图、后面再进行合并。
注意:在绘图时,我们的要提前想好自己的文章或这张图的颜色设置,以及图形的色调是属于什么类型的。理论上,一整篇文章图形色调和类型要保持一致。
如果,在后期的调整中。图形颜色需要重新调整,我们可使用AI进行调整或是重新绘制,少量还是比较方便,但是图形又大有多,重画是很奔溃的事情。
三、IA图形美化
美化,我罗列出单个章节进行讲解。表明,是很重要的。以及,图形的美化,需要不断学习和模范大牌期刊的图形类型,以及自己要时刻进行总结和创新。对于创新,这个就比较玄学,每个人的审美不同,逻辑不同,关注点不同.......导致最终看到的点也不同。因此,我们在不是很离谱的创作中,结合自己的审美进行美化即可。我们要坚信:审美。首先要符合自己,其次,再考虑别人。只有自己先认同,你才有可能让其他人也认同!
3.1 使用工具
1.推荐使用的工具:AI、PS
如果不知道类似软件的,自己百度。
- 如何安装
- 有钱人:购买正版
- 穷人(和我一样):薅羊毛,使用破解版
- 如何获取安装包
在本公众号中回复关键词获得。
- PS安装包关键词:
PS
- AI安装包关键词:
AI
或是你自己寻找相关版本的安装包即可。
提示:请自己输入正确的关键词(每次看到有些同学们的关键词,真的很无语......)
3.2 实际操作
- 打开AI,新建图层
A4
-
导入进化树,适当调整进化树的宽度和字体大小
- 依次导入的目的基因的motif、基因结构域等图形。并依次按进化树基因名进行排序即可。
- 为后期的图形的整齐性,我们使用参考线进行对齐,便于后期的调整。
注意:这里看到我们的motif的图形颜色很难看,这就是前期没有考虑颜色的结果。因此,我一直强调,文章图形颜色统一的重要性,图形颜色搭配合理,你的论文已经成功1/3了。
换一种颜色就感觉好多了呀。 -
添加基因结构图
添加图形的操作都是一样的,不做多赘述。
- 如何美化
对于美化,每个人的要求不一致,只要符合你的审美即可。我们在这里就直接添加渐变色。 -
新建一个图层
新建图层置于最底层。
-
选择图形工具
-
利用进化树的分支,将其进行分类
-
填充颜色(根据自己的喜好)
-
更改透明图
-
渐变色
- 不透明度:60
-
中间位置:10-50%
结合实际情况调节。
-
最后图形
图形很多细节需要自己耐心调节,这里只是做示范,相对比较粗糙。
四、多物种共线分析
共线分析依旧是使用TBtools,哈哈哈哈,做基因家族TBtools可以帮你完成80%的生信分析。毫不夸张!!!!!
TBtools共线分析的教程很多,我们以零基础多物种间共线性分析教程作为参考(也不是作为参考了,是直接按他的步骤进行操作)。其他参考教程:全基因组共线性分析、无限个!物种共线性分析结果可视化、任何人!一键完成物种间的共线性分析与可视化。
4.1 需要文件
- 参考基因组fa文件
-
注释文件GFF or GTF
TBtools可以对无限个作物进行共线分析,牛!!!
4.2 染色体统一命名
在这个教程中,有这样的一个步骤,如果你需要,你就进行操作。
-
gtf文件进行ID prefix
-
fa文件进行ID prefix
4.3 实操
- 打开
one step MCScanx
小程序
-
输入两个作物的文件信息
- 点击开始
Start
- 如果是多个作物,那么依次进行两两比较。比如:共线结果是以这样的顺序:Tomato-LA-Arabidopsis
比对顺序:
- Tomato-LA
- LA-Arabidopsis
- 比对结果GFF文件合并
- 打开
Text Merge for MCScanX
程序
合并多个的MCScanX的结果文件中的GFF文件
拖拽文件
6.比对结果ChrLayout.tab.xls
文件合并
- 比对结果
geneLinks.tab.xls
文件合并
同上操作! -
合并文件
最终获得以下3个文件,用于绘制图形。
- 要在共线中显色的基因ID
Solyc03g062790.3.1
Solyc10g018590.2.1
Solyc01g104320.4.1
Solyc03g083420.4.1
AT4G22240.1
AT2G35490.1
AT1G51110.1
AT5G53450.3
........
- 绘图。打开
Multiple synteny plot
输入参数
输出图形
注意,在输出图形中,我们可以看到作物染色体位置是有改变的。那么,如何更改呢?回答:直接更改Chr文件即可。
更改这里的顺序即可!
五、同源目标基因元件预测
目标基因的元件预测,我们这里主要介绍使用两个网站进行。
5.1 提取目标基因上游2000bp
参考教程顺式作用元件预测和新的可视化方式、植物启动子-顺式作用元件-批量提取-预测-可视化分析,同样是使用TBtools操作。
- 需要文件
- 作物参考基因组fa文件
- 注释文件GFF or GTF
- 目标基因ID
-
直接使用TBtools中的
Gtf /Gff3 Sequences Extractor
获得每个基因的fa序列
输出文件
点击Initalize
,选择CDS
选择上游2000bp的fa序列
-
目标基因的fa序列,打开
Fasta Extract or Filter (Quick)
输出结果文件:
-
查看信息是否正确,打开
Fasta Stats
-
转换序列(全部为大写),打开
Sequence Manipulate (Rev&Comp
5.2 提交预测网址进行顺式作用预测
预测,这里使用两个网站进行预测,分别是PlantCare
和PLCAE
。
5.2.1 使用Plantcare进行预测
网址:http://bioinformatics.psb.ugent.be/webtools/plantcare/html/
-
上传序列后,Plant可以提供你自己的邮箱,运行结束后,结果直接发送到你的邮箱中。
-
邮箱中获得结果,根据你的序列多少,10分钟以上吧!
-
结果
-
使用execl打开后
1. 基因ID;
2. 顺式作用元件名称;
3. 顺式作用元件序列;
4. 顺式作用元件的起始位置;
5. 顺式作用元件的长度;
6. 顺式作用元件所在的链的方向;
7. 物种名;
8. 顺式作用元件所在的功能分类;
删除某些不需要的结果:
需要删除:
1. 剔除第2列为空的行
2. 剔除第2列为unnamed的行
3. 最后一列,无功能作用的
具体删除的数据,根据自己的分析来做。
最后,可以删除掉1000行以内
--
来自顺式作用元件预测和新的可视化方式,这个意见有重要的参考意义。如果不合并,导致元件的作用太多,绘制出的图形颜色太杂,且不好看。
-
绘图
绘图前还需要准备基因的长度文件
输入数据,设置参数
结果:
在TBtools中也可以输入进化树文件。
我们这里也可以使用的起那么AI中的呢进化树进行模板进行美化。
5.2.2 PLACE进行预测
网址:https://www.dna.affrc.go.jp/PLACE/?action=newplace
- 缺点:PLACE一次最大只能输入20条基因序列,有一定的限制性。获得结果为网页版,如要整理,只能手动整理或使用脚本进行整理。
- 优点:速度快!
-
获得结果
每个基因为单独的,需要自己整理。
- 只给元件名称、开始位置、序列、功能(SITE,需要点击进去才可以看到)
- 整理,单独粘贴复制到execl中,并使用脚本进行整理。
选择哪个网站进行预测,取决于自己。只要结果符合我们自己的预期结果即可!!!
5.2.3 热图可视化
输入数据格式如下(可以根据自己的情况筛选):
脚本:
install.packages('tidyverse')
intall.packages('RColorBrewer')
# 加载包
library(tidyverse)
library(RColorBrewer)
# 1.读取数据
df <- read_tsv('data.txt', col_names = F) %>% select(1,2)
# 2.整理数据
tidy <- df %>%
group_by(X1, X2) %>%
summarise(number = n()) %>%
arrange(desc(number))
# 3.查看数量分布,确定配色个数
summary(tidy$number)
# 最大值为9,所以下面的代码 hcl.colors(9, "RdYlGn")中为9
# 4.画图
ggplot(tidy, aes(x = X2, y = X1, fill = number)) +
geom_tile(color = 'black') +
geom_text(aes(label = number),col='black',cex = 1.5) +
scale_fill_gradientn(colors = rev(hcl.colors(9, "RdYlGn"))) +
scale_x_discrete(position = "top")+
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 0),
axis.title = element_blank(),
axis.text = element_text(size = 7, color = 'black'))
# 通过修改 scale_fill_gradientn参数给每一个值指定颜色
cc <- c('#d9d9d9', '#f7fcb9', '#d9f0a3', '#addd8e', '#78c679', '#feb24c', '#fd8d3c', '#fc4e2a', '#b10026')
ggplot(tidy, aes(x = X2, y = X1, fill = number)) +
geom_tile(color = 'black') +
geom_text(aes(label = number),col='black',cex = 2.5) +
scale_fill_gradientn(colors = cc) +
scale_x_discrete(position = "top")+
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 0),
axis.title = element_blank(),
axis.text = element_text(size = 7, color = 'black'))
5.2.4 美化
基于AI进行美化,方法同上
六 ENDING
说实话,基因家族的文章分析确实消耗的时间和精力不算是很多。生信部分就差不多这些吧!再加上一些组学的数据来验证即可。除了生信的部分,剩余就是实验来验证,将两者进行结合,好一点的文章也可以发。我自己前面没有接触过基因家族的分析,因此,本次就是现学现做,做的还是比较简单。
本次来接触基因家族的分析,感触最深的就是,TBtools真的很强大。基因家族的分析、画图都可以使用它来完成。不得了啊,真的是将做生信的门槛一降再降,点赞点赞
本期内容是自己的做了一个整理,算是“教程搬运工”,也是自己在做分析后做的总结。自己不知道,这次分析后,多久以后还能涉及基因家族的分析。总结总结!! 但是,说实话!这个总结也花费自己很长的时间,如果你想获得这个教程的文本文档,可以“喜欢点赞,支持”,我在后台看到后会第一时间将文档链接发给你!!
小杜的生信筆記,主要发表或收录生物信息学的教程,以及基于R的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!
网友评论