-
Notifications
You must be signed in to change notification settings - Fork 1
/
phyloseq.Rmd
56 lines (42 loc) · 1.64 KB
/
phyloseq.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
---
title: "Phyloseq"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library("phyloseq")
library("ggplot2")
library("biomformat")
```
#adapted from https://github.com/BikLab/ME_microbiome_2017/blob/80691956fe65fb51221ac25539908c0d5e9a8897/scripts/betadiv-and-filtering/PCoAs-in-R.sh
#https://joey711.github.io/phyloseq-demo/import-biom-sd-example.html
#https://github.com/joey711/phyloseq/issues/167
#could use OPEN OTU
#could also have chosen intermediates/OTU_SILVA/CSS_normalized_otu_table.biom
#./SILVA_119_QIIME_release/97_FastTree_trees/Silva_119_rep_set97_aligned_16S_only_pfiltered.tre
```{r renametaxalevels}
if(!exists("biomfile")){biomfile<-"intermediates/OTU_SILVA_CLOSED/otu_table.biom"}
if(!exists("treefile")){treefile<-"SILVA_119_QIIME_release/97_FastTree_trees/Silva_119_rep_set97_aligned_16S_only_pfiltered.tre"}
ps<-import_biom(BIOMfilename=biomfile,treefilename = treefile)
#"intermediates/OTU_SILVA_CLOSED/rep_set.tre")
#"intermediates/OTU_SILVA_CLOSED/otu_table_mc2_w_tax.biom"
#https://github.com/joey711/phyloseq/issues/436
data(GlobalPatterns)
t <- data.frame(tax_table(ps))
colnames(t) <- rank_names(GlobalPatterns)
tax_table(ps) <- tax_table(as.matrix(t))
```
```{r process}
#should we use this?
ps.log <- transform_sample_counts(ps, function(x) log(1 + x))
bmsd <- import_qiime_sample_data("mapfile.txt")
p0 <- merge_phyloseq(ps, bmsd)
#3 min
out.uf.log <- ordinate(p0, method = "MDS", distance = "unifrac")
evals <- out.uf.log$values$Eigenvalues
```
```{r pressure, echo=FALSE}
plot_ordination(p0, out.uf.log, color = "Type") +
labs(col = "Type") +
coord_fixed(sqrt(evals[2] / evals[1]))
```