qiime 1.9.0 biom files and phyloseq

[Clayton Collings]

Qiime 1.8.0 biom files work in Phyloseq, but Qiime 1.9.0 biom files do not.  

Here is the error.

Error in fromJSON(content, handler, default.size, depth, allowComments,  : 

  invalid JSON input

I suppose that I can use differential_abundance.py now, but I like the figures that phyloseq/ggplot2 makes

as well as my MA plots and DESeq2's 2D PCoA plots based solely on OTU abundance (euclidean distances).

I wonder if QIIME could add some figure options to this module.  Also, I are the taxonomies included in the ouput OTU tables from differential_abundance.py?

I am going to look at this today.

Here is my R code that I use when implementing phyloseq.  See additional comments and additional R code below for multifactor DESeq2.

> data <- import_biom("IAvsIP.biom")

> x <- read.table("IAvsIP.txt", header=TRUE, sep="\t", row.names=1)

> x

                         Treatment  Type Animal MergedReads

Phage-Cecum-18-F11           Phage Cecum     18       65853

Phage-Cecum-14-E11           Phage Cecum     14      106799

Phage-Cecum-11-B09           Phage Cecum     11      106897

Phage-Cecum-23-H11           Phage Cecum     23      109586

Phage-Cecum-12-C10           Phage Cecum     12      113183

Antibiotic-Cecum-7-B12  Antibiotic Cecum      7       37891

Antibiotic-Cecum-17-F02 Antibiotic Cecum     17      131735

Antibiotic-Cecum-20-G02 Antibiotic Cecum     20      149945

Antibiotic-Cecum-22-H04 Antibiotic Cecum     22      293886

Antibiotic-Cecum-15-E02 Antibiotic Cecum     15      470301

> sampledata = sample_data(data.frame(x))

> data2 = merge_phyloseq(data, sampledata)

> data2

phyloseq-class experiment-level object

otu_table()   OTU Table:         [ 28733 taxa and 8 samples ]

sample_data() Sample Data:       [ 8 samples by 4 sample variables ]

tax_table()   Taxonomy Table:    [ 28733 taxa by 7 taxonomic ranks ]

> diagdds = phyloseq_to_deseq2(data2, ~ Treatment)

converting counts to integer mode

> diagdds = DESeq(diagdds, test="Wald", fitType="parametric")

estimating size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

> res = results(diagdds, cooksCutoff = FALSE)

> alpha = 0.05

> sigtab = res[which(res$padj < alpha), ]

> sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(data2)[rownames(sigtab), ], "matrix"))

> dim(sigtab)

[1] 75 13

> res = cbind(as(res, "data.frame"), as(tax_table(data2)[rownames(res), ], "matrix"))

> res = cbind(as(res, "data.frame"), as(otu_table(data2)[rownames(res), ], "matrix"))

> write.table(res, file = "DESeq2_Output_IAvsIP.txt", sep = "\t", col.names=NA)

> theme_set(theme_bw())

> scale_fill_discrete <- function(palname = "Set1", ...) {

+     scale_fill_brewer(palette = palname, ...)

+ }

> x = tapply(sigtab$log2FoldChange, sigtab$Rank2, function(x) max(x))

> x = sort(x, TRUE)

> sigtab$Phylum = factor(as.character(sigtab$Rank2), levels=names(x))

> x = tapply(sigtab$log2FoldChange, sigtab$Rank6, function(x) max(x))

> x = sort(x, TRUE)

> sigtab$Genus = factor(as.character(sigtab$Rank6), levels=names(x))

> pdf(file = "IAvsIP.pdf")

> ggplot(sigtab, aes(x=Genus, y=log2FoldChange, color=Phylum)) + geom_point(size=6) + 

+   theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))

> dev.off()

null device 

          1 

> pdf(file="IAvsIP_MA.pdf")

> plot(res$baseMean,res$log2FoldChange,ylim=c(-12,12),main="IAvsIP", log = "x", pch = 19, cex = ifelse(res$padj < .05, 0.6, 0.4), col = ifelse(res$padj < .05, "red", "black"), xlab="mean OTU abundance", ylab = "log fold change", xlim = c(1,10000))

Warning message:

In xy.coords(x, y, xlabel, ylabel, log) :

  22544 x values <= 0 omitted from logarithmic plot

> dev.off()

null device 

          1 

> rld <- rlog(diagdds)

> pdf(file="IAvsIP_PCA.pdf")

> plotPCA(rld, intgroup="Treatment")

> dev.off()

null device 

          1 

***** 

I have also done mutlifactor DESeq2.  It is also important to know what group the log2 ratios are relative to.  Note code in bold Red.

*****

library(DESeq2)

library("phyloseq")

library("ggplot2")

library("plyr")

data <- import_biom("AllCvsAllP.biom")
x <- read.table("AllCvsAllP.txt", header=TRUE, sep="\t", row.names=1)
sampledata = sample_data(data.frame(x))
data2 = merge_phyloseq(data, sampledata)
data2
diagdds = phyloseq_to_deseq2(data2, ~ Type + Treatment)
diagdds$Treatment <- relevel(diagdds$Treatment, "Control")
diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(data2)[rownames(sigtab), ], "matrix"))
dim(sigtab)
res = cbind(as(res, "data.frame"), as(tax_table(data2)[rownames(res), ], "matrix"))
res = cbind(as(res, "data.frame"), as(otu_table(data2)[rownames(res), ], "matrix"))
write.table(res, file = "DESeq2_Output_AllCvsAllP.txt", sep = "\t", col.names=NA)
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Rank2, function(x) max(x))
x = sort(x, TRUE)
sigtab$Phylum = factor(as.character(sigtab$Rank2), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Rank6, function(x) max(x))
x = sort(x, TRUE)
sigtab$Genus = factor(as.character(sigtab$Rank6), levels=names(x))
pdf(file = "AllCvsAllP.pdf")
ggplot(sigtab, aes(x=Genus, y=log2FoldChange, color=Phylum)) + geom_point(size=6) +
theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()
pdf(file="AllCvsAllP_MA.pdf")
plot(res$baseMean,res$log2FoldChange,ylim=c(-12,12),main="AllCvsAllP", log = "x", pch = 19, cex = ifelse(res$padj < .05, 0.6, 0.4), col = ifelse(res$padj < .05, "red", "black"), xlab="mean OTU abundance", ylab = "log fold change", xlim = c(1,10000))
dev.off()
rld <- rlog(diagdds)
pdf(file="AllCvsAllP_PCA.pdf")
plotPCA(rld, intgroup=c("Treatment", "Type"))
dev.off()

[Jose Carlos Clemente]

Hi Clayton,

you might want to convert your biom file to JSON, as detailed here:

http://biom-format.org/documentation/biom_conversion.html

Jose

--

---
You received this message because you are subscribed to the Google Groups "Qiime Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qiime-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Claudia]

Hi there,

I am trying to convert my biom file to a tab-delimited file using: biom convert -i out_table_filtered.biom -o table.from_biom_w_consensuslineage.txt --to-tsv --header-key taxonomy --output-metadata-id "ConsensusLineage" http://biom-format.org/documentation/biom_conversion.html

Unfortunately this command returns the following message: pyqi.core.exception.CommandError: Must specify --table-type!

I am using the MacQIIME version 1.9 which includes the biom-format version 2.1.3.

Is it possible that the command I used does not work with the biom format HDF5? I believe this is the format of my biom file http://biom-format.org/documentation/format_versions/biom-2.1.html

 

Any help?

 

Thanks!

[Jose Carlos Clemente]

Hi Claudia,

how did you generate the biom table you are trying to convert?

Jose

[Claudia]

Hi Jose,

I used the following commands before I got my filtered biom file

First I used Mothur to extract fasta and quality files: sffinfo(sff=1J_28F.sff-2J_28F.sff-3J_28F.sff-4J_28F.sff-5J_28F.sff-6J_28F.sff)

Then in Qiime I used:

split_libraries.py -b 8 -m RemovedPrimersMappingFile.txt -f 1J_28F.fasta,2J_28F.fasta,3J_28F.fasta,4J_28F.fasta,5J_28F.fasta,6J_28F.fasta -q 1J_28F.qual,2J_28F.qual,3J_28F.qual,4J_28F.qual,5J_28F.qual,6J_28F.qual -z truncate_only -o 4_Output/

 

pick_de_novo_otus.py -i seqs.fna -o 5_OTUs

 

filter_taxa_from_otu_table.py -i otu_table.biom -o otu_table_filtered.biom -n Unassigned,c__Chloroplast,f__mitochondria

 

Claudia

On Thursday, March 12, 2015 at 1:58:39 PM UTC-3, Clayton Collings wrote:

[Colin Brislawn]

Hello Claudia,

Like the error said, you "must specify --table-type". I think you can solve the problem by adding this to your command: --table-type="OTU table"

Converting tables is a little differenet in qiime 1.9.0. I found the answer to your questions here:

http://biom-format.org/documentation/biom_conversion.html

I hope that helps,

Colin

On Thursday, March 12, 2015 at 4:12:10 PM UTC-7, Claudia wrote:

Hi there,

I am trying to convert my biom file to a tab-delimited file using: biom convert -i out_table_filtered.biom -o table.from_biom_w_consensuslineage.txt --to-tsv --header-key taxonomy --output-metadata-id "ConsensusLineage" http://biom-format.org/documentation/biom_conversion.html

Unfortunately this command returns the following message: pyqi.core.exception.CommandError: Must specify --table-type!

I am using the MacQIIME version 1.9 which includes the biom-format version 2.1.3.

Is it possible that the command I used does not work with the biom format HDF5? I believe this is the format of my biom file http://biom-format.org/documentation/format_versions/biom-2.1.html

 

Any help?

 

...

[Claudia]

Dear Colin,

Thanks for your reply. Yes, your comment helped me, thanks a lot!

Claudia

On Thursday, March 12, 2015 at 1:58:39 PM UTC-3, Clayton Collings wrote: