0.参考
https://github.com/zhaotao1987/SynNet-Pipeline
Github地址， 部分文本我是直接翻译的
https://www.pnas.org/cgi/doi/10.1073/pnas.1801757116
18年发表在PNAS上; 文章的引言提到其最大的优势在于处理大型复杂的数据集
以我的测试来看，构建10个物种的进化树只需要1-2小时，并且占用计算资源很少
类似地还有wgdi，wgdi文章引言提到wgdi的这种基于共线性构建系统发育
可以避免不完全谱系分选和基因流对系统发育带来的影响
wgdi有缘再说, (https://wgdi.readthedocs.io/en/latest/)
安装DIAMOND和MCScanX, 略
DIAMOND: https://github.com/bbuchfink/diamond
MCScanX: http://chibba.pgml.uga.edu/mcscan2/
git clone https://github.com/zhaotao1987/SynNet-Pipeline.git
2.数据准备
数据是物种的gff文件、cds/pep文件皆可
当下载基因组时，我们将看到基因名称以各种不同的方式命名。比如:
Aradu.20JM2 genotype-assembly-annot=V14167.a1.M1
DTZ79_01g11390
AT1G01010.1 | NAC domain containing protein 1 | Chr1:3760-5630 FORWARD LENGTH=429 | 201606
Bv1_000040_cpku.t1 cDNAEvidence=88.9
除了第二个基因名(可能还需要一个有意义的物种前缀)，其他基因名都需要缩短。例如，删除所有空格后面的字符，也更好地取代名字中的点，并为物种添加一个唯一的3-5个字母长的前缀。比如>ath_AT1G01010， >Aradu_20JM2， ..看起来不错还要注意其他几件事
。检查fasta中的序列总数是否与GFF文件中的记录相匹配(以第3列的“基因”计数)
。确保fasta文件中单个编码区域没有替代转录本或蛋白质序列
。删除每个序列末尾的'*'(这可能会在Diamond稍后构建序列数据库时出现问题)
。gff文件中的第一列也就是染色体名，也需要处理，不同物种的染色体名要能相互区分将gff文件转为以下格式(4列,第一列染色体名, 第二列基因名,第三列/第四列是起始/终止bp):
athChr1 ath_AT1G01010 3631 5899
athChr1 ath_AT1G01020 5928 8737
athChr1 ath_AT1G01030 11649 13714最后，我们需要提取最长转录本
以下参考https://blog.csdn.net/u012110870/article/details/113544271
此外，同时，对于每个基因我们只需要一个转录本，我通常使用最长的转录本作为该基因的代表。之前我写过一个脚本 (get_the_longest_transcripts.py) 提取每个基因的最长转录本，有一个新的脚本用来根据参考基因组和注释的GFF文件生成wgdi的两个输入文件，脚本地址为: https://github.com/xuzhougeng/myscripts/blob/master/comparative/generate_conf.py
我的注释文件格式, 注意到红色箭头的地方都是分号
3 获得输入文件
./iTools_Code/iTools Fatools getCdsPep -Ref spp1.fa -Gff spp1.gff -OutPut spp1
gunzip -c spp1.pep.fa.gz > spp1.pep
#提取pep，获得第一个输入文件
awk '{if($3=="mRNA")print}' spp1.gff | awk -F ";" '{print $1}' | sed 's/ID=//g' | awk 'BEGIN{ FS="\t";OFS="\t" }{print $1,$9,$4,$5}' > spp1.bed
#获得第二个输入文件
#每个基因组都得到了一个.cds/.pep文件和一个.bed文件
#spp1.bed, spp2.bed, spp3.bed, spp4.bed ......
4 构建共线性网络
首先把SynetBuild-X.sh的第70行, 改成我们自己的物种缩写
原来是 array=(Umar Vpac Mmus)，改成 array=(spp1 spp2 spp3 spp4)
更改如图所示的地方
./SynetBuild-X.sh 6 5 25 10
#参数照抄的, 感兴趣参阅github和原文章，有对不同参数组合进行测序
the four variables are: tophits, -s, -m, -threads
tophits (k or b): number of target sequences for each query sequence to keep from Diamond blast.
-s: number of anchors to call a synteny block in MCScanX, default is 5, okay for angiosperms.
-m: number of genes upstream and downstream to search for anchors, default is 25. Lower this value for a stricter detection.
-threads: this is used for both Diamond and MCScanX. No need to set this super high, 10-30 is already quite good, depending on your CPU.
会生成一个文件夹，文件夹中有最终的输出文件 “SynNet-k5s5m25”
其表头大概是这种格式:
Alyr0-0b 936.0 AlyrAL1G19310 AlyrAL1G65100
Alyr0-1b 936.0 AlyrAL1G19350 AlyrAL1G65150
Alyr0-2b 936.0 AlyrAL1G19400 AlyrAL1G65290
Alyr0-3b 936.0 AlyrAL1G19420 AlyrAL1G65430
Alyr0-4b 936.0 AlyrAL1G19430 AlyrAL1G65470
Till now, synteny network construction was done, which belong to the upstream part of the analysis. Downstream analysis and applications can be highly dynamic. For example you can..
- Filter the edgelist of your interested genes/ gene family, and analyze the sub-networks.
- Cluster the entire network, and analyze the patterns of the resulting clusters.
- Build phylogeny based on synteny clusters.
- 构建系统发育
#提取最终的输出文件 “SynNet-k6s5m25”的最后两列
awk -F "\t" '{print $3"\t"$4}' SynNet-k6s5m25 > SynNet.2cols
Rscript infomap.r SynNet.2cols SynNet.infoclusters
Rscript Phylogenomic_Profiling.r SynNet.infoclusters profiled profiled_clustered
#这个脚本的22行需要修改，修改为自己的物种顺序
在如图所示地方修改, 和前文的修改相同
会得到第一个输出的图, 这个地方很像泛基因组
最后
Rscript BinaryDataandTranspose.r SynNet.infoclusters.profiled.clustered 1
会输出一个后缀为 _cluster_binary_transposed的文件
在该文件的第一行添加物种的数量和cluster的数量
eg: https://github.com/zhaotao1987/SynNet-Pipeline/blob/master/sample-profiled_size_1_19751_clusters_binary_transposed
也就是这一步会弹出的2个数字
改完长这一
vim SynNet.infoclusters.profiled.clustered_size_1_24950_clusters_binary_transposed
iqtree -s SynNet.infoclusters.profiled.clustered_size_1_24950_clusters_binary_transposed
-bb 1000 -alrt 1000 -nt 10 -m MK+R+FO -redo -st MORPH
就拿到树文件了
6.共线性网络的其他分析 - 待续
网友评论