基因家族分析

基因家族分析是一种常见的生物信息学分析套路，也是生物信息学数据挖掘发表小文章常用的分析方法，和GEO挖掘等类似。基因家族的分析鉴定可以用pfam上的hmm文件进行基因家族的检索鉴定，也可以用blast的方法进行比对鉴定，通常是适用拟南芥对应的基因家族进行比对鉴定。
植物转录因子数据库，点击访问。

基因家族分析思路及文章撰写思路

基因家族分析是继GEO数据挖掘后，一种新的生物信息学挖掘策略。

如何做基因家族研究，可以参考这个帖子：http://www.planttech.com.cn/blog/58882464a46。

数据准备

需要准备的数据主要是参考基因组数据，包括fasta格式的序列文件、gff或gtf格式的基因组注释文件、蛋白质序列文件（通常是每个转录本的蛋白序列）、cDNA序列等文件。如果有转录组数据的话进行对应的转录组分析即可。除开上述这些文件外，还需要适用的文件还有.hmm格式的文件。

软件准备

只需要会使用Linux系统，会安装Docker即可，然后下载组学大讲堂的镜像即可。点击浏览镜像地址。Docker的安装适用方法参考\@ref(WSL4Docker)。

分析过程

  • mRNA_ID与基因ID的提取
    由于一个基因会对应多个转录本，因此一个基因下会对应多个mRNA的编号。在后续的分析中，每个基因只需要选择一个转录本的编号进行分析即可，因为每个基因不同的转录本的序列差异是较小的。提取的代码：
1
2
perl code/script/mRNAid_to_geneid.pl data/unzip_data/*.gff* results/step.1.get.mRNA.and.gene.ID/mRNA2geneID.txt
perl code/script/geneid_to_mRNAid.pl data/unzip_data/*.gff* results/step.1.get.mRNA.and.gene.ID/geneID2mRNAid.txt
  • 检索结构域
    这一步主要是以.hmm文件为基础检索该物种蛋白序列中含有该结构域的序列。输入文件包括.hmm文件和蛋白文件，输出hmmsearch的检索结果。其中用于后续筛选的是evalue这个参数，部分文章以0.001为阈值。of那一列表示的是某个基因对应的这个结构域有几个。
1
hmmsearch --domtblout results/step.2.domain.search/hmm.txt --cut_tc data/unzip_data/*.hmm data/unzip_data/*.pep*
  • 选择结构域
    由于一个基因的单个转录本可能会比对到多个结构域，因此需要对比对到的结构域进行选择。默认选择的是第一个结构域。下面代码的最后一个参数是hmmsearch输出文件里面的E-value,如果需要全部的第一个结构域，将阈值设置为1即可。
1
perl code/script/domain_xulie.pl results/step.2.domain.search/hmm.txt data/unzip_data/*.pep* results/step.2.domain.search/domain.fa 1.2e-28
  • 多序列比对
    之所以要进行多序列比对，是因为下载的.hmm文件是来自很多物种的这个结构域组成的隐马尔科夫模型，进行多序列比对后将该物种检索到的结构域序列进行比对，再次构建该物种该基因家族的隐马尔科夫模型，会更加准确。
1
echo -e '1\nresults/step.2.domain.search/domain.fa\n2\n1\nresults/step.2.domain.search/out.aln\nr.domain.search/out.dnd\nX\n\n\nX\n' |clustalw
  • 重构隐马尔科夫模型
1
hmmbuild results/step.2.domain.search/new.hmm results/step.2.domain.search/out.aln
  • 重新进行检索
    利用构建得到的新的隐马尔科夫模型重新进行检索结构域。
1
hmmsearch --domtblout results/step.2.domain.search/new.out.txt --cut_tc results/step.2.domain.search/new.hmm data/unzip_data/*.pep*
  • 筛选输出结果
    对重新检索后的结果进行筛选，也是对E-value进行筛选。
1
grep -v "^#" results/step.2.domain.search/new.out.txt|awk '$7<0.001 {print}' > results/step.2.domain.search/domain.new.out.selected.txt
  • 去除重复的ID
    上一步筛选得到的是该种中哪些基因是潜在的目标基因家族成员，而一个基因对应了多个mRNA，因此，只需要在筛选后的每个基因中选择一个具有代表性的mRNA进行后续的分析即可。这个提取唯一ID的步骤需要手动完成（PS：手动完成也很快）。手动挑选完mRNA的ID放在第一列，另存为文件uniqueID.txt。
1
perl code/script/select_redundant_mRNA.pl results/step.1.get.mRNA.and.gene.ID/mRNA2geneID.txt results/step.2.domain.search/domain.new.out.selected.txt results/step.2.domain.search/remove_redundant_IDlist.txt
  • 提取蛋白序列
    在得到基因ID后需要提取蛋白序列进行后续的分析。在SMART或者Pfam或NCBI CDD确认这些基因是真真正正含有该结构域，没有的基因要剔除！在SMART中没有检索到结构域的基因在gene.not.in.SMART.txt中；在Pfam中全都是WRKY结构域，对应文件为Pfam.results.txt。
1
2
perl code/script/get_fa_by_id.pl results/step.2.domain.search/uniqueID.txt data/unzip_data/*.pep* results/step
.2.domain.search/pep.need.confirm.fa
  • 计算蛋白质分子量等
1
perl code/script/stat_protein_fa.pl results/step.2.domain.search/pep.need.confirm.fa results/step.2.domain.search/pep.MW.txt
  • 构建进化树
    选择利用软件CLUSTALW进行多序列比对，然后利用MEGA构建进化树。CLUSTALW输出结果转换成.fasta格式的方法参考\@ref(pac4xiang)。
  • Motif分析
1
meme results/step.3.seq.and.tree/pep_confirmed.fa -protein -oc results/step.4.motif/ -nostatus -time 18000 -maxsize 6000000 -mod anr -nmotifs 10 -minw 6 -maxw 100
  • 基因结构分析
1
perl code/script/get_gene_exon_from_gff.pl -in1 results/step.2.domain.search/uniqueID.txt -in2 data/unzip_data/*.gff* -out results/step.5.gene.structure/gene_exon_info.gff
  • 基因在染色体上的定位
1
2
3
samtools faidx data/unzip_data/*.dna*
cp data/unzip_data/*.fai results/step.5.gene.structure/
perl code/script/get_gene_weizhi.pl -in1 results/step.2.domain.search/uniqueID.txt -in2 data/unzip_data/*.gff* -out results/step.5.gene.structure/mrna_location.txt
  • 顺式作用元件分析
    脚本默认的启动子长度是1500bp。将提取得到的启动子序列上传到Plane CARE进行分析。
1
perl code/script/get_promoter.pl data/unzip_data/*dna.top* results/step.5.gene.structure/mrna_location.txt results/step.6.cis.acting.element/promoter.txt
  • 基因家族成员的亚细胞定位分析
    两个网站：WolfPsort和Cello。

基因组共线性分析

  • 从geneid2mRNAid.txt文件中每个基因挑一个转录本ID存储到allmRNAID.txt，也就是每个基因挑选一个转录本。

  • 创建文件夹：

1
2
3
4
mkdir MCScanX
cd MCScanX

mkdir mcscan
  • 获取转录本的位置信息：
1
perl ../script/get_mRNA_position.pl allmRNAID.txt ../Arabidopsis_thaliana.TAIR10.41.gff3 AT.gff
1
2
3
4
5
6
7
8
9
10
11
# 分别是染色体位置，转录本ID，起始位置和终止位置
1       Pno01G013682.t2 7993891 8000167
1       Pno01G011207.t1 6521081 6522022
1       Pno01G002381.t1 15982661        15984634
1       Pno01G006728.t1 32747373        32748586
1       Pno01G001298.t1 134750756       134753812
1       Pno01G000565.t1 11549686        11558266
1       Pno01G005598.t1 23981500        23984173
1       Pno01G012140.t1 707333  718266
1       Pno01G006114.t1 27812391        27816942
1       Pno01G001551.t1 14123991        14125021
  • 获取挑选的转录本的蛋白序列：
1
perl ../script/get_fa_by_id.pl allmRNAID.txt ../Arabidopsis_thaliana.TAIR10.pep.all.fa pep.fa
  • 对所有蛋白进行建库：
1
makeblastdb -in pep.fa  -dbtype prot -title pep.fa
  • blast比对：
1
2
blastall -i pep.fa -d pep.fa -e 1e-10  -p blastp  -b 5 -v 5 -m 8 -o mcscan/AT.blast
cp AT.gff mcscan/AT.gff
  • MCScanX查找共线性：
1
MCScanX mcscan/AT

得到的结果是整个基因组的结果：

1
2
3
4
5
-rwxrwxrwx 1 xiang xiang  17M Dec  6 23:59 sanqi.blast
-rwxrwxrwx 1 xiang xiang 339K Dec  7 10:48 sanqi.collinearity # 基因组的所有大片段复制
-rwxrwxrwx 1 xiang xiang 1.7M Dec  7 10:47 sanqi.gff
drwxrwxrwx 1 xiang xiang 4.0K Dec  7 10:48 sanqi.html
-rwxrwxrwx 1 xiang xiang  39K Dec  7 10:48 sanqi.tandem # 基因组的所有串联重复

  • 结果解读：

    • 串联重复：当一行里面的两个基因都是感兴趣的基因家族基因的时候，就是需要的串联重复基因。
    1
    2
    3
    4
    5
    Pno01G001040.t1,Pno01G001042.t1
    Pno01G001410.t1,Pno01G001412.t1
    Pno01G001524.t1,Pno01G001525.t1
    Pno01G001558.t1,Pno01G001560.t1
    Pno01G001628.t1,Pno01G001629.t1
    • 大片段复制：当一行里面的两个基因都是感兴趣的基因家族基因的时候，就是需要的大片段复制基因。
    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
    ############### Parameters ###############
    # MATCH_SCORE: 50
    # MATCH_SIZE: 5
    # GAP_PENALTY: -1
    # OVERLAP_WINDOW: 5
    # E_VALUE: 1e-05
    # MAX GAPS: 25
    ############### Statistics ###############
    # Number of collinear genes: 10371, Percentage: 21.66
    # Number of all genes: 47870
    ##########################################
    ## Alignment 0: score=269.0 e_value=4.6e-10 N=6 1&1 plus
      0-  0:        Pno01G005483.t1 Pno01G005033.t2   8e-54
      0-  1:        Pno01G005544.t1 Pno01G005040.t3  7e-159
      0-  2:        Pno01G005828.t1 Pno01G005055.t1  2e-103
      0-  3:        Pno01G005829.t1 Pno01G005058.t1       0
      0-  4:        Pno01G005897.t1 Pno01G005088.t1  1e-164
      0-  5:        Pno01G006072.t2 Pno01G005092.t1   1e-80
    ## Alignment 1: score=362.0 e_value=2.8e-14 N=8 1&1 minus
      1-  0:        Pno01G000527.t2 Pno01G004765.t1   2e-75
      1-  1:        Pno01G000535.t1 Pno01G004763.t1       0
      1-  2:        Pno01G000542.t2 Pno01G004760.t2  4e-152
      1-  3:        Pno01G000545.t2 Pno01G004757.t1   7e-68
      1-  4:        Pno01G000558.t1 Pno01G004741.t6   8e-20
      1-  5:        Pno01G000568.t2 Pno01G004737.t2  4e-144
      1-  6:        Pno01G000599.t1 Pno01G004735.t1   2e-68
      1-  7:        Pno01G000636.t1 Pno01G004716.t1  1e-144
    ## Alignment 2: score=290.0 e_value=3.9e-11 N=7 1&1 minus
      2-  0:        Pno01G006165.t1 Pno01G002002.t1   1e-23
      2-  1:        Pno01G006238.t1 Pno01G001895.t1       0
      2-  2:        Pno01G006400.t1 Pno01G001832.t2   3e-87
      2-  3:        Pno01G006416.t2 Pno01G001818.t1   3e-48
      2-  4:        Pno01G006555.t1 Pno01G001725.t1  2e-130
      2-  5:        Pno01G006745.t1 Pno01G001684.t1   2e-46
      2-  6:        Pno01G006864.t1 Pno01G001662.t1       0

  • 找打目标基因家族中存在串联重复的基因：

1
2
3
perl /biosoft/MCScanX/MCScanX/downstream_analyses/detect_collinearity_within_gene_families.pl -i WRKY_family.txt -d mcscan/AT.collinearity -o WRKY.collinear.pairs

perl ../script/get_tandem_gene.pl -id gene_list.txt -tandem mcscan/AT.tandem -od ./ -name WRKY
  • 基因加倍绘图：
1
2
3
4
5
6
7
8
9
10
11
mkdir circos
cd circos

cp ../MCScanX/mcscan/AT.collinearity .   
cp ../MCScanX/WRKY.collinear.pairs .

perl ../script/colline_v3.pl -gff ../MCScanX/AT.gff -list WRKY.collinear.pairs -colline AT.collinearity -od data -name WRKY

cd data

circos -conf config3.txt -outputdir ./ -outputfile WRKY

配置文件：

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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
chromosomes_units=1000000    #刻度单位Mb
karyotype=./chr.info  #染色体信息配置文件
show_tick_labels=yes
show_ticks=yes
spacing=10u
<ticks>  #设置染色体刻度
    color=black
    format=%d
    multiplier=1e-6
    radius=1r
    thickness=2p
    <tick>
        size=10p
        spacing=5u
    </tick>
    <tick>
        color=black
        format=%d
        label_offset=10p
        label_size=25p
        show_label=yes
        size=15p
        spacing=25u
        thickness=4p
    </tick>
</ticks>

<ideogram> #染色体绘制设置
    fill=yes   #是否填充颜色
    label_font=default
    label_parallel=yes
    label_radius=dims(image,radius)-60p
    label_size=45
    radius=0.8r  #设置半径，以免基因名称过长超出显示范围
    show_label=yes
    <spacing>
        default=0.005r
    </spacing>
    stroke_color=dgrey
    stroke_thickness=2p
    thickness=0.03r
</ideogram>

<links>  #设置连线
    bezier_radius=0r
    bezier_radius_purity=0.75
    color=black
    crest=0.5
    <link>  #基因家族共线性
        bezier_radius=0r
        bezier_radius_purity=0.75
        color=set2-8-qual-1
        crest=0.5
        file=./nlr.link.txt    #输入文件
        radius=0.98r
        <rules>
            <rule>
                color=red
                condition=var(intrachr)
            </rule>
            <rule>
                color=red
                condition=var(interchr)
            </rule>
        </rules>
        thickness=8
        z=20   #层数
    </link>
    <link>  #全基因组共线性
        bezier_radius=0r
        bezier_radius_purity=0.75
        color=230,230,230,0.2
        crest=0.5
        ribbon=yes
        file=./genome.blocklink.txt   #输入文件
        radius=0.98r
        thickness=1
        z=15
    </link>
    radius=0.40r
    thickness=1
</links>
<plots>
    <plot>
       color=black
	   label_snuggle=yes  #如多个文本文字距离过近，避免重叠  参考：https://www.omicsclass.com/article/678
        file=./nlr.text.txt  #输入文件
        label_font=condensed
        label_size=24p
        link_color=red
		link_dims=10p,10p,50p,20p,10p
        link_thickness=2p
        r0=1r
        r1=1r+500p
        rpadding=0p
        padding=0p
        show_links=yes
        type=text
    </plot>
    type=histogram
</plots>


<colors>
<<include etc/colors.conf>>
<<include etc/brewer.conf>>
#<<include etc/colors_fonts_patterns.conf>>
#<<include colors.ucsc.conf>>
#<<include colors.hsv.conf>>
</colors>

<fonts>
<<include etc/fonts.conf>>
</fonts>

<image>
<<include etc/image.conf>>
</image>
<<include etc/housekeeping.conf>>

软件JCVI：https://github.com/lixiang117423/jcvi

💌lixiang117423@foxmail.com
💌lixiang117423@gmail.com

生物信息学
#生物信息学 #分子生物学 #基因家族
基因家族分析
https://lixiang117423.github.io/article/genefamily/
作者
小蓝哥
发布于
2021年11月25日
许可协议