R包decontam学习笔记

这两天的工作都是从早上打开热心肠日报开始的,不得不说基本上每天干货满满呀,每天一篇给力的综述的节奏。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)")

发表评论