其他
R语言微生物群落数据分析:从原始数据到结果（No1）
微生物群落数据分析工作流程：从原始数据到群落分析
Bioconductor workflow for microbiome data analysis: from raw reads to community analyses
本教程分析所需R包及其识别安装
设定出图主个主题
随机种子设置，和R环境设置
准备OTU表格预处理和变换，并存为函数，方便随时中断和继续
1. 清除环境，设定路径，载入包
2. 下载数据，提取数据和样品名称
3. 测序质量评估，数据初步过滤
4. 去冗余，计算错误模型。
5. dada提取物种序列
6. 序列拼接，双端序列合并，去除嵌合体
物种注释
构建系统发育树及其进化树裁剪
构建phyloseq对象
微生物群落数据分析工作流程：从原始数据到群落分析
Bioconductor workflow for microbiome data analysis: from raw reads to community analyses
地址: https://f1000research.com/articles/5-1492/v1
翻译及其代码测试：文涛
备注：这篇文章发表于:F1000Research，是一个涵盖所有生命科学领域的全球开放获取期刊。在获得编辑部基本的科学性以及完整性审核后，未经审稿人审稿的论文会立即被刊发在网站上。随后，来自受邀审稿人的评议意见也会与论文列在一起公开发布（包括评审专家的姓名以及评议报告）。作者可以上传文章的新版本，以回应审稿者的评议。一旦通过了同行评议，论文就会被编入PubMed,Scopus和其它数据库的索引。作为F1000Research采取的开放科学模式的一部分，每篇文章背后的数据也会被发布，并且可以自由下载，以便于审稿者和其他研究人员进行分析。F1000Research解决了当今困扰科学出版的主要问题，即研究的及时传播，同行评议和数据共享。
备注：作者称OTU表格为 RSVs (Ribosomal Sequence Variants)表格，一下两种叫法通用。
这篇文章发表于2016年。本篇章内我们将在R语言中整合目前出现的大量的基于扩增子测序的生物信息学分析手段。包括参数和非参数方法。使用的R包主要有：dada2, phyloseq, DESeq2, ggplot2 and vegan。主要涉及测序原属数据处理，数据清洗以及可视化。我们还提供了监督分析方法，像随机森林，偏最小二乘，线性模型等；非参数方法有网络分析。
为了更好的熟练以及运用扩增子测序数据的分析流程，于2019年2月15日，我将这篇完整的流程进行一个演绎，和中文的解释说明，一份很好的教程用来学习R语言生物信息平台和微生物组学数据梳理入门。
扩增子生物信息流程：从原始数据到OTU表格
这部分将是一个，从原始测序序列到微生物群落特征表格（OTU）物种注释 进化树的构建完整的生物信息学过程。这里使用的数据：这360份粪便样本是在12只小鼠出生后的第一年按照时间梯度收集的，目的是研究小鼠体内微生物群落系统发育和稳定性。
本教程分析所需R包及其识别安装
.cran_packages <- c("knitr", "phyloseqGraphTest", "phyloseq", "shiny",
"miniUI", "caret", "pls", "e1071", "ggplot2", "randomForest",
"vegan", "plyr", "dplyr", "ggrepel", "nlme",
"reshape2", "devtools", "PMA","structSSI","ade4",
"igraph", "ggnetwork", "intergraph", "scales")
.github_packages <- c("jfukuyama/phyloseqGraphTest")
.bioc_packages <- c("phyloseq", "genefilter", "impute")
# Install CRAN packages (if not already installed)
.inst <- .cran_packages %in% installed.packages()
if (any(!.inst)) {
install.packages(.cran_packages[!.inst], repos = "http://cran.rstudio.com/")
}
.inst <- .github_packages %in% installed.packages()
if (any(!.inst)) {
devtools::install_github(.github_packages[!.inst])
}
.inst <- .bioc_packages %in% installed.packages()
if (any(!.inst)) {
source("http://bioconductor.org/biocLite.R")
biocLite(.bioc_packages[!.inst])
}
设定出图主个主题
## ---- ggplot-theme ----
library("ggplot2")
theme_set(theme_bw())
min_theme <- theme_update(panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks = element_blank(),
legend.title = element_text(size = 8),
legend.text = element_text(size = 6),
axis.text = element_text(size = 6),
axis.title = element_text(size = 8),
strip.background = element_blank(),
strip.text = element_text(size = 8),
legend.key = element_blank())
随机种子设置，和R环境设置
# ---- misc-setup ----
set.seed(10)
width来调整向量，矩阵的输出宽度
options(width=100)
准备OTU表格预处理和变换，并存为函数，方便随时中断和继续
首先将指定rds文件，并读取，在我们通过DADA得到OTU表格之后可以进行分段分析群落数据
## ---- get-raw-data ----
raw_data <- function() {
ps_file <- "D:/metagenomics/No2_phyloseq/a2_phyloseq_mice_full_story/ps.rds"
readRDS(ps_file)
}
## ---- get-preprocessed-data ----
preprocessed_data <- function() {
ps <- raw_data()
ps <- prune_samples(rowSums(otu_table(ps)) > 1000, ps)
sample_data(ps)$librarysize <- log10(rowSums(otu_table(ps)))
sample_data(ps)$age_binned <- cut(sample_data(ps)$age, breaks = c(0, 100, 200, 400))
ps0 <- ps
pslog <- transform_sample_counts(ps, function(x) log(1 + x))
# detected in preprocessing.R
outliers <- c("F5D165", "F6D165", "M3D175", "M4D175", "M5D175", "M6D175")
ps0 <- prune_samples(!(sample_names(ps0) %in% outliers), ps0)
pslog <- prune_samples(!(sample_names(pslog) %in% outliers), pslog)
list(ps0 = ps0, ps = ps, pslog = pslog)
}
## ---- setup-example ----
setup_example <- function(pkgs) {
# clear workspace, except data getters
all_obj <- ls(envir = .GlobalEnv)
all_obj <- setdiff(all_obj, c("raw_data", "preprocessed_data", "setup_example"))
rm(list = all_obj, envir = .GlobalEnv)
# add packages and data to workspace
sapply(pkgs, require, character = TRUE)
# slight problem with phylo object class...
attach(preprocessed_data()) # from data.R
}
1. 清除环境，设定路径，载入包
#清空内存#######
rm(list=ls())
#phyloseq路径不可用中文
path = 'D:/metagenomics/No2_phyloseq/a2_phyloseq_mice_full_story/'
##phyloseq一篇完整的故事
.cran_packages <- c("ggplot2", "gridExtra")
.bioc_packages <- c("dada2", "msa", "phyloseq")
# Load packages into session
sapply(c(.cran_packages, .bioc_packages), require, character.only = TRUE)
set.seed(100)
2. 下载数据，提取数据和样品名称
原始数据下载地址：http://www.mothur.org/MiSeqDevelopmentData/StabilityNoMetaG.tar
这部分命令我们学习到了字符串处理函数。
miseq_path = 'D:/metagenomics/No2_phyloseq/a2_phyloseq_mice_full_story/'
#提取原始文件路径及其样品名称
##列出文件夹内容并排序赋予fns
fns <- sort(list.files(miseq_path, full.names = TRUE))
fnFs <- fns[grepl("R1", fns)]
fnRs <- fns[grepl("R2", fns)]
sample.names <- sapply(strsplit(basename(fnFs), "_S"), `[`, 1)
sample.names
3. 测序质量评估，数据初步过滤
首先构建过滤文件和过滤后各样品文件名。我们通过filterAndTrim命令对正向序列和反向序列分别裁剪至245和160bp（truncLen=c(240,160)），最大N的数量设置为0（maxN=0），maxEE参数表示允许在条reads中期望错误的最大值（参见“一种比平均质量更好的过滤策略”：https://academic.oup.com/bioinformatics/article/31/21/3476/194979）。
#质量检测
#从样品中随机选取三个样品
ii <- sample(length(fnFs), 3)
for(i in ii) { print(plotQualityProfile(fnFs[i]) + ggtitle("Fwd")) }
for(i in ii) { print(plotQualityProfile(fnRs[i]) + ggtitle("Rev")) }
#裁剪和过滤
#新建过滤文件路径
filtpath <- file.path(miseq_path, "filtered")
filtFs <- file.path(miseq_path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
names(filtFs) <- sample.names
filtRs <- file.path(miseq_path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtRs) <- sample.names
for(i in seq_along(fnFs)) {
fastqPairedFilter(c(fnFs[[i]], fnRs[[i]]),
c(filtFs[[i]], filtRs[[i]]),
trimLeft=10, truncLen=c(245, 160),
maxN=0, maxEE=2, truncQ=2,
compress=TRUE)
}
4. 去冗余，计算错误模型。
#去除重复序列,共1.8G，两者同时运行，16G内存够用
derepFs <- derepFastq(filtFs)
derepRs <- derepFastq(filtRs)
sam.names = sample.names
#sam.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
names(derepFs) <- sam.names
names(derepRs) <- sam.names
##计算错误模型
?dada#如果selfinclude = TRUE，该算法将在样本推理和错误率估计之间交替，直到收敛
?learnErrors#之前我们使用的是该函数进行错误模型的学习，与本次学习类似，但是该方法将全部样品pool起来，dada（默认pool=F）可以指定一批样品进行学习
ddF <- dada(derepFs[1:40], err=NULL, selfConsist=TRUE)
ddR <- dada(derepRs[1:40], err=NULL, selfConsist=TRUE)
plotErrors(ddF)
plotErrors(ddR)
5. dada提取物种序列
#开始dada，pool=TRUE指定了将所有样品放到一起做样品估计
dadaFs <- dada(derepFs, err=ddF[[1]]$err_out, pool=TRUE)
dadaRs <- dada(derepRs, err=ddR[[1]]$err_out, pool=TRUE)
6. 序列拼接，双端序列合并，去除嵌合体
#序列拼接
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs)
#将人工添加的群落去除，构建ASV表格
seqtab.all <- makeSequenceTable(mergers[!grepl("Mock", names(mergers))])
#去除嵌合体
seqtab <- removeBimeraDenovo(seqtab.all)
物种注释
#使用RDP数据库注释
ref_fasta <- tempfile()
download.file("http://benjjneb.github.io/dada2/rdp_train_set_14.fa.gz",
destfile = ref_fasta)
taxtab <- assignTaxonomy(seqtab, refFasta = ref_fasta)
colnames(taxtab) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus")
构建系统发育树及其进化树裁剪
这里提供R语言构建进化树的ape包，及其进化树构建方法
#构建系统发育树及其进化树裁剪
seqs <- getSequences(seqtab)
names(seqs) <- seqs
#进行序列比对, msa提供了三种序列比对的方法的接口 ‘ClustalW’, ‘ClustalOmega’, and ‘MUSCLE’
mult <- msa(seqs, method="ClustalW", type="dna", order="input")
#进行序列比对, msa提供了三种序列比对的方法的接口 ‘ClustalW’, #‘ClustalOmega’, and ‘MUSCLE’
mult <- msa(seqs, method="ClustalW", type="dna", order="input")
#phangorn包用来构建发育树
library("phangorn")
phang.align <- as.phyDat(mult, type="DNA", names=getSequence(seqtab))
?as.phyDat
dm <- dist.ml(phang.align)
treeNJ <- NJ(dm) # Note, tip order != sequence order
fit = pml(treeNJ, data=phang.align)
## negative edges length changed to 0!
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
rearrangement = "stochastic", control = pml.control(trace = 0))
detach("package:phangorn", unload=TRUE)
构建phyloseq对象
##数据处理完毕，进行下游分析
#构建mapping文件
mimarks_path <- file.path(tempdir(), "MIMARKS.csv")
download.file("https://www.dropbox.com/s/wpxnvsui4mjj8z3/MIMARKS_Data_combined.csv?raw=1",
destfile = mimarks_path)
samdf <- read.csv(mimarks_path, header=TRUE)
samdf$SampleID <- paste0(gsub("00", "", samdf$host_subject_id), "D", samdf$age-21)
samdf <- samdf[!duplicated(samdf$SampleID),] # Remove dupicate entries for reverse reads
rownames(seqtab) <- gsub("124", "125", rownames(seqtab)) # Fixing an odd discrepancy
all(rownames(seqtab) %in% samdf$SampleID) # TRUE
rownames(samdf) <- samdf$SampleID
keep.cols <- c("collection_date", "biome", "target_gene", "target_subfragment", "host_common_name", "host_subject_id", "age", "sex", "body_product", "tot_mass", "diet", "family_relationship", "genotype", "SampleID")
samdf <- samdf[rownames(seqtab), keep.cols]
##构建phyloseq对象，这里直接进行导入，但是我们可以包存后，在进行读取到R中
ps <- phyloseq(tax_table(taxa),
sample_data(samdf),
otu_table(seqtab, taxa_are_rows = FALSE),
phy_tree(fitGTR$tree))
#或者我们直接将phyloseq进行保存，直接调取
saveRDS(ps, "D:/metagenomics/No2_phyloseq/a2_phyloseq_mice_full_story/ps.rds")
library("phyloseq")
library("ggplot2")
library("gridExtra")
ps0 = readRDS("ps.rds")
ps0