零基础-6小时-完全重现某个基因家族分析文章（的分析部分）
先说明
本推文出发点就是……个人觉得有趣
本推文已尽可能地保证零基础的朋友能在windows下完全重现，但不保证每个人都能重现。毕竟我没有义务。所以最好不要针对本文的步骤向我提问（星球的朋友除外），近期赶课题，木有时间。
（另，个人时间和精力有限，大群可以学习交流，但我不一定会回复（也没有义务），直接私信我讨论交流的朋友，请先微信转账或附图支付宝转账-）。
出发点
前几天某公众号放出基因家族分析服务，有朋友告知我，那个收费是一个家族三万RMB。对于这个家族，事实上，我个人觉得挺值的。只是，该公众号拿了别人家的基因家族文章（我与文章作者已沟通过啦，作者表示躺着中枪），容易让别人误会（难道是该文作者出来搞？）。可以在该广告文看到TBtools一个比较简单的输出图片。我个人自然是内心不舒服。既然如此，那我就写一篇推文，告诉没做过基因家族分析的朋友，
没有任何基础，照样可以在短时间内（一天之内），完成一个基因家族的分析
如此这般，你可以省下三万RMB。
Anyway，感觉大家都爱做基因家族。既然如此，我就直接重现那个用了TBtools输出图片却没有引用TBtools的文章。
重现过程
整个文章的分析比较简单，我们尽可能用TBtools来实现其中的各个步骤，能简化的就简化，力求结果一致
序列下载
在弹出的窗口中，选中需要的文件，然后点击Download Selected Files
在随后弹出的下载窗口，保存下载文件
下载完毕，解压压缩文件，并进入annotation，即可看到这个对应的注释信息
扫描对应基因家族的模式下载
不管，那就做做看先
看到文章中，做的是一个糖转运基因家族-sugar transporter (STP) gene family，好嘛，这个其实应该找下拟南芥的数据来看看，就知道有什么保守domain了。也可以输入到pfam
使用HMMER构建数据库，输入对应的命令
hmmpress Sugar_tr.hmm
使用TBtools，只提取我们关注的pfam模式
使用Text Block Extractor And Filter工具
使用文本编辑器（如notepad++）打开这个文件，发现扫描出来的序列非常多，且evalue也很低，（同样的情况，同样的操作用在拟南芥也是一样，而使用pfam数据库这个模式对应的5000+个植物序列进行blastp，也是一样的结果）。如下图
这就很麻烦，我本来以为这是一个很好搞的家族，但是这明显是一个坑，而我也没看到这个家族多少报道，然而~，既然要搞他，那就要搞完他。既然如此，我们能走的一条路就是一样，建树拆分支，用什么软件?
当然是我们的好基友MEGA，
MST家族建树-拆STP亚家族
建完之后，可以直接看到只有一个分支是有STP的，无论如何，虽然和论文写得是不同的，但是我们还是拆出来了
于是我们得到这个列表
Manes.03G025000.1.p_pacid_32363352_transcript_Manes.03G025000.1_locus_Manes.03G025000_ID_Manes.03G025000.1.v6.1_annot-version_v6.1
Manes.16G111500.1.p_pacid_32343266_transcript_Manes.16G111500.1_locus_Manes.16G111500_ID_Manes.16G111500.1.v6.1_annot-version_v6.1
AT5G23270.1 | Symbols:STP11ATSTP11 | SUGAR TRANSPORTER 11sugar transporter 11 | Chr5:7839132-7840874 FORWARD LENGTH 514
Manes.03G180400.1.p_pacid_32363574_transcript_Manes.03G180400.1_locus_Manes.03G180400_ID_Manes.03G180400.1.v6.1_annot-version_v6.1
Manes.04G053100.1.p_pacid_32328102_transcript_Manes.04G053100.1_locus_Manes.04G053100_ID_Manes.04G053100.1.v6.1_annot-version_v6.1
Manes.14G139800.1.p_pacid_32361595_transcript_Manes.14G139800.1_locus_Manes.14G139800_ID_Manes.14G139800.1.v6.1_annot-version_v6.1
AT3G05960.1 | Symbols:ATSTP6STP6 | SUGAR TRANSPORTER 6sugar transporter 6 | Chr3:1783587-1785334 REVERSE LENGTH 507
Manes.01G067100.1.p_pacid_32357903_transcript_Manes.01G067100.1_locus_Manes.01G067100_ID_Manes.01G067100.1.v6.1_annot-version_v6.1
Manes.15G027100.1.p_pacid_32350816_transcript_Manes.15G027100.1_locus_Manes.15G027100_ID_Manes.15G027100.1.v6.1_annot-version_v6.1
Manes.02G122100.1.p_pacid_32332522_transcript_Manes.02G122100.1_locus_Manes.02G122100_ID_Manes.02G122100.1.v6.1_annot-version_v6.1
AT1G77210.1 | Symbols:STP14AtSTP14 | sugar transport protein 14 | Chr1:29009036-29010980 REVERSE LENGTH 504
AT3G19930.1 | Symbols:ATSTP4STP4 | SUGAR TRANSPORTER 4sugar transporter 4 | Chr3:6935048-6936841 FORWARD LENGTH 514
AT1G07340.1 | Symbols:ATSTP2STP2 | sugar transporter 2 | Chr1:2254873-2256712 FORWARD LENGTH 498
Manes.11G102900.1.p_pacid_32355740_transcript_Manes.11G102900.1_locus_Manes.11G102900_ID_Manes.11G102900.1.v6.1_annot-version_v6.1
Manes.17G033300.1.p_pacid_32365841_transcript_Manes.17G033300.1_locus_Manes.17G033300_ID_Manes.17G033300.1.v6.1_annot-version_v6.1
AT1G11260.1 | Symbols:STP1ATSTP1 | sugar transporter 1SUGAR TRANSPORTER 1 | Chr1:3777460-3780133 FORWARD LENGTH 522
Manes.18G067900.1.p_pacid_32350241_transcript_Manes.18G067900.1_locus_Manes.18G067900_ID_Manes.18G067900.1.v6.1_annot-version_v6.1
Manes.06G132700.1.p_pacid_32347460_transcript_Manes.06G132700.1_locus_Manes.06G132700_ID_Manes.06G132700.1.v6.1_annot-version_v6.1
Manes.03G180600.1.p_pacid_32364574_transcript_Manes.03G180600.1_locus_Manes.03G180600_ID_Manes.03G180600.1.v6.1_annot-version_v6.1
Manes.01G164600.1.p_pacid_32359497_transcript_Manes.01G164600.1_locus_Manes.01G164600_ID_Manes.01G164600.1.v6.1_annot-version_v6.1
Manes.01G182200.1.p_pacid_32359132_transcript_Manes.01G182200.1_locus_Manes.01G182200_ID_Manes.01G182200.1.v6.1_annot-version_v6.1
AT1G50310.1 | Symbols:STP9ATSTP9 | sugar transporter 9SUGAR TRANSPORTER 9 | Chr1:18635984-18638110 FORWARD LENGTH 517
AT5G26340.1 | Symbols:MSS1ATSTP13STP13 | SUGAR TRANSPORT PROTEIN 13 | Chr5:9243851-9246994 REVERSE LENGTH 526
Manes.15G030700.1.p_pacid_32351597_transcript_Manes.15G030700.1_locus_Manes.15G030700_ID_Manes.15G030700.1.v6.1_annot-version_v6.1
Manes.16G020100.1.p_pacid_32343239_transcript_Manes.16G020100.1_locus_Manes.16G020100_ID_Manes.16G020100.1.v6.1_annot-version_v6.1
Manes.15G027300.1.p_pacid_32351659_transcript_Manes.15G027300.1_locus_Manes.15G027300_ID_Manes.15G027300.1.v6.1_annot-version_v6.1
Manes.02G122000.1.p_pacid_32332038_transcript_Manes.02G122000.1_locus_Manes.02G122000_ID_Manes.02G122000.1.v6.1_annot-version_v6.1
Manes.06G132600.1.p_pacid_32347668_transcript_Manes.06G132600.1_locus_Manes.06G132600_ID_Manes.06G132600.1.v6.1_annot-version_v6.1
这说明什么？说明，即使作者很不小心，并没有在文章中写清楚；使用了TBtools出图却没有告诉我们他用了TBtools的情况下，我们还是能重复出他们的结果，这个研究，做得挺好的哈，起码我重复出来了，只是以一种很诡异的方式
绘制一些图片
恩，这个图片没什么信息含量，我们需要加上结构域！！！
将序列直接提交到NCBI CDD数据库，然后下载对应的信息
使用TBtools输出一张图片，挺丑的。稍微整理一下对应的信息，那就把结构域展示在基因结构上！！
随便配的颜色，感觉比较一般，先这样吧
做一张与作者论文中几乎一样的motif图片
太简单了，因为这是重点，作者就用的TBtools，那么我们做一张几乎一样的，如下（注意配色！配色！配色！是一套系统出来，TBtools的配色）
作者估计控制了motif个数，anyway，我们设置了20个，所以出来的motif数目多了一些，没关系。继续做剩下的部分
整理出一个motif序列信息！！！
既然我们又有MEME的结果，自然就会有一个MAST的结果，整理一下就可以得到这个报个啊，怎么整理？TBtools啊
物种间基因家族分类
作者拿了木薯的20个序列，拟南芥的14个STP以及蓖麻的STP一起，建了一棵树…依旧是套路操作。建树的操作，太过于简单，不需要重复了
家族内部基因的重复事件分析
这个有点麻烦，一般这种分析有两个策略：
基于家族成员之间的相似性，覆盖率，以及在染色体上的距离，确定出串联重复的家族成员（粗暴！！）
这类操作一般是，提取所有蛋白序列为一个文件，文件比对到自身，设定一致性和覆盖率的指标，确定出哪些基因是高度相似的，随后基于这些基因的相对位置做主观判断，结束。基于共线性分析（直接使用MCScanX软件，需要在linux下运行）
物种进化过程中，大部分邻近的基因在短时间内高度连锁，亦即在不同物种中，我们可以发现染色体之间存在一些大片段，这些片段上的基因相对位置也相对保守，亦即所谓的共线性区块。基于这一特点，以共线性区块作为参考（家族成员侧翼的基因为参考），我们不仅可以快速确定出成员间的串联重复时间，也可以确定出片段重复或全基因组复制事件，进而更好地说明家族成员之间的进化关系。
鉴于以上，这里我们直接使用我认为目前在比较基因组中使用最为广泛的一个软件MCScanX。这个软件需要在Linux下面运行，而其实他几乎不耗费计算资源，所以我们直接在windows下面开一个虚拟机（参考博文, https://www.jianshu.com/p/6caf60f58869）
首先，我们要做一个Blast比对，直接在使用TBtools，这个过程，普通电脑2个线程的话，可能两个小时左右，等呗。
于是我们会得到这个文件
将其移动到VirtualBox的共享目录下面，然后我们进入VitrualBox，进入共享目录所在路径（请参考千上文提到的windows下如何安装VitrualBox虚拟机，并提高性能和共享文件夹）
进入之后可以看到对应的文件
在Ubuntu虚拟机中安装MCscanX
打开终端，输入以下命令
首先是下载MCScanX安装包
wget http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip
随后解压与编译，
unzip MCscanXcd MCscanX make
调整一下源码，重新编译
(简书有Bug，一些字符组合可能输入不了，贴图)
perl -i -pe 见上图 msa.h dissect_multiple_alignment.h detect_collinear_tandem_arrays.hmake
实在不行，我们直接重启虚拟机，似乎还是解决不了，那么我们换几个命令
for file in msa.h dissect_multiple_alignment.h detect_collinear_tandem_arrays.h;do perl -pe 见上图 $file > $file.bak mv $file.bak $filedonemake
OK，如果只是出现以下问题，那么我们就基本搞定了
这个报错比较好解决，环境中找不到javac，事实上，MCScanX的出图还是相对粗糙的，用java码写的几个出图脚本，怎么说呢，我们还是要用一下
那就安装java，注意需要安装的是 jdk
sudo apt-get install default-jdk
OK，重新编译一次MCScanX的绘图脚本
javac downstream_analyses/*.java
搞定了，那么就开始进行共线性分析，因为我们是一次性工作，直接在程序目录下搞，搞完删掉这个程序目录好，整个分析其实只需要5分钟
cp ../Mesculenta_305_v6.1.blast.table.txt Me.blast perl -F'\t' -lane 'next unless $F[2] eq "mRNA";print join qq{ \t},$F[0],$F[-1]=~s/ID=([^;]+).v6.1.*$/$1/r.".p",$F[3],$F[4]' ../Mesculenta_305_v6.1.gene_exons.gff3 > Me.gff
查看重复事件
./duplicate_gene_classifier Me
搞定了，前后用了我3个多小时，完成这个教程（虚拟机的另外花了可能有接近1个小时）
回到windows下面，我们就可以使用TBtools，看看我们鉴定的20个STP都是有什么基因复制关系
然后就可以直接批量提取了
打开输出的文件，然后对应一下ID，就搞定了，基本可以确定下来哪些成员是串联重复，哪些成员是片段复制或者是全基因组复制的结果
其实，还有另外几项关于基因家族的分析，但是我不想写了….
给一张图
./MCScanX Mecd downstream_analyses (echo 800;cut -f1 ../Me.gff |grep 'Chr'|u niq|paste -sd,) > Me.STP.family.ctl cut -f1 ../STP.gene.type.txt |paste -sd'\t' > Me.STP.ids java family_circle_plotter -g ../Me.gff -s ../Me.collinearity -c Me.STP.family.ctl -f Me.STP.ids -o STP.circle.PNG
其实重新用circos画，还是很方便的。懒得讲了
就到这里吧
这个事情也告一段落了
再声明
本推文出发点就是……个人觉得有趣
参与讨论交流