美文网首页
[R数据]提取序列,绘制fasta进化树

[R数据]提取序列,绘制fasta进化树

作者: 花生学生信 | 来源:发表于2024-09-09 14:17 被阅读0次
提取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)

相关文章

网友评论

      本文标题:[R数据]提取序列,绘制fasta进化树

      本文链接:https://www.haomeiwen.com/subject/vezeljtx.html