MCScanX使用小记

简介: 一篇不知道拖了多久的笔记……欠下的债早晚还是要还的……0.软件下载及安装官网:http://chibba.pgml.uga.edu/mcscan2/参考教程:其实之前已经有两位前辈写过教程了的,我这只是再给需要的朋友举个例子打个样。

一篇不知道拖了多久的笔记……欠下的债早晚还是要还的……

0.软件下载及安装

官网:

http://chibba.pgml.uga.edu/mcscan2/

参考教程:
其实之前已经有两位前辈写过教程了的,我这只是再给需要的朋友举个例子打个样。详情请戳下面两个链接

http://blog.biochen.com/archives/878

https://www.jianshu.com/p/740cb9eccf2b

下载软件:

wget http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip

解压:

unzip MCScanX.zip

cd到MCScanX后,按照

http://www.genek.tv/article/45

介绍的方法,给

  • msa.h

  • dissect_multiple_alignment.h

  • detect_collinear_tandem_arrays.h

这三个文件前面添加上#include <unistd.h>

即,用vim将文件打开后在开头加上上述的参数段。这个错误的原因是,MCScanX 不支持64位系统。
然后make一下就安装好了

1.文件准备

从TAIR上下载TAIR10的pep和gff3数据,要先用notepad++的正则替换功能(vim的批量替换如果会用的话也是可以的)将蛋白质的基因ID后面的注释信息给删掉,否则会影响后续的工作。接下来处理gff3文件。
先看一眼原始的数据长啥样:


img_74abd6cae77d940b99a93cdb3de5aacc.png
TAIR上下载的拟南芥的gff3数据

我们需要的是第1、4、5列和第9列的ID部分
用awk、grep和vim的全局替换功能可以得到所需要的简化版gff3文件:

 cat TAIR10.gff3 |awk '/gene/ {print $1  "\t" $3 "\t"  $4 "\t" $5 "\t" $9 }'> tmp.gff3
$3是用来标记用的,我们只取含有gene的列
img_58839230b2665a2ae1af0d332c5e0384.png
得到的结果

但是会有一些问题,比如:


img_a5ca79accc7b8c9944f1a436ca13943c.png
搜索gene时的副产物

在搜索gene的时候只要包含gene的词就会被检索到,于是有以下三个混进来的:

  • pseudogene

  • transposable_element_gene

  • mRNA_TE_gene
    于是我就:

grep -v 'mRNA_TE_gene' tmp2.gff3 > tmp3.gff3
grep -v 'pseudogene' tmp3.gff3 > tmp4.gff3 
grep -v 'transposable_element_gene' tmp4.gff3 > tmp5.gff3

由于当时学艺不精比较菜就只能先暂时这样用一下了……


img_8637deffd29b0482fd4f33da4377af46.jpe
就是!又不是不能用

现在想想好像在上一步awk的时候可以:

 cat TAIR10.gff3 |awk '/^gene/ {print $2 "\t" $0}' > tmp.gff3

加上一个“^”表示以gene开头,应该就可以免去上面那繁琐的几步了,最后生成的时候再倒一下回来就好了。
可以用

grep -v "^gene" tmp3.gff3

来检测是否还包含有非gene的行。
经过一系列的awk操作之后获得最终的结果:


img_c13af3b7dea853f0e3f07c392c5eb80e.png
AT.gff3

另一个物种如法炮制即可。
updates:也可以这么写,三步搞定:

awk '$3 == "gene" {print $1 "\t" $3 "\t" $4 "\t" $5 "\t" $9}' TAIR10.gff3 > new.gff3
用vim全局替换处理ID列的冗余信息
用awk将信息位置调换正确

2.建库

makeblastdb -in AT.pep -dbtype prot -parse_seqids -out ATdb

因为我需要做的是maca对拟南芥的共线性比较,所以只需要对拟南芥建库即可。(如果需要物种内也进行共线性比较的话这里需要把两个物种的蛋白质文件cat到一起来建库)

3.比对

nohup blastp -query maca.pep -db ATdb -out AT_maca.blast -evalue 1e-10 -num_threads 30 -outfmt 6 -num_alignments 5 &

参数稍微解释一下。

-query 被比对的物种的蛋白质文件
-db 上一步建的库
-out 输出结果
-evalue 设置输出结果的e-value值,值越小相似度越高。
-num_threads 比对所用的线程数
-outfmt 输出的格式
-num_alignments 输出比对上的序列的最大值条目数

这一步还是要点时间的。记得挂nohup。

4.运行MCScanX

将之前的两个物种的gff3文件cat成一个文件后命名为AT_maca.gff
将上一步生成的AT_maca.blast和AT_maca.gff放到新建的AT_maca文件夹内

MCScanX ./AT_maca

之后就会生成一个叫AT_maca.collinearity的文件和一个AT_maca.html文件夹。


img_0f70136a2e003e716e6e704a374de701.png
AT_maca.collinearity

AT_maca.html里面是这样滴:


img_373217ed3ccb544ab7e582a1a47ea567.png
AT和maca的对应关系

5.下游分析

java circle_plotter -g AT_maca.gff -s AT_maca.collinearity -c circle.ctl -o AT_maca_circle.png

这一步需要修改修改circle.ctl文件


img_e04740d57c2b9296a3b4ab5e9e094e0e.png
ctl文件

只要把你希望画出来的序列的ID给写进去就好了
可以画dot图,bar图和circle图,个人觉得circle图最好看了~


img_c0e55ead20bc1f24522fceb7e97b7441.png
结果展示

大概就是这样了。如果有什么错误的地方欢迎讨论指出哟。
相关文章
|
数据采集 芯片
GWAS全基因组关联分析入门教程
GWAS全基因组关联分析入门教程
|
C++
如何使用MACS进行peak calling
MACS2是peak calling最常用的工具。 callpeak用法 这是MACS2的主要功能,因为MACS2的目的就是找peak,其他功能都是可有可无,唯独callpeak不可取代。
3874 0
|
Linux Python
竟然还有这么好用的工具,xshell你可以滚蛋了~~
竟然还有这么好用的工具,xshell你可以滚蛋了~~
251 0
|
机器学习/深度学习 缓存 网络安全
服务器中的conda环境
服务器中的conda环境
884 1
|
6月前
|
存储 数据可视化 数据挖掘
R语言绘制圈图、环形热图可视化基因组实战:展示基因数据比较
R语言绘制圈图、环形热图可视化基因组实战:展示基因数据比较
|
6月前
hexo 博客插入本地图片时遇到的坑
在使用Hexo博客时,作者遇到了将本地图片插入文章的问题。安装了`hexo-asset-image`插件(v1.0.0)并启用资源文件夹功能后,图片并未正确显示。问题在于图片路径未被转换。通过卸载`hexo-asset-image`,改用`hexo-asset-img`插件,并调整图片引用格式为`![alt](文章名/图片名.png)`,成功解决了图片显示问题。此问题源于`hexo-asset-image`在Hexo 5及更高版本中已失效。
|
6月前
|
Ubuntu Shell
ubuntu卸载conda
ubuntu卸载conda
422 1
|
6月前
基因组组装:Hifiasm 使用教程
基因组组装:Hifiasm 使用教程
225 1
|
数据采集 设计模式 存储
全基因组重测序流程【超细致!!】
全基因组重测序流程【超细致!!】
|
数据可视化 数据挖掘 Go
RNA-seq丨转录组分析标准流程与常用工具
RNA-seq丨转录组分析标准流程与常用工具