-
Notifications
You must be signed in to change notification settings - Fork 0
/
02.metagenome_pred.R
71 lines (55 loc) · 2.79 KB
/
02.metagenome_pred.R
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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
library(phyloseq)
library(ShortRead)
#################################################
#1. Prepare filies for picrust2 and run pipeline#
#################################################
#load phyloseq object with 16S rRNA data
load("/media/kreising/DATA/data/VLASTOVKY/PHYLOSEQ_rerun/HIRUNDO_16S.taxo.phylo.R")
#save ASVs sequences as fasta
REFS<-refseq(HIRUNDO_16S.taxo.phylo)
writeFasta(REFS,file = "/media/kreising/DATA/data/VLASTOVKY/PICRUST2/REF.fasta")
#save ASVs abundance matrix as a tab delim file
OTU_TAB<-t(otu_table(HIRUNDO_16S.taxo.phylo))
OTU_TAB<-data.frame("#OTU ID"=taxa_names(HIRUNDO_16S.taxo.phylo),OTU_TAB)
write.table(OTU_TAB,file = "/media/kreising/DATA/data/VLASTOVKY/PICRUST2/OTU_tab.txt",row.names = F,quote = F,sep="\t")
#run picrust2 pipeline a add descriptions to resulting feature abundance matrix
#!!!must be run in terminal!!!####
picrust2_pipeline.py -s REF.fasta -i OTU_tab.txt -o picrust2_out_pipeline -p 6
add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC \
# -o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
# add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC \
# -o pathways_out/path_abun_unstrat_descrip.tsv.gz
##########################################
#Convert EC predictions to phyloseq#######
##########################################
EC<-read.delim("/media/kreising/DATA/data/VLASTOVKY/PICRUST2/picrust2_out_pipeline/EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv")
#convert descriptions to the tax_table
EC.tax<-EC[,1:2]
EC.tax<-tax_table(EC.tax)
EC.names<-EC[,1]
taxa_names(EC.tax)<-EC.names
#convert the rest of the table to the otu_table()
EC.otu<-EC[,3:dim(EC)[2]]
EC.otu<-round(EC.otu,0)
EC.otu<-otu_table(EC.otu,taxa_are_rows = T)
taxa_names(EC.otu)<-EC.names
#merge these objects with original metadata and save
HIRUNDO_picrust<-merge_phyloseq(EC.otu,EC.tax,sample_data(HIRUNDO_16S.taxo.phylo))
save(HIRUNDO_picrust,file = "/media/kreising/DATA/data/VLASTOVKY/PHYLOSEQ_rerun/HIRUNDO_picrust.R")
##########################################
#Convert pathways predictions to phyloseq#
##########################################
EC<-read.delim("/media/kreising/DATA/data/VLASTOVKY/PICRUST2/picrust2_out_pipeline/pathways_out/path_abun_unstrat_descrip.tsv")
#convert descriptions to the tax_table
EC.tax<-data.frame(Path=EC[,1:2])
EC.tax<-tax_table(EC.tax)
EC.names<-EC[,1]
taxa_names(EC.tax)<-EC.names
#convert the rest of the table to the otu_table()
EC.otu<-EC[,3:dim(EC)[2]]
EC.otu<-round(EC.otu,0)
EC.otu<-otu_table(EC.otu,taxa_are_rows = T)
taxa_names(EC.otu)<-EC.names
#merge these objects with original metadata and save
HIRUNDO_picrust_path<-merge_phyloseq(EC.otu,EC.tax,sample_data(HIRUNDO_16S.taxo.phylo))
save(HIRUNDO_picrust_path,file = "/media/kreising/DATA/data/VLASTOVKY/PHYLOSEQ_rerun/HIRUNDO_picrust_path.R")