解析具有两列详细说明行的父子关系的文件答案

【问题标题】：Parsing a file that has two columns detailing a parent-child relationship of rows解析具有两列详细说明行的父子关系的文件
【发布时间】：2015-06-17 18:22:09
【问题描述】：

我花了很多时间环顾四周，但找不到解决我的具体问题的方法。非常感谢任何帮助。

我有一个包含 5 列的 data.frame，其中每一行详细说明了一个特定的分类级别（即超级王国细菌）。这些级别中的每一个都通过父子关系与另一行（和级别）相关。 taxon_id 是一个唯一值（即超级界细菌的taxon_id 为 2），并通过 parent_id 变量链接到其他行（即变形杆菌门嵌套在超级界细菌中）。快速回顾一下，taxon_id 是唯一值（和父值），而 parent_id 是每个子项与父项的关联方式。还有一个丰度变量，将所有子行的所有丰度相加，使其成为累积总和）。

以下是数据框的示例：

abundance <- c(0.8157, 0.8153, 0.4947, 0.4807, 0.3444, 0.3444, 0.3444, 0.3444, 0.3444, 0.2104, 0.1466, 0.1401, 0.1299, 0.1299, 0.123, 0.1112, 0.1034, 0.1034, 0.1034, 0.1034, 0.07469, 0.07469, 0.06833, 0.06833, 0.06384, 0.06384, 0.06314, 0.06314, 0.06298, 0.06271, 0.0625, 0.0625, 0.0625, 0.06184, 0.06184, 0.05218, 0.05218, 0.05218, 0.04547, 0.04547, 0.04547, 0.04532, 0.03121, 0.02122, 0.01438, 0.01438, 0.01438, 0.01438, 0.01438, 0.01411, 0.01329, 0.01329, 0.01329, 0.01329, 0.0121, 0.0121, 0.0121, 0.01116, 0.01071, 0.008966, 0.008088, 0.008067, 0.007638, 0.007638, 0.00761, 0.00761, 0.00761, 0.00761, 0.00761, 0.007381, 0.006763, 0.00676, 0.00676, 0.00676, 0.006413, 0.004669, 0.004669, 0.004669, 0.004669, 0.004669, 0.004637, 0.004637, 0.004637, 0.004637, 0.004637, 0.003418, 0.003346, 0.002479, 0.002479, 0.002479, 0.002479, 0.001818, 0.001818, 0.001818, 0.001818)
taxon_id <- c(131567, 2, 1224, 1236, 135619, 224372, 1177179, 519051, 59753, 1239, 91061, 186826, 1301, 1300, 1313, 869308, 2037, 85003, 1760, 201174, 543, 91347, 570, 573, 186802, 186801, 31979, 1485, 1491, 1221327, 85006, 85023, 33882, 1160710, 36807, 72407, 1185419, 1460422, 662, 641, 135623, 672, 1219061, 85007, 135622, 53246, 267888, 1307437, 176102, 1229202, 85025, 1817, 1206725, 37327, 1743, 31957, 85009, 1747, 1159092, 28211, 1284824, 1236101, 1716, 1653, 1883, 2062, 85011, 1160718, 114687, 81852, 38284, 1350, 1352, 525279, 1311, 198251, 54526, 82117, 198252, 314261, 90371, 59201, 28901, 590, 1218143, 525260, 862512, 204441, 41295, 191, 95605, 356, 41294, 1073, 95607)
parent_id <- c(1, 131567, 2, 1224, 1236, 135619, 519051, 59753, 224372, 2, 1239, 91061, 1300, 186826, 1301, 1313, 85003, 1760, 201174, 2, 91347, 1236, 543, 570, 186801, 1239, 186802, 31979, 1485, 1491, 2037, 85006, 85023, 36807, 33882, 573, 1460422, 72407, 641, 135623, 1236, 662, 672, 2037, 1236, 267888, 135622, 176102, 53246, 672, 85007, 85025, 37327, 1817, 31957, 85009, 2037, 1743, 1313, 1224, 573, 573, 1653, 85007, 2062, 85011, 2037, 114687, 1883, 186826, 1716, 81852, 1350, 1352, 1301, 54526, 82117, 28211, 198251, 198252, 59201, 28901, 590, 543, 90371, 38284, 38284, 28211, 204441, 41295, 191, 28211, 356, 41294, 1073)
rank <- c("no_rank", "superkingdom", "phylum", "class", "order", "family", "strain", "species", "genus", "phylum", "class", "order", "genus", "family", "species", "strain", "order", "subclass", "class", "phylum", "family", "order", "genus", "species", "order", "class", "family", "genus", "species", "strain", "suborder", "family", "genus", "strain", "species", "subspecies", "no rank", "strain", "genus", "family", "order", "species", "strain", "suborder", "order", "genus", "family", "strain", "species", "strain", "family", "genus", "strain", "species", "genus", "family", "suborder", "species", "strain", "class", "strain", "strain", "genus", "family", "genus", "family", "suborder", "strain", "species", "family", "species", "genus", "species", "strain", "species", "genus", "no rank", "no rank", "species", "strain", "strain", "subspecies", "species", "genus", "no rank", "strain", "strain", "order", "family", "genus", "species", "order", "family", "genus", "species") 
taxonomy_name <- c("cellular_organisms", "Bacteria", "Proteobacteria", "Gammaproteobacteria", "Oceanospirillales", "Alcanivoracaceae", "Alcanivorax_hongdengensis_A-11-3", "Alcanivorax_hongdengensis", "Alcanivorax", "Firmicutes", "Bacilli", "Lactobacillales", "Streptococcus", "Streptococcaceae", "Streptococcus_pneumoniae", "Streptococcus_pneumoniae_SPN021198", "Actinomycetales", "Actinobacteridae", "Actinobacteria", "Actinobacteria", "Enterobacteriaceae", "Enterobacteriales", "Klebsiella", "Klebsiella_pneumoniae", "Clostridiales", "Clostridia", "Clostridiaceae", "Clostridium", "Clostridium_botulinum", "Clostridium_botulinum_CDC66177", "Micrococcineae", "Microbacteriaceae", "Microbacterium", "Microbacterium_laevaniformans_OR221", "Microbacterium_laevaniformans", "Klebsiella_pneumoniae_subsp._pneumoniae", "Klebsiella_pneumoniae_subsp._pneumoniae_ST258-K26BO", "Klebsiella_pneumoniae_subsp._pneumoniae_ST258", "Vibrio", "Vibrionaceae", "Vibrionales", "Vibrio_vulnificus", "Vibrio_vulnificus_NBRC_15645_=_ATCC_27562", "Corynebacterineae", "Alteromonadales", "Pseudoalteromonas", "Pseudoalteromonadaceae", "Pseudoalteromonas_agarivorans_S816", "Pseudoalteromonas_agarivorans", "Vibrio_vulnificus_B2", "Nocardiaceae", "Nocardia", "Nocardia_brevicatena_NBRC_12119", "Nocardia_brevicatena", "Propionibacterium", "Propionibacteriaceae", "Propionibacterineae", "Propionibacterium_acnes", "Streptococcus_pneumoniae_PNI0010", "Alphaproteobacteria", "Klebsiella_pneumoniae_VAKPC252", "Klebsiella_pneumoniae_JHCK1", "Corynebacterium", "Corynebacteriaceae", "Streptomyces", "Streptomycetaceae", "Streptomycineae", "Streptomyces_auratus_AGR0001", "Streptomyces_auratus", "Enterococcaceae", "Corynebacterium_accolens", "Enterococcus", "Enterococcus_faecium", "Enterococcus_faecium_TX1330", "Streptococcus_agalactiae", "Candidatus_Pelagibacter", "SAR11_cluster", "unclassified_Alphaproteobacteria", "Candidatus_Pelagibacter_ubique", "Candidatus_Pelagibacter_ubique_HTCC1002", "Salmonella_enterica_subsp._enterica_serovar_Typhimurium", "Salmonella_enterica_subsp._enterica", "Salmonella_enterica", "Salmonella", "Salmonella_enterica_subsp._enterica_serovar_Typhimurium_str._STm1", "Corynebacterium_accolens_ATCC_49725", "Corynebacterium_accolens_ATCC_49726", "Rhodospirillales", "Rhodospirillaceae", "Azospirillum", "Azospirillum_sp._B4", "Rhizobiales", "Bradyrhizobiaceae", "Rhodopseudomonas", "Rhodopseudomonas_sp._B29")      
mydata <- data.frame(abundance, taxon_id, parent_id, rank, taxonomy_name)

我想将此 data.frame 解析为更易于分析的不同格式。这看起来像这样：

Taxa <- c("cellular_organisms(no_rank)", "cellular_organisms(no_rank)_Bacteria(superkingdom)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)_Oceanospirillales(order)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)_Oceanospirillales(order)_Alcanivoracaceae(family)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)_Oceanospirillales(order)_Alcanivoracaceae(family)_Alcanivorax(genus)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)_Oceanospirillales(order)_Alcanivoracaceae(family)_Alcanivorax(genus)_Alcanivorax_hongdengensis(species)", "cellular_organisms(no_rank)_Bacteria(superkingdom)_Proteobacteria(phylum)_Gammaproteobacteria(class)_Oceanospirillales(order)_Alcanivoracaceae(family)_Alcanivorax(genus)_Alcanivorax_hongdengensis(species)_Alcanivorax_hongdengensis_A-11-3(strain)")
Proportion <- c(0, 0.0103, 0.01396, 0.00176, 0, 0, 0, 0, 0.3444)
OutcomeData <- data.frame(Taxa, Proportion)

从原始数据可以看出，Bacteria 的丰度为 0.8153，但减去所有子行的丰度后，我们得到的余数为 0.0103。我想为所有完整的父子关系生成一个这样的表。关于这个示例表的附加快速说明，为了简单起见，我排除了一些子样本，因此这里的数学不会加回之前的父样本，数据分支很多，但我想要一个更清晰的示例。

【问题讨论】：

比例应该是节点的丰度减去所有孩子的丰度吗？所以父母和祖父母的财富会减去孩子的财富……
@cr1msonB1ade 是的。父母和祖父母的财富等会减去孩子的财富。一旦减去所有子丰度，然后将余数用作该级别丰度的丰度。

标签： r

【解决方案1】：

这是我草拟的一个快速函数 - 它写得不是很好，但看起来很有效。在相信它之前你应该多测试一下：

mydata$rank<-as.character(mydata$rank)
deparsefile<-function(df){
  df$newlabel<-paste(df$rank,df$taxonomy_name)
  df$remainder<-df$abundance
  parentstodo=c(131567)
  while(length(parentstodo)>=1){
    parentstodo1<-df$taxon_id[df$parent_id %in% parentstodo]
    for(i in parentstodo){
      if(sum(df$parent_id==i)==0){
        next()
      }
      df$newlabel[df$parent_id==i]<-paste(df$newlabel[df$parent_id==i],df$newlabel[df$taxon_id==i])
      df$remainder[df$taxon_id==i]<-df$abundance[df$taxon_id==i]-sum(df$abundance[df$parent_id==i])
    }
    parentstodo<-parentstodo1
  }
  df
}

deparsefile(mydata)

【讨论】：

感谢您的回复。当我在完整数据集上尝试此操作时，它运行没有错误，但余数列给出了奇怪的数字。我得到了负数，并且在应该有数字的地方得到了零（比如示例中的应变水平比例应该是 0.3444）。
我刚刚手动检查了几条测试线——它们看起来很合理。那么对于底部的那些分类群（即没有后代），您希望其余部分等于总数吗？（我刚刚改变了）我认为负数主要是由于四舍五入。
太棒了！我也通过它来检查数字和由于四舍五入导致的负数对我来说是有意义的。对于位于父子树底部的分类群（其中 taxon_id 从未出现在 parent_id 列中），余数应等于丰度。
好的，所以将 df$remainder