之前看基因家族相关的文章绘制 sequences logo 基本会使用 meme+webologo 的组合,前一阵子在公众号 R语言中文社区 推送的文章中发现了R语言的 ggseqlogo 可以用来绘制 sequences logo,好多种颜色可供选择,结果非常漂亮,帮助文档十分详细,使用也相对简单。
简单记录使用meme鉴定保守基序,简单处理meme结果文件成R语言可读取的格式,最后用ggseqlogo绘图的过程
先上结果图
本地meme所用到的的命令
/usr/yanming/Bioinformatic/MEME/meme_4.12.0/bin/meme PgLHC_candidate_16.fasta -protein -oc PgLHC_motif
-nostatus -mod anr -nmotifs 5 -minw 6 -maxw 50 -maxsize 600000000
-oc 后接输出的结果文件
-mod 后接的应该是方法,有三个参数可选 oops zoops anr 具体是什么意思自己也不太懂,好像是anr参数比较好,但是相对运行时间会比较长
-nmotif 后接要搜索的motif的数量
-minw -maxw 后接motif的最短和最长的长度
-maxsize 如果蛋白的fasta文件中序列数量较多就需要相应增到后面接的数字
一直都会遇到报错,因为可以输出结果,也就没有仔细探究报错的原因
这是结果文件中的保守基序之一,将其复制到一个新的文件中简单处理成R可以读入的格式,我用的是linux的cut命令
cut -d ")" -f 2 | cut -d " " -f 2 > ggseq_logo_input.txt
首先用括号作为分隔符提取第二列,用竖线将上一步的结果传递给下一步作为输入,再用空格作为分隔符提取第二列,再用大于号重定向到结果文件中
meme在线用法
结果
第一步
PS:今天Xftp链接服务器莫名其妙的就没有了文件的下载权限,什么原因???搞不太明白,正常拖拽上传,却不能拖拽下载