Mercurial > repos > simon-gladman > phyloseq_nmds

comparison phyloseq_nmds.R @ 1:e376a618eb9f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Updated phyloseq_nmds.R to allow for non standard BIOM files.
author simon-gladman
date Sat, 16 Jun 2018 05:03:43 -0400
parents b4606394e7ec
children 20adf95eb758
comparison
equal deleted inserted replaced
0:b4606394e7ec 1:e376a618eb9f
1 library('getopt') 1 library('getopt')
2 library('ape') 2 library('ape')
3 suppressPackageStartupMessages(library('phyloseq')) 3 suppressPackageStartupMessages(library('phyloseq'))
4 library(biomformat)
5 library(plyr)
4 Sys.setenv("DISPLAY"=":1") 6 Sys.setenv("DISPLAY"=":1")
5 library("ggplot2") 7 library("ggplot2")
6 8
7 options(warn=-1) 9 options(warn=-1)
8 10
44 kingdom_str<-options$kingdom 46 kingdom_str<-options$kingdom
45 distance<-options$distance 47 distance<-options$distance
46 keep<-options$keep 48 keep<-options$keep
47 49
48 50
51
52 ### This function accepts different two different type of BIOM file format
53 readBIOM<-function(inBiom){
54 tryCatch({
55 phyloseq_obj<-import_biom(inBiom)
56 return(phyloseq_obj)
57 },
58 error=function(e){
59 biom_obj<-read_biom(inBiom)
60
61 otu_matrix = as(biom_data(biom_obj), "matrix")
62 OTU_TABLE = otu_table(otu_matrix, taxa_are_rows=TRUE)
63
64 taxonomy_matrix = as.matrix(observation_metadata(biom_obj), rownames.force=TRUE)
65 TAXONOMY_TABLE = tax_table(taxonomy_matrix)
66
67 metadata.temp<-sample_metadata(biom_obj)
68 METADATA_TABLE<-plyr::ldply(metadata.temp, rbind)
69 rownames(METADATA_TABLE)<-as.character(METADATA_TABLE$.id)
70
71 phyloseq_obj = phyloseq(OTU_TABLE, TAXONOMY_TABLE,sample_data(METADATA_TABLE))
72 return(phyloseq_obj)
73 }
74 )
75 }
76
77
78
49 if(!is.null(options$biom)){ 79 if(!is.null(options$biom)){
50 80
51 physeq<-import_biom(options$biom) 81 #physeq<-import_biom(options$biom)
52 82 physeq<-readBIOM(options$biom)
83
53 if(length(rank_names(physeq)) == 8){ 84 if(length(rank_names(physeq)) == 8){
54 tax_table(physeq) <- tax_table(physeq)[,-1] 85 tax_table(physeq) <- tax_table(physeq)[,-1]
55 colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") 86 colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
56 } 87 }
57
58 ### select column name from sample table for nmds plot 88 ### select column name from sample table for nmds plot
59 category_type<-colnames(sample_data(physeq))[options$subset] 89 category_type<-colnames(sample_data(physeq))[options$subset]
60 90
61 ### obtain the unique value in the selected column from sample table 91 ### obtain the unique value in the selected column from sample table
62 category_option<-unique(sample_data(physeq))[,options$subset] 92 category_option<-unique(sample_data(physeq))[,options$subset]
93 #kingdom_str<-options$kingdom 123 #kingdom_str<-options$kingdom
94 124
95 125
96 ### Remove OTUs that do not appear more than 5 times in more than half the samples 126 ### Remove OTUs that do not appear more than 5 times in more than half the samples
97 ### filtering OTUs based on cutoff value (e.g., 5) 127 ### filtering OTUs based on cutoff value (e.g., 5)
98 physeq_temp =genefilter_sample(physeq, filterfun_sample(function(x) x > cutoff_value), A=0.5*nsamples(physeq)) 128 physeq_temp =genefilter_sample(physeq, filterfun_sample(function(x) x > cutoff_value), A=0.1*nsamples(physeq))
99 129
100 ### phyloseq object after filtered 130 ### phyloseq object after filtered
101 physeq_filter = prune_taxa(physeq_temp, physeq) 131 physeq_filter = prune_taxa(physeq_temp, physeq)
102 132
103 ## Transform to even sampling depth 133 ## Transform to even sampling depth