首先导入输入文件物种某水平的分类表
gene <- read.delim('../nr/Phylum.txt',row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
图片.png
百分比转化
方式一
gene_precent1 <- as.data.frame(apply(gene, 2, function(x){x/sum(x)}))
方式二
gene_precent2 <- as.data.frame(t(t(gene)/colSums(gene,na=T))*100)
colSums(gene_precent2)
图片.png
在微生物组数据分析中，样品分析之前我们经常需要对微生物组的丰度进行筛选,
- 过滤在任何样本中百分比小于1%的物种
gene_filter <-data.frame(gene_precent1[apply(gene_precent1,1,max)>0.01,])
- 保留在任何样品中百分比大于 1%的物种
gene_filter <-data.frame(gene_precent1[apply(gene_precent1,1,min)>0.001,])
- 过滤样品平均相对丰度小于1%的物种
gene_filter2 <- data.frame(gene_precent1[which(apply(gene_precent1, 1, function(x){mean(x)})
>0.01),], check.names=F)
另一种方法
gene_filter <- gene_precent1[which(rowMeans(gene_precent1) >= 0.01), ]
- 只保留相对丰度总和高于 0.005 的属，换成rowSums即可
gene_filter <- gene_precent1[which(rowSums(gene_precent1) >= 0.005), ]
- 过滤在一半或者大于一半样品中丰度为0的物种
cutoff = .5
gene_filter <- data.frame(gene_precent1[which(apply(gene_precent1, 1, function(x){length(which
(x!= 0))/length(x)}) >= cutoff),])
提一点, 代码中x!=0其实可以换为x大于等于某个值，就代表过滤在一半或者大于一半样品中丰度大于等于多少的物种，注意变通~~~
图片.png
另外，在有的文献中还有是过滤每组中至少一半的样品丰度丰度大于0.1%，也就是说当你有俩组每组6个样品的情况下，你得保证每组都是至少有3个样品的丰度大于0.1%。其实也很简单，我们分别在各组中去执行上述代码，最后筛选到的物种再合并一下就OK了,用到了union函数，它的功能是会整合出现在x数据框中或y数据框中的数据，同时去除了两个数据框中重复的部分。
cutoff = 0.5
gene1 <- data.frame(gene_precent1[,1:6][which(apply(gene_precent1[,1:6], 1, function(x){length(which
(x>=0.001))/length(x)}) > cutoff),])
gene2 <- data.frame(gene_precent1[,7:12][which(apply(gene_precent1[,7:12], 1, function(x){length(which
(x>=0.001))/length(x)}) > cutoff),])
gene_filter1 <- gene_precent1[union(rownames(gene1),rownames(gene2)),]
提供另一种方法过滤在一半或者大于一半样品中丰度为0的物种
gene_filter <- gene_precent1
gene_filter[gene_filter >0] <- 1
gene_filter <- gene_precent1[which(rowSums(gene) >= ncol(gene_precent1)/2), ]
若对代码有问题或者有更便捷的方式可以留言回复，现在简书不能发送公众号二维码，欢迎大家搜索BioparaMeta关注，会定时分享一些生信小技巧，大家一起交流进步~~
