泛基因家族分析R语言计算核苷酸多样性

泛基因家族分析通常会分析 kaks 和 核苷酸多样性，之前的推文介绍过使用wgd这个软件计算dnds，
泛基因组分析计算核心可变基因家族的kaks
这个软件调用的是codeml。计算dnds的时候可以选择把比对好的文件保存下来，直接用这个保存好的比对文件计算核苷酸多样性。codeml计算kaks的代码
wgd ksd –n_threads 8 02.mcl 02.cds/ZH13.v2.CDS.fasta –preserve
02.mcl 是某个基因家族所有序列id 每行是一个基因家族，序列id用制表符分割 ZH13.v2.CDS.fasta 是所有的cds序列
保存的比对文件

image.png
这个比对格式忘记叫什么格式了
R语言里 pegas 包里有一个函数 nuc.div()函数可以直接计算核苷酸多样性
这个函数的输入是 ape包 read.dna()这个函数读取fasta格式的比对序列，读取以后是个DNAbin格式
如何用R语言的函数把wgd 保存的比对文件读入并保存成 DNAbin这个格式 没有找到现成的函数
ape 这个R包里有函数 as.DNAbin() 可以把矩阵转换成 DNAbin的格式
比如
library(tidyverse)
c(“A”,”T”,”C”,”G”,”A”,”A”,”C”,”G”) %>%
matrix(nrow = 2,byrow = TRUE) %>%
ape::as.DNAbin()
library(tidyverse)
c(“A”,”T”,”C”,”G”,”A”,”A”,”C”,”G”) %>%
matrix(nrow = 2,byrow = TRUE) %>%
ape::as.DNAbin() %>%
pegas::nuc.div()
把上面提到的格式数据读取进来转换成一个矩阵就行了
library(tidyverse)
read_lines(“D:/Jupyter/panGeneFamily/GF_000001.fasta.msa.nuc”,
skip = 1) %>%
matrix(ncol=2,byrow=TRUE) %>%
as.data.frame() %>%
pull(V2) %>%
str_split(pattern = “”) %>%
unlist() %>%
matrix(nrow = 5,byrow = TRUE) %>%
as.DNAbin()

欢迎大家关注我的公众号