物种丰度处理

利用shell对物种丰度表进行处理：

1.原始丰度表如下：

image.png
每列由空格隔开，物种由分号隔开
处理如下：（注：处理完需要在物种分类手动加入界门纲目科属种）

#因为物种分类也有空格，因此需要和丰度分开处理，首先去掉最前面的空格，然后以分号分隔取丰度这一边，然后将空格转换为制表符
cat tax.xls | sed 's/^\ //g' | awk -F ";" '{print $1}' | sed 's/\ /\t/g' > tax11.txt 
#然后取物种分类这边，将分号替换为制表符
cat tax.xls | sed 's/^\ //g' | awk -F " " '{for (i=37;i<=NF;i++)printf("%s ", $i);print ""}' | sed 's/\;/\t/g' > tax_z.txt

处理完物种分类如下：

image.png

然后将两个文件导入R语言：

ab <- read.delim("tax11.txt",sep = "\t",header = T,row.names= NULL)
species <- read.delim("tax_z.txt",sep = "\t",header = T,row.names= NULL)
species_ab <- cbind(ab,species)
species_ab <- species_ab[, ！colnames(species_ab) == "Taxonomy"]
#按照不同的界提取子表格
k__Archaea <- subset(species_ab,species_ab$Kingdom == 'k__Archaea')
k__Bacteria <- subset(species_ab,species_ab$Kingdom == 'k__Bacteria')
k__Eukaryota <- subset(species_ab,species_ab$Kingdom == 'k__Eukaryota')
k__unidentified <- subset(species_ab,species_ab$Kingdom == 'k__unidentified')
k__Viruses <- subset(species_ab,species_ab$Kingdom == 'k__Viruses')

然后就能对这些界的物种进行进一步的处理，比如，想得到细菌界所有纲的累积丰度：
之前在diamond比对以及结果数据中写过一个叠加小代码。
然而，发现这代码运行在R中运行非常慢，不适合大数据处理。
今天换一个思路写下代码，一般几百Mb大小的丰度文件，也就一分钟不到就运行完了。
附上代码：（注：文件处理格式，前几列为样本丰度，后几列是物种分类）

#参数data指数据需要处理的数据框，lie指需要对哪个分类等级进行叠加，sample指样本数目。
diejia <- function(data = data.frame(),lie = "Class" ,sample = 1){
  list1 <- list()
  x2 <- data.frame()
  sum_ab <- data.frame()
  xx <- data.frame()
  xxx <- data.frame()
  engame <- data.frame()
  z <- levels(factor(data[,lie])) #提取指定列的不同的分类名
  for(i in z){
    x2 <- subset(data,data[,lie] %in% i)
    list1[[i]] <- x2
  }    #第一个循环将相同名字的所有行形成新的数据框，并加入到列表中
  for (o in 1:length(z)) {
    specie <- list1[[o]][1,(sample+1):ncol(list1[[o]])]
    xxx <- list1[[o]][,1:sample]
    sum_ab <- colSums(xxx)
    sum_ab <- as.data.frame(t(as.data.frame(sum_ab)))
    xx <- cbind(sum_ab,specie)
    engame <- rbind(engame,xx)
  }   #第二个循环，从列表中取每一个数据框，对他们求每列的和，也就是该物种在样本中的所有丰度。然后再将这个丰度与每一个数据框所处理 的丰度合并成一个新的数据框。
  return(engame)
} 
#写完代码就可以运行了
Bacteria_class <- diejia(data = k__Bacteria , lie = "Class" ,sample = 36)

最后得到每个纲的丰度表格

image.png