基于TBtools做基因家族分析 | 生信部分






基于TBtools做基因家族分析 | 生信部分

BioinfoDu  生信果  2023-05-10 19:01:07

一、 写在前面

2023年4月中旬自己开始做基因家族的分析,对于这块自己没有接触过,因此也是一个挑战,没事!!!(安慰自己),对于基因家族的分析网上的教程很多,跟着步骤走就可以。在这部分,我自己主要是做生信这块,实验验证是师姐在做,所以论文结构自己不用操心。此外,可视化的工具很多,也很方便,不需要自己特意去学。我们这里就80%使用TBtools软件进行可视化和分析。

此外,本次分析80%的内容都是基于TBtools。确实牛X!!自己开始接触TBtools是在2019年吧,也是通过一个师兄的推荐才知道的。2019年CJ还没将TBtools发表在MP上,那时还是预印版本吧。但是,引用已经有了很多,了不起哦。后面TBtools一直在开发新的的小“软件” or“程序”,将生信分析的门槛一降再降。点赞点赞!!!

–Du


「如想获得本文文档,可看文末!!」


「注意:此教程有些话语可能会带有自己的方言,读不通时也不要在意!![泪目!]」



一,在Pfam数据中获得基因家族

我们这里预测作物中某一个基因家族的基因,目前在此作物中未报道。因此,使用Pfam数据库中一致的基因进行同源搜索(其实,你也可以使用已知作物中的基因进行同源搜索,获得结果基本一致)。那么我们就根据文章中和报道的Pfam数据库中的基因作为基序,进行同源搜索。

  1. 在Pfam数据库中下载FBNs基因家族(Pfam 04755),Pfam网址:https://pfam-legacy.xfam.org/

  2. 打开网址:http://www.ebi.ac.uk/interpro/entry/pfam/?search=0477#table

  • 点击进入PF04775,下载所有的Proteins序列

以上只是其中的一种方法,但为获得FBN基因家族的蛋白序列。下面使用Pfam数据中搜索

  1. 打开网页。https://pfam-legacy.xfam.org/
  2. 搜索
  3. 进入
  4. 搜索后获得PAP_fibrillin下载Reviewed的PF04755序列

二、同源序列检索预测

对于同源基因的搜索,很多基因家族的文章都使用HMMER进行检索,也有一些文章是使用BLAST。你任选其中一个即可,都能获得你想要的结果同源基因。在做分析的时候,我将使用Hmmer寻找同源基因的文章分享在公众号中,在评论区有一个大佬对HMMER和BLAST之间的差异给出回答。

这两个方法原理上区别,balstp是基于序列同源性进行打分的,有打分矩阵,hmm是基于隐马尔可夫模型,对序列结构域进行比对。「来自“泼皮混混”的评论。」

2.1 HMMER同源结构域搜索

2.1.1 Hmmer的安装

安装,主要是使用源码安装或是是使用conda进行安装即可。

  1. 「conda安装」
conda install -y hmmer
  1. 「源码安装:」「官网」:http://www.hmmer.org/任意下载一个版本即可,安装步骤不再做说明。

2.1.2 使用hmmbuild构建.hmm文件

在有些数据库中是有.hmm文件,只需要下载即可。但是,这仅仅只限于有些大数据库。对于我们自己使用,不可能全部都有,这就需要我们自己构建,**很多教程到这步就是让你收费了…….**。

本教程,讲述其中一种方法吧,希望对大家有所帮助。

hmmbuild构建时,需要使用.sto文件进行构建。因此,我们必须获得.sto文件。

  1. 使用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等的分析,这些分析都需使用目的同源基因的序列。

「如何获得同源基因序列??」


  1. 使用脚本获得
  2. 使用ggffead获得,需要获得同源基因的.gtf文件等信息。
  3. 生信工具获得、如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软件的下载

  1. 网址:https://github.com/CJ-Chen/TBtools
  2. 安装。
  3. 动手运行

2.3.2 提取序列

  1. 准备作物所有的蛋白序列文件(or基因文件)
  2. 目的基因的ID
  3. 打开TBtools,Fasta Extract or Filter (Qyick)
  4. 获得结果

##2.4 目的同源基因motif分析

2.4.1 使用MEME进行motif预测

  1. 网址:https://meme-suite.org/meme/tools/meme
  • 上传相关的fa文件,以及修改相关的参数,进行提交
  • 输出结果输出结果很快,有以下几个结果文件。

2.4.2 motif可视化

对于motifi分析可以参考一下文章:

  1. TBtools | 多图合一至强版教程!进化树 + Motifs + 结构域 + 启动子 + 基因结构 + ….,TBtools开发者本人的教程
  2. TBtools | 基因家族分析 (进化树、Motifs、结构域)
  3. 「或是本篇教程」

MEME网址结果可以给我们的seqlogo信息和motif信息。


  1. 「Seqlogo」结果文件中就有seqlogo文件信息。也可以自己的下载后绘制。按以下操作即可下载序列。也可以下载已有的seqlogo图片。下载后所有的motif序列信息。

  1. 使用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的分析

  1. 下载结果文件MAST XML output,使用TBtools软件进行可视化。

  2. 打开TBtools中的Gene Structure View,只需上传MEME中的XML文件即可,上传上去直接点击Start
    操作:结果:保存!!

注意:我们这里保存的时候最好保存为PDF或SVG格式,输出为「矢量图」

如果我们的教程只是到这里,那么就没有什么意义了。因为,类似非常优秀和详细的教程很多。「绘制出图形是一方面、美化可是重头戏」


在MEME输出文件中,也提供了motif的图形,也可直接使用。

2.5 基因家族保守结构域分析

  1. 使用Batch CD-Search进行预测,网址:https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi
  2. 提交序列信息即可
  3. 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.

可以提供你的邮箱,等运行结束后,直接发送到你的邮箱。如果序列较多,建议提供邮箱。

  1. 下载文件结果文件:
  2. 打开TBtools中的Visualize NCBI CDD DOmainPattern
  3. 输入结果文件和fa文件
  4. 根据自己的需求进行调整即可。
  5. 输出文件。

2.6 进化树分析

进化树分析,在基因家族中是必须的,以及在很多图中都是需要的。进化树分析和绘制,也有很多教程,参考iqtree+ggtree绘制进化树教程、或是你也可以使用MEGA来做分析。

2.6.1 iqtree+ggtree绘制进化树教程

「参考」iqtree+ggtree绘制进化树教程

  1. iqtree获得树文件

「所需软件」


    1. mafft

    1. iqtree「mafft安装」我是使用服务器中运行的,安装可以使用conda
conda install mafft

「iqtree官网」

http://www.iqtree.org/

iqtree功能很强大,大家可以查看软件的官方文档。「安装」

conda install iqtree

软件安装好后直接运行即可。

  1. 序列准备

进化树序列可以使用蛋白序列或核酸序列即可,格式按其准备即可。

>B2LU34
MTSIAFWNAFTVNPFPAAARRSPPPLTPFTSGALSPARKPRILEISHPRTLPSFRVQAIAEDEWESEKKALKGVVGSVAL
AEDETTGADLVVSDLKKKLIDQLFGTDRGLKATSETRAEVNELITQLEAKNPNPAPTEALSLLNGRWILAYTSFAGLFPL
LGAESLQQLLKVDEISQTIDSEGFTVQNSVRFVGPFSSTSVTTNAKFEVRSPKRVQIKFEEGIIGTPQLTDSIVIPDKFE
FFGQNIDLSPFKGVISSLQDTASSVAKTISSQPPIKFPISNSNAQSWLLTTYLDDELRISRADGGSVFVLIKEGSPLLT
>B4F6G1
MTSIAFCNAFTVNPFLAAARRSPPPLTPLTSVALSPARKPRILAIFHPRTFPSFRVQAIAEDEWESEKKTLKGVVGSVAL
AEDEKTGADLVVSDLKKKLIDQLFGTDRGLKATSETRAEVNELITQLEAKNPNPAPTEALSLLNGKWILAYTSFVGLFPL
LGAESLQQLLKVDEISQTIDSEGFTVQNSVRFVGPFSSTSVTTNAKFEVRSPKRVQIKFEEGIIGTPQLTDSIVIPDKVE
FFGQNIDLSPFKGVISSLQDTASSVAKTISSQPPIKFPISNSNAQSWLLTTYLDDELRISRADGGSVFVLILESSPLLT
>O49629
MATVQLSTQFSCQTRVSISPNSKSISKPPFLVPVTSIIHRPMISTGGIAVSPRRVFKVRATDTGEIGSALLAAEEAIEDV
EETERLKRSLVDSLYGTDRGLSASSETRAEIGDLITQLESKNPTPAPTEALFLLNGKWILAYTSFVNLFPLLSRGIVPLI
KVDEISQTIDSDNFTVQNSVRFAGPLGTNSISTNAKFEIRSPKRVQIKFEQGVIGTPQLTDSIEIPEYVEVLGQKIDLNP
IRGLLTSVQDTASSVARTISSQPPLKFSLPADNAQSWLLTTYLDKDIRISRGDGGSVFVLIKEGSPLLNP
  1. mafft比对

使用mafft将序列对齐。

mafft test.fa > test.aligend.fa

我们获得对齐后的数据格式。

  1. 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次

最终,我们可以获得以下结果文件。


  1. 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.62.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绘制。

  1. 点击FileOpen,导入数据
  2. 获得进化树
  3. 调整。全部参数可以在左侧调整即可。包括,大小、间距、距离参数等。以上参数,仅仅只是必要调整的参数,具体看自己的分析进行调整即可,无固定模式。

2.7.2 Figtree绘图的模式

我在前面说过Figtree绘制进化树的图类型很少,只有三种大类型。具体如下所示。

  1. 一般的聚类类型
  2. 圆形circular

2.7.3 Figtree绘制进化树美化图形

如何进行美化,是我们一直在追求的方向。在进化树中分支的上色是必须的,在Figtree中依旧可以做。「注意:我们这里只是简单的说明如何上色,具体操作自己进行。」最终图形可以获得如下图所示。

2.7.4 Figtree导出图形

调整好图形参数,如何导出图形呢?操作如下所示。FileExport JPEG/PNG/PDF.....,导出适合的的图形格式即可,但是建议导出的矢量图。后期AI进行调整。(通过上面导出图形,我们可以看到图形的颜色长度是不同的,这个问题要如何解决,暂时没有找到好的方法。在ggtree绘制中自己也遇到这里的问题。如果在的图形软件中无法解决,只能通过后期解决。)

2.7.4 重新文章中图形

那么如何绘制类似的图形呢?根据前期的参数,只需要进一步优化即可。2. 主图(1) 将图形性状选择圆形(2) 调整Root AngleAngle Rangle调整到适合的形状。2.  分类附图在这个图中,我个人将其进化树分为进化树分类附图。这个图也是使用的Figtree进行绘制。具体操作如下所示。

  • 选择分类图形
  • 调整参数
  • 树枝的宽度可以宽1-2个size
  • 调整自己喜欢的Trabsform Branches
  • 继续调整

「注意:」进化树的分支,主图和附图要一致。为了进一步确定明确两个图的一致性,建议直接在附图中,对分支进行填充颜色。操作与上述一致。

2.7.5 AI合并美化

  1. 打开AI
  2. 新建图形
  3. 导入进化树图形
  4. Ctrl + R打开AI中的标尺、拖出x轴或Y轴参考线
  5. 调整半圆进化树,做到“横平竖直”
  6. Ctrl + A全选,选择图形,Ctrl + C进行复制,或直接进行拖拽到新建图形中。
  7. 调整适合的图形大小,调整时,一直按住shuft,避免图形横纵大小改变。
  8. 建议,在图形中如有新的图形产出,建议每个新的图形都新建立一个图层,利于后期的修改。
  9. 随后就进行进化的调整,我们在这里,需要对AI有一定的基础知识,才可以。比如,如何随意修改图形的形状,类似图例所示。这里操作很繁琐,具体操作自己进行。
  10. 导入进化树分支
  11. 如何线条太细,可以进行调整适合粗细。
  12. 分支添加颜色
  • 新建图形

  • 选择椭圆工具

  • 绘制椭圆,调整适合的分支位置和的添加分支颜色

  • 适当的调整颜色

  • 依次绘绘制即可

  • 字体调整(如果在图形中梯子较小,也可以在AI中调整)使用选择工具,选择调整字体,直接进行修改即可。

  • 调整图形大小

  • 最终出图

  • 也可以直接间监矩形进化树进行进行合并,相比育德圆形或半圆,调整颜色柱就很容易,直接拉成一样长度即可。「细节自己调整。」

2.8 目的基因结构可视化

需要文件:

  1. 目的基因注释文件(GFF or GTF)
  2. 进化树文件(可选)

2.8.1 使用ID和基因组注释文件绘制

  1. 使用TBtools直接操作,依次点击:Gene Structure View结果如图所示:

2.8.2 提取目的基因的注释文件(推荐)

我们会发现,输入ID处也是可以输入进化树文件信息。因此,我们推荐直接提取获得目的基因的注释文件信息,单独使用GTF文件信息或是GFF信息进行绘制。

  1. 获得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";
  1. 进化树获得
    同上的方法获得

  2. MEMExml or MAST.xml文件
    同上

  3. 绘图
    依次提交相关的文件即可

2.9 进化树、Motifs、结构域、基因结构合图绘制

以上的操作,都可以获得单张图形,那么如何多图绘制在一起呢?TBtools也提供了相关的教程,TBtools | 多图合一至强版教程!进化树 + Motifs + 结构域 + 启动子 + 基因结构 + ….,我们可以根据此教程进操作。具体如下:获得结果(来自CJ的教程):


2.10 图形美化

到这里,我们的整张图形就可以获得。但是,只是这样的话,我觉得自己的这个教程就没有意义。我前面说过,我的这个教程重点是图形美化。自己是更喜欢,TBtools单张出图的类型,然后进行AI或PS美化的。软件默认的颜色,我自己不是很喜欢,但是也可以自己调整,也是很方便的哦。

2.10.1 TBtools图形颜色的调整

我们这里只是随意进行调整,图形无任何意义。

  1. 步骤一、点击图形中的方块、右键
  2. 调整色块3、选择先要的色块、点击Selecteed4、更改成功
    但是你会大发现,图中所有一样的颜色色块都会改变。

类似的功能、自己逐渐去摸索。

2.10.2 单张出图

如果上面的方式没有很好实现自己想要的效果。那么,我们就只能单张出图、后面再进行合并。

「注意:在绘图时,我们的要提前想好自己的文章或这张图的颜色设置,以及图形的色调是属于什么类型的。理论上,一整篇文章图形色调和类型要保持一致。」

如果,在后期的调整中。图形颜色需要重新调整,我们可使用AI进行调整或是重新绘制,少量还是比较方便,但是图形又大有多,重画是很奔溃的事情。

三、IA图形美化

美化,我罗列出单个章节进行讲解。表明,是很重要的。以及,图形的美化,需要不断学习和模范大牌期刊的图形类型,以及自己要时刻进行总结和创新。对于创新,这个就比较玄学,每个人的审美不同,逻辑不同,关注点不同…….导致最终看到的点也不同。因此,我们在不是很离谱的创作中,结合自己的审美进行美化即可。「我们要坚信:审美。首先要符合自己,其次,再考虑别人。」只有自己先认同,你才有可能让其他人也认同!

3.1 使用工具

1.「推荐使用的工具」:AI、PS

如果不知道类似软件的,自己百度。

  1. 如何安装
  • 有钱人:购买正版
  • 穷人(和我一样):薅羊毛,使用破解版
  1. 如何获取安装包

在本公众号中回复关键词获得。

  • PS安装包关键词:PS
  • AI安装包关键词:AI

或是你自己寻找相关版本的安装包即可。

「提示:」请自己输入正确的关键词(每次看到有些同学们的关键词,真的很无语……)


3.2 实际操作

  1. 打开AI,新建图层A4
  2. 导入进化树,适当调整进化树的宽度和字体大小
  3. 依次导入的目的基因的motif、基因结构域等图形。并依次按进化树基因名进行排序即可。
  4. 为后期的图形的整齐性,我们使用参考线进行对齐,便于后期的调整。「注意:这里看到我们的motif的图形颜色很难看,这就是前期没有考虑颜色的结果。因此,我一直强调,文章图形颜色统一的重要性,图形颜色搭配合理,你的论文已经成功1/3了。」换一种颜色就感觉好多了呀。
  5. 添加基因结构图添加图形的操作都是一样的,不做多赘述。
  6. 如何美化
    对于美化,每个人的要求不一致,只要符合你的审美即可。我们在这里就直接添加渐变色。
  7. 新建一个图层
    新建图层置于最底层。
  8. 选择图形工具
  9. 利用进化树的分支,将其进行分类
  10. 填充颜色(根据自己的喜好)
  11. 更改透明图
  12. 渐变色
  • 不透明度:60
  • 中间位置:10-50%
    结合实际情况调节。
  1. 最后图形图形很多细节需要自己耐心调节,这里只是做示范,相对比较粗糙。

四、多物种共线分析

共线分析依旧是使用TBtools,哈哈哈哈,做基因家族TBtools可以帮你完成80%的生信分析。毫不夸张!!!!!
TBtools共线分析的教程很多,我们以零基础多物种间共线性分析教程作为参考(也不是作为参考了,是直接按他的步骤进行操作)。其他参考教程:全基因组共线性分析无限个!物种共线性分析结果可视化、任何人!一键完成物种间的共线性分析与可视化。

4.1 需要文件

  1. 参考基因组fa文件
  2. 注释文件GFF or GTF

TBtools可以对无限个作物进行共线分析,牛!!!

4.2 染色体统一命名

在这个教程中,有这样的一个步骤,如果你需要,你就进行操作。

  1. gtf文件进行ID prefix

3. fa文件进行ID prefix

4.3 实操

  1. 打开one step MCScanx小程序
  2. 输入两个作物的文件信息
  3. 点击开始Start
  4. 如果是多个作物,那么依次进行两两比较。比如:共线结果是以这样的顺序:Tomato-LA-Arabidopsis

比对顺序:

  • Tomato-LA
  • LA-Arabidopsis
  1. 比对结果GFF文件合并
  • 打开Text Merge for MCScanX程序合并多个的MCScanX的结果文件中的GFF文件拖拽文件

6.比对结果ChrLayout.tab.xls文件合并7. 比对结果geneLinks.tab.xls文件合并同上操作!8. 合并文件
最终获得以下3个文件,用于绘制图形。9. 要在共线中显色的基因ID

Solyc03g062790.3.1
Solyc10g018590.2.1
Solyc01g104320.4.1
Solyc03g083420.4.1
AT4G22240.1
AT2G35490.1
AT1G51110.1
AT5G53450.3
........
  1. 绘图。打开Multiple synteny plot输入参数输出图形

注意,在输出图形中,我们可以看到作物染色体位置是有改变的。那么,如何更改呢?回答:直接更改Chr文件即可。更改这里的顺序即可!

五、同源目标基因元件预测

目标基因的元件预测,我们这里主要介绍使用两个网站进行。

5.1 提取目标基因上游2000bp

参考教程顺式作用元件预测和新的可视化方式植物启动子-顺式作用元件-批量提取-预测-可视化分析,同样是使用TBtools操作。

  1. 需要文件
  • 作物参考基因组fa文件
  • 注释文件GFF or GTF
  • 目标基因ID
  1. 直接使用TBtools中的Gtf /Gff3 Sequences Extractor获得每个基因的fa序列输出文件点击Initalize,选择CDS选择上游2000bp的fa序列

  2. 目标基因的fa序列,打开Fasta Extract or Filter (Quick)输出结果文件:

  3. 查看信息是否正确,打开Fasta Stats

  4. 转换序列(全部为大写),打开Sequence Manipulate (Rev&Comp

5.2 提交预测网址进行顺式作用预测

预测,这里使用两个网站进行预测,分别是PlantCarePLCAE

5.2.1 使用Plantcare进行预测

网址:http://bioinformatics.psb.ugent.be/webtools/plantcare/html/

  1. 上传序列后,Plant可以提供你自己的邮箱,运行结束后,结果直接发送到你的邮箱中。
  2. 邮箱中获得结果,根据你的序列多少,10分钟以上吧!
  3. 结果
  4. 使用execl打开后
1. 基因ID;
2. 顺式作用元件名称;
3. 顺式作用元件序列;
4. 顺式作用元件的起始位置;
5. 顺式作用元件的长度;
6. 顺式作用元件所在的链的方向;
7. 物种名;
8. 顺式作用元件所在的功能分类;

删除某些不需要的结果:「需要删除:」

1. 剔除第2列为空的行
2. 剔除第2列为unnamed的行
3. 最后一列,无功能作用的

具体删除的数据,根据自己的分析来做。最后,可以删除掉1000行以内

— 

来自顺式作用元件预测和新的可视化方式,这个意见有重要的参考意义。如果不合并,导致元件的作用太多,绘制出的图形颜色太杂,且不好看。5. 绘图绘图前还需要准备基因的长度文件输入数据,设置参数结果:在TBtools中也可以输入进化树文件。


我们这里也可以使用的起那么AI中的呢进化树进行模板进行美化。

5.2.2 PLACE进行预测

网址:https://www.dna.affrc.go.jp/PLACE/?action=newplace

  • 缺点:PLACE一次最大只能输入20条基因序列,有一定的限制性。获得结果为网页版,如要整理,只能手动整理或使用脚本进行整理。
  • 优点:速度快!
  1. 获得结果每个基因为单独的,需要自己整理。
  • 只给元件名称、开始位置、序列、功能(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的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!