美文网首页
[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