这两天的工作都是从早上打开热心肠日报开始的，不得不说基本上每天干货满满呀，每天一篇给力的综述的节奏。ps,发现ubiome的科学家和志愿者也在维护着一个国外版的热心肠日报，地址放在这：https://microbiomedigest.com/ 值得点赞，英语好的同鞋可以参考参考，不知道国内的灵感是不是来自这个，但是平心而论，国内的做得已经水平不次于国外的了。
很感谢这些人，为大众做着科普的事情，有许多人把今天的生命科学比做几十年前的计算机产业，当年，计算机还是庞然大特，像供神一样地供在一间特殊的房间，多么像今天的测序仪呀。而牛津纳米孔，直接把测序仪缩小到了手持设备的级别，简直从当年的大型机飞跃到了今天的智能手机呀！相信未来，测序仪会像手机一样人手一部来检测自己的各项基因、蛋白和代谢指标，这才是精准医学嘛。
好了，回到肠道微生物，许多人甚至还不知道他的存在和重要性，甚至包括生命科学的学生，今天，日益增多的研究成果和检测企业让它走进寻常百姓家，速度甚至比基因检测更加快速。但是，仍然存在着变化较大，娱乐性偏高，以及科学性存疑。相信随着研究的深入，大的人群队列研究结果的出炉，肠道微生物检测一定会更加有用，更加深入人心。
扯的有点远，哈哈，再回到这个R包，它是用来去除16S等Marker基因和宏基因组测序中的样本污染序列的工具。热心肠日报上对其编译文章的总结在这：
① 污染序列是指由试剂和实验操作等因素造成的样品中没有的序列，它们严重影响测序结果的准确性；② R包 decontam 可根据两种分类模式识别和去除污染序列：一是通常出现在较低浓度样品中，二是存在于阴性对照中；③ 基于口腔菌群扩增子数据和早产儿-阴道菌群宏基因组数据对该方法进行验证，结果表明可有效去除污染物DNA序列，提高测序质量；④ 该R包可与现有数据分析流程整合，帮助研究人员更高效、低成本的获取准确的微生物群落概况。
下面是它的安装和使用方法学习：
最简单的安装方法来自它的github里的README，使用devtools安装。代码如下：
install.packages("devtools") #如果之前装过就不用运行了
library(devtools)
devtools::install_github("benjjneb/decontam")
后面的代码先简单记录下，R还没入门是个比较愁人的问题。
#4.setting up
library(phyloseq); packageVersion("phyloseq")
library(ggplot2); packageVersion("ggplot2")
library(decontam); packageVersion("decontam")
ps <- readRDS(system.file("extdata", "MUClite.rds", package="decontam"))
ps
#5.Inspect Library Sizes
head(sample_data(ps))
df <- as.data.frame(sample_data(ps)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_or_Control)) + geom_point()
# 6 Identify Contaminants - Frequency
contamdf.freq <- isContaminant(ps, method="frequency", conc="quant_reading")
head(contamdf.freq)
table(contamdf.freq$contaminant)
head(which(contamdf.freq$contaminant))
plot_frequency(ps, taxa_names(ps)[c(1,3)], conc="quant_reading") +
xlab("DNA Concentration (PicoGreen fluorescent intensity)")
set.seed(100)
plot_frequency(ps, taxa_names(ps)[sample(which(contamdf.freq$contaminant),3)], conc="quant_reading") +
xlab("DNA Concentration (PicoGreen fluorescent intensity)")
ps
ps.noncontam <- prune_taxa(!contamdf.freq$contaminant, ps)
ps.noncontam
#7.Identify Contaminants - Prevalence
sample_data(ps)$is.neg <- sample_data(ps)$Sample_or_Control == "Control Sample"
contamdf.prev <- isContaminant(ps, method="prevalence", neg="is.neg")
table(contamdf.prev$contaminant)
head(which(contamdf.prev$contaminant))
contamdf.prev05 <- isContaminant(ps, method="prevalence", neg="is.neg", threshold=0.5)
table(contamdf.prev05$contaminant)
# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(ps, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$Sample_or_Control == "Control Sample", ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$Sample_or_Control == "True Sample", ps.pa)
# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
contaminant=contamdf.prev$contaminant)
ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")