论文
https://www.nature.com/articles/s41559-022-01771-6#code-availability
论文没有权限下载
但是查看数据代码链接的时候发现github主页上提供了论文的下载链接
论文中的图做的都非常好看，而且提供数据和代码，我们可以找来学习
数据代码链接
https://github.com/sebepedroslab/chromatin-evolution-analysis
今天的推文我们学习一下论文中提供的画多序列比对的代码。我没有在论文中找到对应的图，只是github的链接里有数据和代码
论文中国提供的代码
library(stringr) library(seqinr) library(msa) library(Biostrings) list_ali = list( "data/20220620/archaeal_Ntails/tails_archaea-sub-MXKK.g.fasta", "data/20220620/archaeal_Ntails/tails_archaea-sub-mcl1.g.fasta", "data/20220620/archaeal_Ntails/tails_archaea-sub-mcl2.g.fasta", "data/20220620/archaeal_Ntails/tails_archaea-sub-mcl4.g.fasta" ) for (i in 1:length(list_ali)) { fai = list_ali[[i]] ali = Biostrings::readAAMultipleAlignment(fai, format = "fasta") # plot alignment msa::msaPrettyPrint( ali, askForOverwrite=FALSE, shadingMode = "functional", shadingModeArg = "chemical", showConsensus = "none", showLogo = "none", shadingColors = "blues", logoColors = "chemical", psFonts = TRUE, paperWidth = 32, paperHeight = 10, output = "pdf", file = sprintf("%s.colour.pdf",basename(fai))) } 这里作图用到的是msa这个R包
读取多序列比对的数据使用的是Biostrings这个R包
上面的代码是写了一个简单的循环，做了四个数据的图，我试着做其中一个图，但是遇到了报错
Error in texi2dvi(texfile, quiet = !verbose, pdf = identical(output, "pdf"), : unable to run pdflatex on 'tails_archaea-sub-MXKK.g.fasta.colour.tex' LaTeX errors: ! Paragraph ended before \inf@@get was complete. <to be read again> \par l.26 ...cal/Temp/RtmpagEAWD/seq397c6ecb2093.fasta} ! Misplaced alignment tab character &. \msfline ->\par & & & & @ l.26 ...cal/Temp/RtmpagEAWD/seq397c6ecb2093.fasta} ! Misplaced alignment tab character &. \msfline ->\par & & & & @ l.26 ...cal/Temp/RtmpagEAWD/seq397c6ecb2093.fasta} ! Misplaced alignment tab character &. \msfline ->\par & & & & @ l.26 ...cal/Temp/RtmpagEAWD/seq397c6ecb2093.fasta} 搞不懂是什么原因
R语言里做多序列比对的图还有更好的选择使用 ggmsa 这个R包
安装
devtools::install_github("YuLab-SMU/ggmsa") 这里R语言必须是4.1以上的
输入数据是比对好的fasta文件的路径
作图代码
library(ggmsa) fai<-"data/20220620/archaeal_Ntails/tails_archaea-sub-MXKK.g.fasta" pdf(file = "aligned_fasta.pdf",width = 9.4,height = 4) ggmsa(fai)+ geom_seqlogo() + geom_msaBar() dev.off() 
image.png
这个我认为比论文中提供的代码出图要好看的多
这个是论文中的代码出的图
示例数据和代码可以在公众号后台回复20220620获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！