提取16列文件做进化树
# 读取文件
data <- read.csv("1vir.csv", header = TRUE)
# 检查数据列数是否足够
if (ncol(data) < 16) {
stop("数据列数少于 16 列")
}
# 将数据转换为数据框（确保是 data.frame 类型）
data <- as.data.frame(data)
# 创建一个新的数据框只包含第一列和第十六列，同时去除包含缺失值的行
new_data <- data.frame(
col1 = data[, 1][complete.cases(data[, c(1, 16)])],
col16 = data[, 16][complete.cases(data[, c(1, 16)])]
)
# 根据第一列合并第十六列（使用 aggregate 函数）
merged_data <- aggregate(col16 ~ col1, data = new_data, FUN = function(x) paste(x, collapse = ","))
# 重命名列
colnames(merged_data) <- c("第一列", "合并后的第十六列")
# 保存结果为新的 CSV 文件
write.csv(merged_data, "merged_data.csv", row.names = FALSE)
# 读取 CSV 文件
data <- read.csv("merged_data.csv", header = TRUE)
# 提取第一列和第二列
col1 <- data[, 1]
col2 <- data[, 2]
# 构建 FASTA 格式的字符串
fasta_str <- ""
for (i in 1:length(col1)) {
fasta_str <- paste0(fasta_str, paste0(">", col1[i], "\n", col2[i], "\n"))
}
# 保存为 FASTA 文件
write(fasta_str, file = "output.fasta")
###BiocManager::install("msa")
require(msa)
mySequenceFile <- readAAStringSet("output.fasta")
#多序列比对
myFirstAlignment <- msa(mySequenceFile)
head(mySequenceFile)
library(ggplot2)
require(seqinr)
myAlignment <- msaConvert(myFirstAlignment, type="seqinr::alignment")
d <- dist.alignment(myAlignment, "identity")
#构建NJ树
require(ape)
tree <- nj(d)
#画树并输出到PDF文件ggtree.pdf
require(ggtree)
环形进化树
##环状图
p1<-ggtree(tree, layout='circular', ladderize=FALSE, size=0.8, branch.length="none",col="red")+
geom_tiplab2(hjust=-0.3)+
geom_tippoint(size=1.5,col="blue")+
geom_nodepoint(color="black", alpha=1/4, size=2) +
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))
p1
# 图例位置、文字大小
###长方形图
p2<- ggtree(tree, layout='rectangular', size=0.8, col="deepskyblue3") +
geom_tiplab(size=3, color="purple4", hjust=-0.05)+
geom_tippoint(size=1.5, color="deepskyblue3")+
geom_nodepoint(color="pink", alpha=1/4, size=5)+
theme_tree2()
ggsave(p1, file="tree_circular.pdf", width=9, height=9)
ggsave(p2, file="tree_rectangular.pdf", width=9, height=9)
