说在前面:
官网http://chibba.pgml.uga.edu/mcscan2/index.php信息很全面,有示例代码/Manual手册,你值得拥有。以下内容基本出自官网教程,现在你可以自行打开官网学习了~
安装
只能在Mac OS (via X11) or Linux systems上执行。Linux systems,需要JDK
和libpng
(一般都有)
官网下载压缩包
Then simply put MCscanX.zip into a directory and run:
unzip MCscanX.zip
cd MCScanX
make
运行示例文件
分析拟南芥与葡萄的共线性
输入文件需要有:at_vv.gff ; at_vv.blast 【在 data
文件夹中】
In MCScanX package directory, type
./MCScanX data/at_vv
./ 表示当前目录;data是输入文件所在目录,注意gff和blast文件名保持一致
3秒搞定
Pairwise collinear blocks written to data/at_vv.collinearity [17.348 seconds elapsed] Tandem pairs written to data/at_vv.tandem
上述两个文件at_vv.collinearity
at_vv.tandem
就是输出文件,默认输出到data文件夹
Writing multiple syntenic blocks to HTML files
除此之外,还输出个 at_vv.html文件夹,里面就是屏幕上有输出的html文件。
一切顺利,安装成功
准备输入文件
先看看拟南芥-葡萄的示例文件格式
at_vv.gff
提供基因的位置信息,两个4列的 gff文件 合并而来
#head
at3 AT3G19630 6818676 6820674
at5 AT5G11220 3577057 3577854
at2 AT2G29110 12506880 12510552
#tail
vv1 GSVIVT01012140001 1098524 1102524
vv18 GSVIVT01034900001 16234068 16271170
vv8 GSVIVT01025554001 14085153 14089377
at_vv.blast
提供比对信息,有AT-AT;AT-VV;VV-AT;VV-VV。出现4种比对的原因是把at和vv的蛋白文件合并成了一个,然后既做database又做query,进行blastp比对。
注意:blast文件中的基因名要与gff文件中的基因名对应!!gff中有且仅有一个。
基因id提前检查好,不然就全盘皆输从头再来……
xyz.blast文件
(个人需求,一般不用这步)批量修改基因id,区分2个物种
sed -i 's/evm.model.Contig/HHGevm.model.Contig/' HHG.pep
sed -i 's/evm.model.Contig/BHGevm.model.Contig/' BHG.pep
将两个物种的蛋白文件合并
cat HHG.pep.fa BHG.pep.fa >> twoHG.fa
将该蛋白文件既做database,又做query,进行blastp
blastp的m8输出格式,也就是blast+版本的 -outfmt 6
blast版本
formatdb -i database_file -p T
#(F表示核酸库 T表示蛋白库)
blastall -i query_file -d database –p blastp –e 1e-5 –b 5 –v 5 –m 8 –o xyz.blast
解释2个参数——
-b : 显示的比对结果的最大数目,缺省值250
-v : 单行描述(one-line description)的最大数目,缺省值500
blast+版本
makeblastdb -in database.fa -dbtype prot -out db_name -parse_seqids
-dbtype 后接序列类型,nucl为核酸,prot为蛋白
-parse_seqids 推荐加上
-out 后接数据库名
-logfile 日志文件,如果没有默认输出到屏幕 (-logfile xxx.log)
blastp -query query.fa -db db_name -out xyz.blast -evalue 1e-5 -outfmt 6 -num_alignments 5 -num_threads 16
-num_alignments 显示比对数Default = 250 (相当于blast版本的 -b 5)
-num_descriptions:单行描述的最大数目 default=50 -num_threads:线程
得到 HHG_BHG.blast
xyz.gff文件
要求:
sp# gene starting_position ending_position
sp is the two-letter short name for the species; # is the chromosome number. (For example, the second chromosome of Arabidopsis thaliana should be denoted as at2.)
第一列表示物种和染色体号。如拟南芥2号染色体 at2 . (不是2个字母也ok啦)
第二三四列分别是 基因名 起始位点 终止位点
gff格式介绍:
https://www.jianshu.com/p/b26c285bd027
使用awk
修改,根据自己的文件以及自己的需求
awk 'BEGIN{OFS="\t"} $3 == "gene" {split($9,x,";");split(x[1],y,"=");id=y[2];gsub("\"","",id);print "HHG"$1,id,$4,$5}' HHG.gff > HHG1.gff
awk 'BEGIN{OFS="\t"} $3 == "gene" {split($9,x,";");split(x[1],y,"=");id=y[2];gsub("\"","",id);print "BHG"$1,id,$4,$5}' BHG.gff > BHG1.gff
#生成2份整理过的GFF,然后修改id,再合并【个人需求,略】
'>' 直接把内容生成到指定文件,会覆盖源文件中的内容
'>>' 尾部追加,不会覆盖掉文件中原有的内容。
参考:
https://www.plob.org/article/11373.html
https://www.jianshu.com/p/8373e50722f6
MCScanX分析共线性
注:如果是多个物种,将所有需要的blastp文件合并为一个xyz.blast,物种的gff文件合并成一个 xyz.gff
将准备好的 blast gff文件放入data文件夹 使用 mv
移动
进入MCScanx文件夹
接下来就跟示例一样操作了。