一个 R 包带你挖掘宏基因组公共数据库
本文转自“生信菜鸟团”，已获授权背景介绍目前虽然已经有越来越多的宏基因组数据被上传到公共数据库，但由于种种原因（比如，需要消耗更多计算资源，不同研究之间的表型信息定义或格式不同等等），导...
本文转自“生信菜鸟团”，已获授权
背景介绍
目前虽然已经有越来越多的宏基因组数据被上传到公共数据库，但由于种种原因（比如，需要消耗更多计算资源，不同研究之间的表型信息定义或格式不同等等），导致很多数据并没有被充分利用。为了克服这些挑战，curatedMetagenomicData 包应运而生（或者说它更是一个数据库）。目前已有来自 46 个数据集的八千多个样本纳入其中，这里所有的原始宏基因组测序数据都使用一套统一的流程进行上游分析（使用 MetaPhlAn2 流程进行菌群丰度注释，使用 HUMAnN2 流程进行功能分析），所有数据集的表型信息都重新使用一套相同的标准进行注释和重新定义，最后将每一个数据集打包为 ExpressionSet 对象，这样一来就大大降低了挖掘宏基因组公共数据的门槛，即使是从前没接触过宏基因组数据的小伙伴也能进行下游进一步探索。
curatedMetagenomicData 包含哪些数据：
来自 46 个数据集的 8184 个样本，主要是人类肠道的样本，也包括一些人类微生物组计划（HMP）中其他的身体部位
对于这八千多个样本，包括了所有宏基因组数据（菌群分类谱，Marker 信息，基因家族丰度，通路覆盖度及丰度）和经过统一定义的表型数据，并集成为 Bioconductor 中的 ExpressionSet 对象
包含了大约 80 种 metadata 数据，所有的条目都有标准的注释
这些数据使用 MetaPhlAn2 流程进行菌群丰度注释，使用 HUMAnN2 流程进行功能分析
这项工作需要分析大约 100T 的原始数据
使用教程
安装 curatedMetagenomicData 包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ExperimentHubData")
library(curatedMetagenomicData)
查看所有可用的样本和 metadata
所有样本的 metadata 都已经包含在 combined_metadata 这一个表格中：
?combined_metadata
View(combined_metadata)
table(combined_metadata$antibiotics_current_use)
##
## no yes
## 3029 661
table(combined_metadata$disease)
##
## AD
## 10
## AD;AR
## 16
## AD;AR;asthma
## 8
## AD;asthma
...
## skininf
## 2
## stomatitis
## 2
## suspinf
## 1
## tonsillitis
## 3
下载数据集
挑选好自己感兴趣的样本就可以开始下载数据啦~ 使用 curatedMetagenomicData() 函数下载数据集，会返回一个 ExpressionSet 对象。
下载单个数据集：
loman <- curatedMetagenomicData("LomanNJ_2013.metaphlan_bugs_list.stool", dryrun = FALSE)
## Working on LomanNJ_2013.metaphlan_bugs_list.stool
## snapshotDate(): 2019-04-29
## see ?curatedMetagenomicData and browseVignettes('curatedMetagenomicData') for documentation
## downloading 0 resources
## loading from cache
## 'EH1816 : 1816'
loman
## List of length 1
## names(1): LomanNJ_2013.metaphlan_bugs_list.stool
loman.eset <- loman[[1]]
loman.eset
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 736 features, 43 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: OBK1122 OBK1196 ... OBK5066 (43 total)
## varLabels: subjectID body_site ... NCBI_accession (19 total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## pubMedIds: 23571589
## Annotation:
同时下载多个数据集：
oral <- c("BritoIL_2016.metaphlan_bugs_list.oralcavity",
"Castro-NallarE_2015.metaphlan_bugs_list.oralcavity")
esl <- curatedMetagenomicData(oral, dryrun = FALSE)
esl
esl[[1]]
esl[[2]]
使用 dryrun = TRUE 参数查看所有可用的数据集：
curatedMetagenomicData("*metaphlan_bugs_list.stool*", dryrun = TRUE)
## Dry run: see return values for datasets that would be downloaded. Run with `dryrun=FALSE` to actually download these datasets.
## [1] "AsnicarF_2017.metaphlan_bugs_list.stool"
## [2] "BackhedF_2015.metaphlan_bugs_list.stool"
## [3] "Bengtsson-PalmeJ_2015.metaphlan_bugs_list.stool"
## [4] "BritoIL_2016.metaphlan_bugs_list.stool"
## [5] "ChengpingW_2017.metaphlan_bugs_list.stool"
...
## [37] "VogtmannE_2016.metaphlan_bugs_list.stool"
## [38] "WenC_2017.metaphlan_bugs_list.stool"
## [39] "XieH_2016.metaphlan_bugs_list.stool"
## [40] "YuJ_2015.metaphlan_bugs_list.stool"
## [41] "ZellerG_2014.metaphlan_bugs_list.stool"
合并多个数据集
接下来的操作会把上面下载的两个口腔微生态的宏基因组数据集合并到一个单独的 ExpressionSet对象中。
eset <- mergeData(esl)
eset
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 1914 features, 172 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: BritoIL_2016.metaphlan_bugs_list.oralcavity:M1.1.SA
## BritoIL_2016.metaphlan_bugs_list.oralcavity:M1.10.SA ...
## Castro-NallarE_2015.metaphlan_bugs_list.oralcavity:ES_211 (172
## total)
## varLabels: subjectID body_site ... studyID (19 total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:
使用 phyloseq 包进行菌群分析
phyloseq 包是一个集菌群丰度数据导入，存储，分析和可视化于一体的工具。它不但利用了 R 中许多经典的工具进行生态学和系统发育分析（例如：vegan，ade4，ape， picante 等），同时还结合 ggplot2 以轻松生成发表级别的可视化图表，可谓是分析宏基因组数据的利器。所以，对于菌群丰度信息 curatedMetagenomicData 包十分贴心地提供了从 ExpressionSet 对象转为 phyloseq 对象的函数 ExpressionSet2phyloseq()，这样我们就可以直接使用 phyloseq 包进行菌群多样性下游分析了。
suppressPackageStartupMessages(library(phyloseq))
loman.pseq = ExpressionSet2phyloseq( loman.eset, phylogenetictree = TRUE)
loman.pseq
phyloseq-class experiment-level object
otu_table() OTU Table: [ 237 taxa and 43 samples ]
sample_data() Sample Data: [ 43 samples by 19 sample variables ]
tax_table() Taxonomy Table: [ 237 taxa by 8 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 237 tips and 236 internal nodes ]
更多关于 phyloseq 包的使用教程可参见我之前的教程~
使用 ExpressionSet 对象
由于 S4 类的集成性及其绑定数据和元数据的能力，因此所有数据集都被表示为 ExpressionSet 对象。
使用 experimentData() 查看实验设计相关的数据：
experimentData(loman.eset)
## Experiment data
## Experimenter name: Loman NJ, Constantinidou C, Christner M, Rohde H, Chan JZ, Quick J, Weir JC, Quince C, Smith GP, Betley JR, Aepfelbacher M, Pallen MJ
## Laboratory: [1] Institute of Microbiology and Infection, University of Birmingham, Birmingham, England.
## Contact information:
## Title: A culture-independent sequence-based metagenomics approach to the investigation of an outbreak of Shiga-toxigenic Escherichia coli O104:H4.
## URL:
## PMIDs: 23571589
##
## Abstract: A 308 word abstract is available. Use 'abstract' method.
## notes:
## Sequencing platform:
## IlluminaHiSeq
使用 pData() 查看样本信息数据：
head(pData(loman.eset))
## subjectID body_site antibiotics_current_use study_condition disease
## OBK1122 OBK1122 stool <NA> STEC STEC
## OBK1196 OBK1196 stool <NA> STEC STEC
## OBK1253 OBK1253 stool <NA> STEC STEC
## OBK2535 OBK2535 stool <NA> STEC STEC
## OBK2535b OBK2535b stool <NA> STEC STEC
## OBK2638 OBK2638 stool <NA> STEC STEC
## age_category country non_westernized DNA_extraction_kit
## OBK1122 adult DEU no Qiagen
## OBK1196 adult DEU no Qiagen
## OBK1253 adult DEU no Qiagen
## OBK2535 adult DEU no Qiagen
## OBK2535b adult DEU no Qiagen
## OBK2638 adult DEU no Qiagen
## number_reads number_bases minimum_read_length median_read_length
## OBK1122 464060 69609000 150 150
## OBK1196 30181380 4527207000 150 150
## OBK1253 65704818 9855722700 150 150
## OBK2535 43597736 6539660400 150 150
## OBK2535b 22253314 3337997100 150 150
## OBK2638 1105492 165823800 150 150
## days_after_onset stec_count shigatoxin_2_elisa stool_texture
## OBK1122 NA <NA> <NA> <NA>
## OBK1196 NA <NA> <NA> <NA>
## OBK1253 NA <NA> <NA> <NA>
## OBK2535 3 moderate positive smooth
## OBK2535b NA <NA> <NA> <NA>
## OBK2638 5 high positive bloody
## curator NCBI_accession
## OBK1122 Paolo_Manghi ERR262939;ERR262938
## OBK1196 Paolo_Manghi ERR262941;ERR262940
## OBK1253 Paolo_Manghi ERR262942;ERR260476
## OBK2535 Paolo_Manghi ERR260478;ERR260477
## OBK2535b Paolo_Manghi ERR260479
## OBK2638 Paolo_Manghi ERR262943;ERR260480
更多关于 ExpressionSet 对象的内容可以看看曾老师的教程：ExpressionSet 对象简单讲解 http://www.bio-info-trainee.com/1510.html 。
Reference
http://bioconductor.org/packages/release/data/experiment/vignettes/curatedMetagenomicData/inst/doc/curatedMetagenomicData.html
猜你喜欢
phyloseq | 用 R 分析微生物组数据及可视化（一）
phyloseq | 用 R 分析微生物组数据及可视化（二）
phyloseq | 用 R 分析微生物组数据及可视化（三）
猜你喜欢
10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
文献阅读 热心肠 SemanticScholar Geenmedical
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。
更多推荐
4
0- 0
已为社区贡献5条内容
所有评论(0)