Phyloseq + MicrobiotaProcess + PICRUSt2

gene_x 0 like s 1121 view s

Tags: R, metagenomics, 16S, tool, pipeline

  1. Phyloseq.Rmd 

    1. author: ""
    2. date: '`r format(Sys.time(), "%d %m %Y")`'
    3. header-includes:
    4.   - \usepackage{color, fancyvrb}
    5. output:
    6.   rmdformats::readthedown:
    7.     highlight: kate
    8.     number_sections : yes    
    9.   pdf_document: 
    10.     toc: yes
    11.     toc_depth: 2
    12.     number_sections : yes
    13. ---
    15. ```{r, echo=FALSE, warning=FALSE}
    16. ## Global options
    17. # TODO: reproduce the html with the additional figure/SVN-files for editing.
    18. # IMPORTANT NOTE: needs before "mkdir figures"
    19. #rmarkdown::render('Phyloseq.Rmd',output_file='Phyloseq.html')
    20. ```
    22. ```{r load-packages, include=FALSE}
    23. library(knitr)
    24. library(rmdformats)
    25. library(readxl)
    26. library(dplyr)
    27. library(kableExtra)
    29. options(max.print="75")
    30. knitr::opts_chunk$set(fig.width=8, 
    31.                       fig.height=6, 
    32.                       eval=TRUE, 
    33.                       cache=TRUE,
    34.                       echo=TRUE,
    35.                       prompt=FALSE,
    36.                       tidy=TRUE,
    37.                       comment=NA,
    38.                       message=FALSE,
    39.                       warning=FALSE)
    40. opts_knit$set(width=85)
    41. # Phyloseq R library
    42. #* Phyloseq web site : https://joey711.github.io/phyloseq/index.html
    43. #* See in particular tutorials for
    44. #    - importing data: https://joey711.github.io/phyloseq/import-data.html
    45. #    - heat maps: https://joey711.github.io/phyloseq/plot_heatmap-examples.html
    46. ```
    48. # Data
    50. Import raw data and assign sample key:
    52. ```{r, echo=TRUE, warning=FALSE}
    53. #extend map_corrected.txt with Diet and Flora
    54. #setwd("~/DATA/Data_Laura_16S_2/core_diversity_e4753")
    55. map_corrected <- read.csv("map_corrected.txt", sep="\t", row.names=1)
    56. knitr::kable(map_corrected) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    57. ```
    59. # Prerequisites to be installed
    61. * R : https://pbil.univ-lyon1.fr/CRAN/
    62. * R studio : https://www.rstudio.com/products/rstudio/download/#download
    64. ```R
    65. install.packages("dplyr")     # To manipulate dataframes
    66. install.packages("readxl")    # To read Excel files into R
    67. install.packages("ggplot2")   # for high quality graphics
    68. install.packages("heatmaply")
    69. source("https://bioconductor.org/biocLite.R")
    70. biocLite("phyloseq")
    71. ```
    73. ```{r libraries, echo=TRUE, message=FALSE}
    74. library("readxl") # necessary to import the data from Excel file
    75. library("ggplot2") # graphics
    76. library("picante")
    77. library("microbiome") # data analysis and visualisation
    78. library("phyloseq") # also the basis of data object. Data analysis and visualisation
    79. library("ggpubr") # publication quality figures, based on ggplot2
    80. library("dplyr") # data handling, filter and reformat data frames
    81. library("RColorBrewer") # nice color options
    82. library("heatmaply")
    83. library(vegan)
    84. library(gplots)
    85. ```
    87. # Read the data and create phyloseq objects
    89. Three tables are needed
    91. * OTU
    92. * Taxonomy
    93. * Samples
    95. ```{r, echo=TRUE, warning=FALSE}
    96.     #Change your working directory to where the files are located
    97.     ps.ng.tax <- import_biom("./table_even42369.biom", "./clustering/rep_set.tre")
    98.     sample <- read.csv("./map_corrected.txt", sep="\t", row.names=1)
    99.     SAM = sample_data(sample, errorIfNULL = T)
    100.     rownames(SAM) <-
    101.     c("1","2","3","5","6","7","8","9","10","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","46","47","48","49","50","51","52","53","55")
    102.     ps.ng.tax <- merge_phyloseq(ps.ng.tax, SAM)
    103.     print(ps.ng.tax) 
    104.     colnames(tax_table(ps.ng.tax)) <- c("Domain","Phylum","Class","Order","Family","Genus","Species")
    105.     saveRDS(ps.ng.tax, "./ps.ng.tax.rds")
    106. ```
    108. Visualize data
    109. ```{r, echo=TRUE, warning=FALSE}
    110.   sample_names(ps.ng.tax)
    111.   rank_names(ps.ng.tax)
    112.   sample_variables(ps.ng.tax)
    113. ```
    115. Normalize number of reads in each sample using median sequencing depth.
    116. ```{r, echo=TRUE, warning=FALSE}
    117. # RAREFACTION
    118. #set.seed(9242)  # This will help in reproducing the filtering and nomalisation. 
    119. #ps.ng.tax <- rarefy_even_depth(ps.ng.tax, sample.size = 42369)
    120. #total <- 42369
    122. # NORMALIZE number of reads in each sample using median sequencing depth.
    123. total = median(sample_sums(ps.ng.tax))
    124. #> total
    125. #[1] 42369
    126. standf = function(x, t=total) round(t * (x / sum(x)))
    127. ps.ng.tax = transform_sample_counts(ps.ng.tax, standf)
    128. ps.ng.tax_rel <- microbiome::transform(ps.ng.tax, "compositional")
    130. saveRDS(ps.ng.tax, "./ps.ng.tax.rds")
    131. hmp.meta <- meta(ps.ng.tax)
    132. hmp.meta$sam_name <- rownames(hmp.meta)
    133. ```
    135. # Heatmaps
    137. ```{r, echo=TRUE, warning=FALSE}
    138. #MOVE_FROM_ABOVE: The number of reads used for normalization is **`r sprintf("%.0f", total)`**. 
    139. #A basic heatmap using the default parameters.
    140. #  plot_heatmap(ps.ng.tax, method = "NMDS", distance = "bray")
    141. #NOTE that giving the correct OTU numbers in the text (1%, 0.5%, ...)!!!
    142. ```
    144. We consider the most abundant OTUs for heatmaps. For example one can only take OTUs that represent at least 1% of reads in at least one sample. Remember we normalized all the sampples to median number of reads (total).  We are left with only 168 OTUS which makes the reading much more easy.
    145. ```{r, echo=TRUE, warning=FALSE}
    147. # Custom function to plot a heatmap with the specified sample order
    148. #plot_heatmap_custom <- function(ps, sample_order, method = "NMDS", distance = "bray") {
    149. ps.ng.tax_abund <- phyloseq::filter_taxa(ps.ng.tax, function(x) sum(x > total*0.01) > 0, TRUE)
    150. kable(otu_table(ps.ng.tax_abund)) %>%
    151. kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    152. # Calculate the relative abundance for each sample
    153. ps.ng.tax_abund_rel <- transform_sample_counts(ps.ng.tax_abund, function(x) x / sum(x))
    155. datamat_ = as.data.frame(otu_table(ps.ng.tax_abund))
    156. #datamat <- datamat_[c("1","2","5","6","7",  "8","9","10","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","51",    "40","41","42","43","44","46",  "47","48","49","50","52","53","55")]
    157. datamat <- datamat_[c("8","9","10","12","13","14",    "21","22","23","24","25","26","27","28",    "33","34","35","36","37","38","39","51",    "47","48","49","50","52","53","55")]
    158. hr <- hclust(as.dist(1-cor(t(datamat), method="pearson")), method="complete")
    159. hc <- hclust(as.dist(1-cor(datamat, method="spearman")), method="complete")
    160. mycl = cutree(hr, h=max(hr$height)/1.08)
    161. mycol = c("YELLOW", "DARKBLUE", "DARKORANGE", "DARKMAGENTA", "DARKCYAN", "DARKRED",  "MAROON", "DARKGREEN", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN","LIGHTGREEN", "BLUE", "ORANGE", "CYAN", "RED", "GREEN");
    163. mycol = mycol[as.vector(mycl)]
    164. sampleCols <- rep('GREY',ncol(datamat))
    165. #names(sampleCols) <- c("Group1", "Group1", "Group1", "Group1", "Group1",   "Group2", "Group2", "Group2", "Group2", "Group2","Group2",    "Group3", "Group3", "Group3", "Group3",   "Group1", "Group1", "Group1", "Group1", "Group1",   "Group1", "Group1", "Group1", "Group1", "Group1",  "Group1", "Group1", "Group1", "Group1", "Group1",  "Group1", "Group1", "Group1", "Group1", "Group1",  "Group1", "Group1", "Group1", "Group1", "Group1",    "Group1", "Group1", "Group1", "Group1", "Group1",  "Group1", "Group1", "Group1", "Group1", "Group1")
    167. #sampleCols[colnames(datamat)=='1'] <- '#a6cee3'
    168. #sampleCols[colnames(datamat)=='2'] <- '#a6cee3'
    169. #sampleCols[colnames(datamat)=='5'] <- '#a6cee3'
    170. #sampleCols[colnames(datamat)=='6'] <- '#a6cee3'
    171. #sampleCols[colnames(datamat)=='7'] <- '#a6cee3'
    173. sampleCols[colnames(datamat)=='8'] <- '#1f78b4'
    174. sampleCols[colnames(datamat)=='9'] <- '#1f78b4'
    175. sampleCols[colnames(datamat)=='10'] <- '#1f78b4'
    176. sampleCols[colnames(datamat)=='12'] <- '#1f78b4'
    177. sampleCols[colnames(datamat)=='13'] <- '#1f78b4'
    178. sampleCols[colnames(datamat)=='14'] <- '#1f78b4'
    180. #sampleCols[colnames(datamat)=='15'] <- '#b2df8a'
    181. #sampleCols[colnames(datamat)=='16'] <- '#b2df8a'
    182. #sampleCols[colnames(datamat)=='17'] <- '#b2df8a'
    183. #sampleCols[colnames(datamat)=='18'] <- '#b2df8a'
    184. #sampleCols[colnames(datamat)=='19'] <- '#b2df8a'
    185. #sampleCols[colnames(datamat)=='20'] <- '#b2df8a'
    187. sampleCols[colnames(datamat)=='21'] <- '#33a02c'
    188. sampleCols[colnames(datamat)=='22'] <- '#33a02c'
    189. sampleCols[colnames(datamat)=='23'] <- '#33a02c'
    190. sampleCols[colnames(datamat)=='24'] <- '#33a02c'
    191. sampleCols[colnames(datamat)=='25'] <- '#33a02c'
    192. sampleCols[colnames(datamat)=='26'] <- '#33a02c'
    193. sampleCols[colnames(datamat)=='27'] <- '#33a02c'
    194. sampleCols[colnames(datamat)=='28'] <- '#33a02c'
    196. #sampleCols[colnames(datamat)=='29'] <- '#fb9a99'
    197. #sampleCols[colnames(datamat)=='30'] <- '#fb9a99'
    198. #sampleCols[colnames(datamat)=='31'] <- '#fb9a99'
    199. #sampleCols[colnames(datamat)=='32'] <- '#fb9a99'
    201. sampleCols[colnames(datamat)=='33'] <- '#e31a1c'
    202. sampleCols[colnames(datamat)=='34'] <- '#e31a1c'
    203. sampleCols[colnames(datamat)=='35'] <- '#e31a1c'
    204. sampleCols[colnames(datamat)=='36'] <- '#e31a1c'
    205. sampleCols[colnames(datamat)=='37'] <- '#e31a1c'
    206. sampleCols[colnames(datamat)=='38'] <- '#e31a1c'
    207. sampleCols[colnames(datamat)=='39'] <- '#e31a1c'
    208. sampleCols[colnames(datamat)=='51'] <- '#e31a1c'
    210. #sampleCols[colnames(datamat)=='40'] <- '#cab2d6'
    211. #sampleCols[colnames(datamat)=='41'] <- '#cab2d6'
    212. #sampleCols[colnames(datamat)=='42'] <- '#cab2d6'
    213. #sampleCols[colnames(datamat)=='43'] <- '#cab2d6'
    214. #sampleCols[colnames(datamat)=='44'] <- '#cab2d6'
    215. #sampleCols[colnames(datamat)=='46'] <- '#cab2d6'
    217. sampleCols[colnames(datamat)=='47'] <- '#6a3d9a'
    218. sampleCols[colnames(datamat)=='48'] <- '#6a3d9a'
    219. sampleCols[colnames(datamat)=='49'] <- '#6a3d9a'
    220. sampleCols[colnames(datamat)=='50'] <- '#6a3d9a'
    221. sampleCols[colnames(datamat)=='52'] <- '#6a3d9a'
    222. sampleCols[colnames(datamat)=='53'] <- '#6a3d9a'
    223. sampleCols[colnames(datamat)=='55'] <- '#6a3d9a'
    224. #bluered(75)
    225. #color_pattern <- colorRampPalette(c("blue", "white", "red"))(100)
    226. library(RColorBrewer)
    227. custom_palette <- colorRampPalette(brewer.pal(9, "Blues"))
    228. heatmap_colors <- custom_palette(100)
    229. #colors <- heatmap_color_default(100)
    230. png("figures/heatmap.png", width=1200, height=2400)
    231. #par(mar=c(2, 2, 2, 2))  , lwid=1    lhei=c(0.7, 10)) # Adjust height of color keys   keysize=0.3, 
    232. heatmap.2(as.matrix(datamat),Rowv=as.dendrogram(hr),Colv = NA, dendrogram = 'row',
    233.             scale='row',trace='none',col=heatmap_colors, cexRow=1.2, cexCol=1.5,
    234.             RowSideColors = mycol, ColSideColors = sampleCols, srtCol=15, labRow=row.names(datamat), key=TRUE, margins=c(10, 15), lhei=c(0.7, 15), lwid=c(1,8))
    235. dev.off()
    236. ```
    237. ```{r, echo=TRUE, warning=FALSE, fig.cap="Heatmap", out.width = '100%', fig.align= "center"}
    238. knitr::include_graphics("./figures/heatmap.png")
    239. ```
    241. ```{r, echo=FALSE, warning=FALSE}
    242.   #It is possible to use different distances and different multivaraite methods. Many different built-in distances can be used.
    243.   #dist_methods <- unlist(distanceMethodList)
    244.   #print(dist_methods)
    245. ```
    247. \pagebreak
    249. # Taxonomic summary
    251. ## Bar plots in phylum level
    253. ```{r, echo=FALSE, warning=FALSE}
    254. #Make the bargraph nicer by removing OTUs boundaries. This is done by adding ggplot2 modifier.
    255. # 1: uniform color. Color is for the border, fill is for the inside
    256. #ggplot(mtcars, aes(x=as.factor(cyl) )) +
    257. #  geom_bar(color="blue", fill=rgb(0.1,0.4,0.5,0.7) )
    258. # 2: Using Hue
    259. #ggplot(mtcars, aes(x=as.factor(cyl), fill=as.factor(cyl) )) + 
    260. #  geom_bar( ) +
    261. #  scale_fill_hue(c = 40) +
    262. #  theme(legend.position="none")
    263. # 3: Using RColorBrewer
    264. #ggplot(mtcars, aes(x=as.factor(cyl), fill=as.factor(cyl) )) + 
    265. #  geom_bar( ) +
    266. #  scale_fill_brewer(palette = "Set1") +
    267. #  theme(legend.position="none")
    268. # 4: Using greyscale:
    269. #ggplot(mtcars, aes(x=as.factor(cyl), fill=as.factor(cyl) )) + 
    270. #  geom_bar( ) +
    271. #  scale_fill_grey(start = 0.25, end = 0.75) +
    272. #  theme(legend.position="none")
    273. # 5: Set manualy
    274. #ggplot(mtcars, aes(x=as.factor(cyl), fill=as.factor(cyl) )) +  
    275. #  geom_bar( ) +
    276. #  scale_fill_manual(values = c("red", "green", "blue") ) +
    277. #  theme(legend.position="none")
    278. #NOT SUCCESSFUL!
    279. #allGroupsColors<- c(
    280. #  "grey0", "grey50", "dodgerblu", "deepskyblue",
    281. #  "red", "darkred", "green", "green4")
    282. #  plot_bar(ps.ng.tax_rel, fill="Phylum") + 
    283. #  geom_bar(stat="identity", position="stack") + scale_color_manual(values = allGroupsColors) #, fill=Phylum   + scale_fill_brewer(palette = "Set1")
    284.   # ##### Keep only the most abundant phyla and 
    285.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Actinobacteria")) #1.57
    286.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Bacteroidetes"))  #27.27436
    287.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Cyanobacteria"))  #0.02244249
    288.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Epsilonbacteraeota"))  #0.01309145
    289.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Euryarchaeota"))  #0.1210024
    290.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Firmicutes"))     #32.50589
    291.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Lentisphaerae"))  #0.0001870208
    292.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Patescibacteria")) #0.008789976
    293.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Planctomycetes")) #0.01365252
    294.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Proteobacteria")) #6.769216
    295.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Synergistetes"))  #0.005049561 
    296.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Tenericutes"))    #0.0005610623
    297.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Verrucomicrobia"))  #2.076304
    298.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c(NA))  #sum(otu_table(ps.ng.tax_most)) = 2.619413
    299. ```
    300. ```{r, echo=TRUE, warning=FALSE}
    301.   library(ggplot2)
    302.   geom.text.size = 6
    303.   theme.size = 8 #(14/5) * geom.text.size
    304.   #ps.ng.tax_most <- subset_taxa(ps.ng.tax_rel, Phylum %in% c("D_1__Actinobacteria", "D_1__Bacteroidetes", "D_1__Firmicutes", "D_1__Proteobacteria", "D_1__Verrucomicrobia", NA))
    305.   ps.ng.tax_most = phyloseq::filter_taxa(ps.ng.tax_rel, function(x) mean(x) > 0.001, TRUE)
    306.   #CONSOLE(OPTIONAL): for sampleid in 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 71 72 73; do
    307.   #echo "otu_table(ps.ng.tax_most)[,${sampleid}]=otu_table(ps.ng.tax_most)[,${sampleid}]/sum(otu_table(ps.ng.tax_most)[,${sampleid}])" done
    308.   #OR
    309.   ps.ng.tax_most_ = transform_sample_counts(ps.ng.tax_most, function(x) x / sum(x))
    310. ```
    312. ```{r, echo=FALSE, warning=FALSE}
    313.   ##--Creating 100% stacked bar plots with less abundant taxa in a sub-category #901--
    314.   ##https://github.com/joey711/phyloseq/issues/901
    315.   ##ps.ng.tax_most_df <- psmelt(ps.ng.tax_most_)  #5986x19
    316.   #glom <- tax_glom(ps.ng.tax_most_, taxrank = 'Phylum')
    317.   #tax_table(glom) # should list # taxa as # phyla
    318.   #data <- psmelt(glom) # create dataframe from phyloseq object
    319.   #data$Phylum <- as.character(data$Phylum) #convert to character
    320.   ##simple way to rename phyla with < 1% abundance
    321.   #data$Phylum[data$Abundance < 0.001] <- "< 0.1% abund."
    322.   #
    323.   #library(plyr)
    324.   #medians <- ddply(data, ~Phylum, function(x) c(median=median(x$Abundance)))
    325.   #remainder <- medians[medians$median <= 0.001,]$Phylum
    326.   #data[data$Phylum %in% remainder,]$Phylum <- "Phyla < 0.1% abund."
    327.   #data$Phylum[data$Abundance < 0.001] <- "Phyla < 0.1% abund."
    328.   ##--> data are not used!
    329.   #
    330.   ##in class level
    331.   #glom <- tax_glom(ps.ng.tax_most_, taxrank = 'Class')
    332.   #tax_table(glom) # should list # taxa as # phyla
    333.   #data <- psmelt(glom) # create dataframe from phyloseq object
    334.   #data$Class <- as.character(data$Class) #convert to character
    335.   #
    336.   ##simple way to rename phyla with < 1% abundance
    337.   #data$Class[data$Abundance < 0.001] <- "< 0.1% abund."
    338.   #Count = length(unique(data$Class))
    339.   # 
    340.   ##unique(data$Class)
    341.   ##data$Class <- factor(data$Class, levels = c("Bacilli", "Bacteroidia", "Verrucomicrobiae", "Clostridia", "Gammaproteobacteria", "Alphaproteobacteria", "Actinobacteria", "Negativicutes", "Erysipelotrichia", "Methanobacteria", "< 0.1% abund."))
    342.   ##------- Creating 100% stacked bar plots END --------
    344.   library(stringr)
    345. #FITTING1: 
    346. # tax_table(ps.ng.tax_most_)[1,"Domain"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Domain"], "__")[[1]][2]
    347. # ... ...
    348. # tax_table(ps.ng.tax_most_)[167,"Species"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Species"], "__")[[1]][2]
    349. #ps.ng.tax_most_
    350. #in total [ 89 taxa and 55 samples ]
    351. #otu_table()   OTU Table:         [ 166 taxa and 54 samples ]
    352. #otu_table()   OTU Table:         [ 168 taxa and 50 samples ]
    353. #for id in 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 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100  101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166; do   
    354. #echo "tax_table(ps.ng.tax_most_)[${id},\"Domain\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Domain\"], \"__\")[[1]][2]"
    355. #echo "tax_table(ps.ng.tax_most_)[${id},\"Phylum\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Phylum\"], \"__\")[[1]][2]"
    356. #echo "tax_table(ps.ng.tax_most_)[${id},\"Class\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Class\"], \"__\")[[1]][2]"
    357. #echo "tax_table(ps.ng.tax_most_)[${id},\"Order\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Order\"], \"__\")[[1]][2]"
    358. #echo "tax_table(ps.ng.tax_most_)[${id},\"Family\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Family\"], \"__\")[[1]][2]"
    359. #echo "tax_table(ps.ng.tax_most_)[${id},\"Genus\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Genus\"], \"__\")[[1]][2]"
    360. #echo "tax_table(ps.ng.tax_most_)[${id},\"Species\"] <- str_split(tax_table(ps.ng.tax_most_)[${id},\"Species\"], \"__\")[[1]][2]"
    361. #done
    363. tax_table(ps.ng.tax_most_)[1,"Domain"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Domain"], "__")[[1]][2]
    364. tax_table(ps.ng.tax_most_)[1,"Phylum"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Phylum"], "__")[[1]][2]
    365. tax_table(ps.ng.tax_most_)[1,"Class"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Class"], "__")[[1]][2]
    366. tax_table(ps.ng.tax_most_)[1,"Order"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Order"], "__")[[1]][2]
    367. tax_table(ps.ng.tax_most_)[1,"Family"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Family"], "__")[[1]][2]
    368. tax_table(ps.ng.tax_most_)[1,"Genus"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Genus"], "__")[[1]][2]
    369. tax_table(ps.ng.tax_most_)[1,"Species"] <- str_split(tax_table(ps.ng.tax_most_)[1,"Species"], "__")[[1]][2]
    370. #... ...
    371. tax_table(ps.ng.tax_most_)[167,"Domain"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Domain"], "__")[[1]][2]
    372. tax_table(ps.ng.tax_most_)[167,"Phylum"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Phylum"], "__")[[1]][2]
    373. tax_table(ps.ng.tax_most_)[167,"Class"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Class"], "__")[[1]][2]
    374. tax_table(ps.ng.tax_most_)[167,"Order"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Order"], "__")[[1]][2]
    375. tax_table(ps.ng.tax_most_)[167,"Family"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Family"], "__")[[1]][2]
    376. tax_table(ps.ng.tax_most_)[167,"Genus"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Genus"], "__")[[1]][2]
    377. tax_table(ps.ng.tax_most_)[167,"Species"] <- str_split(tax_table(ps.ng.tax_most_)[167,"Species"], "__")[[1]][2]
    378. ```
    380. ```{r, echo=TRUE, warning=FALSE}
    381.   #aes(color="Phylum", fill="Phylum") --> aes()
    382.   #ggplot(data=data, aes(x=Sample, y=Abundance, fill=Phylum)) 
    383.   plot_bar(ps.ng.tax_most_, fill="Phylum") + geom_bar(aes(), stat="identity", position="stack") +
    384.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2))                                  #6 instead of theme.size
    385. ```
    386. ```{r, echo=FALSE, warning=FALSE}
    387.   #png("abc.png")
    388.   #knitr::include_graphics("./Phyloseq_files/figure-html/unnamed-chunk-7-1.png")
    389.   #dev.off()
    390. ```
    392. \pagebreak
    394. Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.
    396. ```{r, echo=TRUE, warning=FALSE}
    397.   ps.ng.tax_most_pre_post_stroke <- merge_samples(ps.ng.tax_most_, "pre_post_stroke")
    398.   ps.ng.tax_most_pre_post_stroke_ = transform_sample_counts(ps.ng.tax_most_pre_post_stroke, function(x) x / sum(x))
    399.   #plot_bar(ps.ng.tax_most_SampleType_, fill = "Phylum") + geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack")
    400.   plot_bar(ps.ng.tax_most_pre_post_stroke_, fill="Phylum") + geom_bar(aes(), stat="identity", position="stack") +
    401.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = theme.size, colour="black"))
    402. ```
    404. \pagebreak
    406. Use color according to phylum. Do separate panels Stroke and Sex_age.
    407. ```{r, echo=TRUE, warning=FALSE}  
    408.   ps.ng.tax_most_copied <- data.table::copy(ps.ng.tax_most_)
    409.   #FITTING6: regulate the bar height if it has replicates: 5+6+6+8+4+8+6+7=25+25=50
    410.   otu_table(ps.ng.tax_most_)[,c("1")] <- otu_table(ps.ng.tax_most_)[,c("1")]/5
    411.   otu_table(ps.ng.tax_most_)[,c("2")] <- otu_table(ps.ng.tax_most_)[,c("2")]/5
    412.   otu_table(ps.ng.tax_most_)[,c("5")] <- otu_table(ps.ng.tax_most_)[,c("5")]/5
    413.   otu_table(ps.ng.tax_most_)[,c("6")] <- otu_table(ps.ng.tax_most_)[,c("6")]/5
    414.   otu_table(ps.ng.tax_most_)[,c("7")] <- otu_table(ps.ng.tax_most_)[,c("7")]/5
    416.   otu_table(ps.ng.tax_most_)[,c("8")] <- otu_table(ps.ng.tax_most_)[,c("8")]/6
    417.   otu_table(ps.ng.tax_most_)[,c("9")] <- otu_table(ps.ng.tax_most_)[,c("9")]/6
    418.   otu_table(ps.ng.tax_most_)[,c("10")] <- otu_table(ps.ng.tax_most_)[,c("10")]/6
    419.   otu_table(ps.ng.tax_most_)[,c("12")] <- otu_table(ps.ng.tax_most_)[,c("12")]/6
    420.   otu_table(ps.ng.tax_most_)[,c("13")] <- otu_table(ps.ng.tax_most_)[,c("13")]/6
    421.   otu_table(ps.ng.tax_most_)[,c("14")] <- otu_table(ps.ng.tax_most_)[,c("14")]/6
    423.   otu_table(ps.ng.tax_most_)[,c("15")] <- otu_table(ps.ng.tax_most_)[,c("15")]/6
    424.   otu_table(ps.ng.tax_most_)[,c("16")] <- otu_table(ps.ng.tax_most_)[,c("16")]/6
    425.   otu_table(ps.ng.tax_most_)[,c("17")] <- otu_table(ps.ng.tax_most_)[,c("17")]/6
    426.   otu_table(ps.ng.tax_most_)[,c("18")] <- otu_table(ps.ng.tax_most_)[,c("18")]/6
    427.   otu_table(ps.ng.tax_most_)[,c("19")] <- otu_table(ps.ng.tax_most_)[,c("19")]/6
    428.   otu_table(ps.ng.tax_most_)[,c("20")] <- otu_table(ps.ng.tax_most_)[,c("20")]/6
    430.   otu_table(ps.ng.tax_most_)[,c("21")] <- otu_table(ps.ng.tax_most_)[,c("21")]/8
    431.   otu_table(ps.ng.tax_most_)[,c("22")] <- otu_table(ps.ng.tax_most_)[,c("22")]/8
    432.   otu_table(ps.ng.tax_most_)[,c("23")] <- otu_table(ps.ng.tax_most_)[,c("23")]/8
    433.   otu_table(ps.ng.tax_most_)[,c("24")] <- otu_table(ps.ng.tax_most_)[,c("24")]/8
    434.   otu_table(ps.ng.tax_most_)[,c("25")] <- otu_table(ps.ng.tax_most_)[,c("25")]/8
    435.   otu_table(ps.ng.tax_most_)[,c("26")] <- otu_table(ps.ng.tax_most_)[,c("26")]/8
    436.   otu_table(ps.ng.tax_most_)[,c("27")] <- otu_table(ps.ng.tax_most_)[,c("27")]/8
    437.   otu_table(ps.ng.tax_most_)[,c("28")] <- otu_table(ps.ng.tax_most_)[,c("28")]/8
    439.   otu_table(ps.ng.tax_most_)[,c("29")] <- otu_table(ps.ng.tax_most_)[,c("29")]/4
    440.   otu_table(ps.ng.tax_most_)[,c("30")] <- otu_table(ps.ng.tax_most_)[,c("30")]/4
    441.   otu_table(ps.ng.tax_most_)[,c("31")] <- otu_table(ps.ng.tax_most_)[,c("31")]/4
    442.   otu_table(ps.ng.tax_most_)[,c("32")] <- otu_table(ps.ng.tax_most_)[,c("32")]/4
    444.   otu_table(ps.ng.tax_most_)[,c("33")] <- otu_table(ps.ng.tax_most_)[,c("33")]/8
    445.   otu_table(ps.ng.tax_most_)[,c("34")] <- otu_table(ps.ng.tax_most_)[,c("34")]/8
    446.   otu_table(ps.ng.tax_most_)[,c("35")] <- otu_table(ps.ng.tax_most_)[,c("35")]/8
    447.   otu_table(ps.ng.tax_most_)[,c("36")] <- otu_table(ps.ng.tax_most_)[,c("36")]/8
    448.   otu_table(ps.ng.tax_most_)[,c("37")] <- otu_table(ps.ng.tax_most_)[,c("37")]/8
    449.   otu_table(ps.ng.tax_most_)[,c("38")] <- otu_table(ps.ng.tax_most_)[,c("38")]/8
    450.   otu_table(ps.ng.tax_most_)[,c("39")] <- otu_table(ps.ng.tax_most_)[,c("39")]/8
    451.   otu_table(ps.ng.tax_most_)[,c("51")] <- otu_table(ps.ng.tax_most_)[,c("51")]/8
    453.   otu_table(ps.ng.tax_most_)[,c("40")] <- otu_table(ps.ng.tax_most_)[,c("40")]/6
    454.   otu_table(ps.ng.tax_most_)[,c("41")] <- otu_table(ps.ng.tax_most_)[,c("41")]/6
    455.   otu_table(ps.ng.tax_most_)[,c("42")] <- otu_table(ps.ng.tax_most_)[,c("42")]/6
    456.   otu_table(ps.ng.tax_most_)[,c("43")] <- otu_table(ps.ng.tax_most_)[,c("43")]/6
    457.   otu_table(ps.ng.tax_most_)[,c("44")] <- otu_table(ps.ng.tax_most_)[,c("44")]/6
    458.   otu_table(ps.ng.tax_most_)[,c("46")] <- otu_table(ps.ng.tax_most_)[,c("46")]/6
    460.   otu_table(ps.ng.tax_most_)[,c("47")] <- otu_table(ps.ng.tax_most_)[,c("47")]/7
    461.   otu_table(ps.ng.tax_most_)[,c("48")] <- otu_table(ps.ng.tax_most_)[,c("48")]/7
    462.   otu_table(ps.ng.tax_most_)[,c("49")] <- otu_table(ps.ng.tax_most_)[,c("49")]/7
    463.   otu_table(ps.ng.tax_most_)[,c("50")] <- otu_table(ps.ng.tax_most_)[,c("50")]/7
    464.   otu_table(ps.ng.tax_most_)[,c("52")] <- otu_table(ps.ng.tax_most_)[,c("52")]/7
    465.   otu_table(ps.ng.tax_most_)[,c("53")] <- otu_table(ps.ng.tax_most_)[,c("53")]/7
    466.   otu_table(ps.ng.tax_most_)[,c("55")] <- otu_table(ps.ng.tax_most_)[,c("55")]/7
    468.   #plot_bar(ps.ng.tax_most_swab_, x="Phylum", fill = "Phylum", facet_grid = Patient~RoundDay) + geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") + theme(axis.text = element_text(size = theme.size, colour="black"))
    469.   plot_bar(ps.ng.tax_most_, x="Phylum", fill="Phylum", facet_grid = pre_post_stroke~Sex_age) + geom_bar(aes(), stat="identity", position="stack") +
    470.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black"), axis.text.x=element_blank(), axis.ticks=element_blank()) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2))
    471. ```
    472. ```{r, echo=FALSE, warning=FALSE}  
    473.   #knitr::include_graphics("./Phyloseq_files/figure-html/unnamed-chunk-10-1.png")
    474.   #> tax_table(carbom)
    475.   #Taxonomy Table:     [205 taxa by 7 taxonomic ranks]:
    476.   #     Domain      Supergroup       Division      Class                 
    477.   #Otu001 "Eukaryota" "Archaeplastida" "Chlorophyta" "Mamiellophyceae"    
    478.   #       Order                      Family                 Genus     
    479.   #Otu001 "Mamiellales"              "Bathycoccaceae"       "Ostreococcus"
    480.   #sample_data(ps.ng.tax)
    481. ```
    483. ## Bar plots in class level
    484. ```{r, echo=TRUE, warning=FALSE}
    485.   plot_bar(ps.ng.tax_most_copied, fill="Class") + geom_bar(aes(), stat="identity", position="stack") +
    486.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=3))
    487. ```
    489. Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.
    490. ```{r, echo=TRUE, warning=FALSE}
    491.   plot_bar(ps.ng.tax_most_pre_post_stroke_, fill="Class") + geom_bar(aes(), stat="identity", position="stack") +
    492.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = theme.size, colour="black"))
    493. ```
    494. \pagebreak
    496. Use color according to class. Do separate panels Stroke and Sex_age.
    497. ```{r, echo=TRUE, warning=FALSE}
    498.   #-- If existing replicates, to be processed as follows --
    499.   plot_bar(ps.ng.tax_most_, x="Class", fill="Class", facet_grid = pre_post_stroke~Sex_age) + geom_bar(aes(), stat="identity", position="stack") +
    500.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black"), axis.text.x=element_blank(), axis.ticks=element_blank()) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=3)) 
    501. ```
    503. ## Bar plots in order level
    505. ```{r, echo=TRUE, warning=FALSE}
    506.   plot_bar(ps.ng.tax_most_copied, fill="Order") + geom_bar(aes(), stat="identity", position="stack") +
    507.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=4))                                 
    508. ```
    510. Regroup together pre vs post stroke and normalize number of reads in each group using median sequencing depth.
    511. ```{r, echo=TRUE, warning=FALSE}
    512.   plot_bar(ps.ng.tax_most_pre_post_stroke_, fill="Order") + geom_bar(aes(), stat="identity", position="stack") +
    513.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = theme.size, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=4))
    514. ```
    515. \pagebreak
    517. Use color according to order. Do separate panels Stroke and Sex_age.
    518. ```{r, echo=TRUE, warning=FALSE}
    519.   #FITTING7: regulate the bar height if it has replicates
    520.   plot_bar(ps.ng.tax_most_, x="Order", fill="Order", facet_grid = pre_post_stroke~Sex_age) + geom_bar(aes(), stat="identity", position="stack") +
    521.   scale_fill_manual(values = c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")) + theme(axis.text = element_text(size = 5, colour="black"), axis.text.x=element_blank(), axis.ticks=element_blank()) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=4))
    522. ```
    524. ## Bar plots in family level
    526. ```{r, echo=TRUE, warning=FALSE}
    527.   plot_bar(ps.ng.tax_most_copied, fill="Family") + geom_bar(aes(), stat="identity", position="stack") +
    528.   scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = 5, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=8))
    529. ```
    531. Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.
    532. ```{r, echo=TRUE, warning=FALSE}
    533.   plot_bar(ps.ng.tax_most_pre_post_stroke_, fill="Family") + geom_bar(aes(), stat="identity", position="stack") +
    534.   scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = theme.size, colour="black")) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=8))
    535. ```
    536. \pagebreak
    538. Use color according to family. Do separate panels Stroke and Sex_age.
    539. ```{r, echo=TRUE, warning=FALSE}
    540.   #-- If existing replicates, to be processed as follows --
    541.   plot_bar(ps.ng.tax_most_, x="Family", fill="Family", facet_grid = pre_post_stroke~Sex_age) + geom_bar(aes(), stat="identity", position="stack") +
    542.   scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = 5, colour="black"), axis.text.x=element_blank(), axis.ticks=element_blank()) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=8)) 
    543. ```
    545. ```{r, echo=FALSE, warning=FALSE}
    546.   #MOVE_FROM_ABOVE: ## Bar plots in genus level
    547.   #MOVE_FROM_ABOVE: Regroup together pre vs post stroke samples and normalize number of reads in each group using median sequencing depth.
    548.   #plot_bar(ps.ng.tax_most_pre_post_stroke_, fill="Genus") + geom_bar(aes(), stat="identity", position="stack") +
    549.   #scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = theme.size, colour="black")) + theme(legend.position="bottom")
    550. ```
    551. \pagebreak
    553. ```{r, echo=FALSE, warning=FALSE}
    554.   #MOVE_FROM_ABOVE: Use color according to genus. Do separate panels Stroke and Sex_age.
    555.   ##-- If existing replicates, to be processed as follows --
    556.   #plot_bar(ps.ng.tax_most_, x="Genus", fill="Genus", facet_grid = pre_post_stroke~Sex_age) + geom_bar(aes(), stat="identity", position="stack") +
    557.   #scale_fill_manual(values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7")) + theme(axis.text = element_text(size = 6, colour="black"), axis.text.x=element_blank(), axis.ticks=element_blank()) + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=18)) 
    558. ```
    560. \pagebreak
    562. # Alpha diversity
    563. Plot Chao1 richness estimator, Observed OTUs, Shannon index, and Phylogenetic diversity. 
    564. Regroup together samples from the same group.
    565. ```{r, echo=FALSE, warning=FALSE}
    566. # using rarefied data
    567. #FITTING2: CONSOLE: 
    568. #gunzip table_even4753.biom.gz
    569. #alpha_diversity.py -i table_even42369.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering/rep_set.tre
    570. #gunzip table_even4753.biom.gz
    571. #alpha_diversity.py -i table_even4753.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering_stool/rep_set.tre
    572. #gunzip table_even4753.biom.gz
    573. #alpha_diversity.py -i table_even4753.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering_swab/rep_set.tre
    574. ```
    576. ```{r, echo=TRUE, warning=FALSE}
    577. hmp.div_qiime <- read.csv("adiv_even.txt", sep="\t") 
    578. colnames(hmp.div_qiime) <- c("sam_name", "chao1", "observed_otus", "shannon", "PD_whole_tree")
    579. row.names(hmp.div_qiime) <- hmp.div_qiime$sam_name
    580. div.df <- merge(hmp.div_qiime, hmp.meta, by = "sam_name")
    581. div.df2 <- div.df[, c("Group", "chao1", "shannon", "observed_otus", "PD_whole_tree")]
    582. colnames(div.df2) <- c("Group", "Chao-1", "Shannon", "OTU", "Phylogenetic Diversity")
    583. #colnames(div.df2)
    584. options(max.print=999999)
    585. #27     H47 830.5000 5.008482 319               10.60177
    586. #FITTING4: if occuring "Computation failed in `stat_signif()`:not enough 'y' observations"
    587. #means: the patient H47 contains only one sample, it should be removed for the statistical p-values calculations. 
    588. #delete H47(1)
    589. #div.df2 <- div.df2[-c(3), ] 
    590. #div.df2 <- div.df2[-c(55,54, 45,40,39,27,26,25,1), ] 
    591. knitr::kable(div.df2) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    593. #https://uc-r.github.io/t_test
    594. #We can perform the test with t.test and transform our data and we can also perform the nonparametric test with the wilcox.test function.
    595. stat.test.Shannon <- compare_means(
    596. Shannon ~ Group, data = div.df2,
    597. method = "t.test"
    598. )
    599. knitr::kable(stat.test.Shannon) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    601. div_df_melt <- reshape2::melt(div.df2)
    602. #head(div_df_melt)
    604. #https://plot.ly/r/box-plots/#horizontal-boxplot
    605. #http://www.sthda.com/english/wiki/print.php?id=177
    606. #https://rpkgs.datanovia.com/ggpubr/reference/as_ggplot.html
    607. #http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/82-ggplot2-easy-way-to-change-graphical-parameters/
    608. #https://plot.ly/r/box-plots/#horizontal-boxplot
    609. #library("gridExtra")
    610. #par(mfrow=c(4,1))
    611. p <- ggboxplot(div_df_melt, x = "Group", y = "value",
    612.               facet.by = "variable", 
    613.               scales = "free",
    614.               width = 0.5,
    615.               fill = "gray", legend= "right")
    616. #ggpar(p, xlab = FALSE, ylab = FALSE)
    617. lev <- levels(factor(div_df_melt$Group)) # get the variables
    618. #FITTING4: delete H47(1) in lev
    619. #lev <- lev[-c(3)]
    620. # make a pairwise list that we want to compare.
    621. #my_stat_compare_means
    622. #https://stackoverflow.com/questions/47839988/indicating-significance-with-ggplot2-in-a-boxplot-with-multiple-groups
    623. L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i]) #%>% filter(p.signif != "ns")
    624. my_stat_compare_means  <- function (mapping = NULL, data = NULL, method = NULL, paired = FALSE, 
    625.     method.args = list(), ref.group = NULL, comparisons = NULL, 
    626.     hide.ns = FALSE, label.sep = ", ", label = NULL, label.x.npc = "left", 
    627.     label.y.npc = "top", label.x = NULL, label.y = NULL, tip.length = 0.03, 
    628.     symnum.args = list(), geom = "text", position = "identity", 
    629.     na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, ...) 
    630. {
    631.     if (!is.null(comparisons)) {
    632.         method.info <- ggpubr:::.method_info(method)
    633.         method <- method.info$method
    634.         method.args <- ggpubr:::.add_item(method.args, paired = paired)
    635.         if (method == "wilcox.test") 
    636.             method.args$exact <- FALSE
    637.         pms <- list(...)
    638.         size <- ifelse(is.null(pms$size), 0.3, pms$size)
    639.         color <- ifelse(is.null(pms$color), "black", pms$color)
    640.         map_signif_level <- FALSE
    641.         if (is.null(label)) 
    642.             label <- "p.format"
    643.         if (ggpubr:::.is_p.signif_in_mapping(mapping) | (label %in% "p.signif")) {
    644.             if (ggpubr:::.is_empty(symnum.args)) {
    645.                 map_signif_level <- c(`****` = 1e-04, `***` = 0.001, 
    646.                   `**` = 0.01, `*` = 0.05, ns = 1)
    647.             } else {
    648.               map_signif_level <- symnum.args
    649.             } 
    650.             if (hide.ns) 
    651.                 names(map_signif_level)[5] <- " "
    652.         }
    653.         step_increase <- ifelse(is.null(label.y), 0.12, 0)
    654.         ggsignif::geom_signif(comparisons = comparisons, y_position = label.y, 
    655.             test = method, test.args = method.args, step_increase = step_increase, 
    656.             size = size, color = color, map_signif_level = map_signif_level, 
    657.             tip_length = tip.length, data = data)
    658.     } else {
    659.         mapping <- ggpubr:::.update_mapping(mapping, label)
    660.         layer(stat = StatCompareMeans, data = data, mapping = mapping, 
    661.             geom = geom, position = position, show.legend = show.legend, 
    662.             inherit.aes = inherit.aes, params = list(label.x.npc = label.x.npc, 
    663.                 label.y.npc = label.y.npc, label.x = label.x, 
    664.                 label.y = label.y, label.sep = label.sep, method = method, 
    665.                 method.args = method.args, paired = paired, ref.group = ref.group, 
    666.                 symnum.args = symnum.args, hide.ns = hide.ns, 
    667.                 na.rm = na.rm, ...))
    668.     }
    669. }
    671. p2 <- p + 
    672. stat_compare_means(
    673.   method="t.test",
    674.   comparisons = list(c("Group1", "Group2"), c("Group1", "Group3"), c("Group1", "Group4"), c("Group1", "Group6"), c("Group1", "Group8"), c("Group2", "Group5"),c("Group4", "Group5"),c("Group4", "Group6"),c("Group4", "Group7"),c("Group6", "Group7")), 
    675.   label = "p.signif",
    676.   symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns"))
    677. )
    678. #comparisons = L.pairs,
    679. #symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
    680. #stat_pvalue_manual
    681. #print(p2)
    682. #https://stackoverflow.com/questions/20500706/saving-multiple-ggplots-from-ls-into-one-and-separate-files-in-r
    683. #FITTING3: mkdir figures
    684. ggsave("./figures/alpha_diversity_Group.png", device="png", height = 10, width = 12)
    685. ggsave("./figures/alpha_diversity_Group.svg", device="svg", height = 10, width = 12)
    687. p3 <- p + 
    688. stat_compare_means(
    689.   method="t.test",
    690.   comparisons = list(c("Group2", "Group4"), c("Group2", "Group6"), c("Group4", "Group8"), c("Group6", "Group8")), 
    691.   label = "p.signif",
    692.   symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
    693. )
    694. #symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
    695. #stat_pvalue_manual
    696. #print(p2)
    697. #https://stackoverflow.com/questions/20500706/saving-multiple-ggplots-from-ls-into-one-and-separate-files-in-r
    698. #FITTING3: mkdir figures
    699. ggsave("./figures/alpha_diversity_Group2.png", device="png", height = 10, width = 12)
    700. ggsave("./figures/alpha_diversity_Group2.svg", device="svg", height = 10, width = 12)
    701. ```
    703. # Selected alpha diversity
    704. ```{r, echo=TRUE, warning=FALSE, fig.cap="Alpha diversity", out.width = '100%', fig.align= "center"}
    705. knitr::include_graphics("./figures/alpha_diversity_Group2.png")
    707. selected_alpha_diversities<-read.csv("selected_alpha_diversities.txt",sep="\t")
    708. knitr::kable(selected_alpha_diversities) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    709. ```
    711. # Beta diversity
    712. ```{r, echo=TRUE, warning=FALSE, fig.cap="Beta diversity", out.width = '100%', fig.align= "center"}
    713. #file:///home/jhuang/DATA/Data_Marius_16S/core_diversity_e42369/bdiv_even42369_Group/unweighted_unifrac_boxplots/Group_Stats.txt
    714. beta_diversity_group_stats<-read.csv("unweighted_unifrac_boxplots_Group_Stats.txt",sep="\t")
    715. knitr::kable(beta_diversity_group_stats) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    716. #NOTE: Run this Phyloseq0.Rmd, then run the code of MicrobiotaProcess.R to manually generate PCoA.png, then run this Phyloseq.Rmd!
    717. #NOTE: AT_FIRST_DEACTIVATE_THIS_LINE: knitr::include_graphics("./figures/PCoA.png")
    718. ```
    720. # Differential abundance analysis
    722. Differential abundance analysis aims to find the differences in the abundance of each taxa between two groups of samples, assigning a significance value to each comparison.
    724. ## Group2 vs Group4
    726. ```{r, echo=TRUE, warning=FALSE}
    727. library("DESeq2")
    728. #ALTERNATIVE using ps.ng.tax_most_copied: ps.ng.tax (40594) vs. ps.ng.tax_most_copied (166)
    729. ps.ng.tax_sel <- ps.ng.tax
    730. #FITTING5: correct the id of the group members, see FITTING6
    731. otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("8","9","10","12","13","14",  "21","22","23","24","25","26","27","28")]
    732. diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
    733. diagdds$Group <- relevel(diagdds$Group, "Group4")
    734. diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
    735. resultsNames(diagdds)
    737. res = results(diagdds, cooksCutoff = FALSE)
    738. alpha = 0.05
    739. sigtab = res[which(res$padj < alpha), ]
    740. sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab), ], "matrix"))
    741. sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)), ]
    742. kable(sigtab) %>%
    743.   kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    744. #rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied))
    745. #subv %in% v
    746. ### returns a vector TRUE FALSE
    747. #is.element(subv, v)
    748. ### returns a vector TRUE FALSE
    750. library("ggplot2")
    751. theme_set(theme_bw())
    752. scale_fill_discrete <- function(palname = "Set1", ...) {
    753.     scale_fill_brewer(palette = palname, ...)
    754. }
    755. x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
    756. x = sort(x)
    757. sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
    758. x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
    759. x = sort(x)
    760. sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))
    761. ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    762.   theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
    764. #Error in checkForExperimentalReplicates(object, modelMatrix) : 
    765. #  The design matrix has the same number of samples and coefficients to fit,
    766. #  so estimation of dispersion is not possible. Treating samples
    767. #  as replicates was deprecated in v1.20 and no longer supported since v1.22.
    768. ```
    770. ## Group2 vs Group6
    772. ```{r, echo=TRUE, warning=FALSE}
    773. ps.ng.tax_sel <- ps.ng.tax
    774. otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("8","9","10","12","13","14",    "33","34","35","36","37","38","39","51")]
    775. diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
    776. diagdds$Group <- relevel(diagdds$Group, "Group6")
    777. diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
    778. resultsNames(diagdds)
    780. res = results(diagdds, cooksCutoff = FALSE)
    781. alpha = 0.05
    782. sigtab = res[which(res$padj < alpha), ]
    783. sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab), ], "matrix"))
    784. sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)), ]
    785. kable(sigtab) %>%
    786.   kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    788. library("ggplot2")
    789. theme_set(theme_bw())
    790. scale_fill_discrete <- function(palname = "Set1", ...) {
    791.     scale_fill_brewer(palette = palname, ...)
    792. }
    793. x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
    794. x = sort(x)
    795. sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
    796. x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
    797. x = sort(x)
    798. sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))
    799. ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    800.   theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
    801. ```
    803. ## Group4 vs Group8
    805. ```{r, echo=TRUE, warning=FALSE}
    806. ps.ng.tax_sel <- ps.ng.tax
    807. otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("21","22","23","24","25","26","27","28",    "47","48","49","50","52","53","55")]
    808. diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
    809. diagdds$Group <- relevel(diagdds$Group, "Group8")
    810. diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
    811. resultsNames(diagdds)
    813. res = results(diagdds, cooksCutoff = FALSE)
    814. alpha = 0.05
    815. sigtab = res[which(res$padj < alpha), ]
    816. sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab), ], "matrix"))
    817. sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)), ]
    818. kable(sigtab) %>%
    819.   kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    821. library("ggplot2")
    822. theme_set(theme_bw())
    823. scale_fill_discrete <- function(palname = "Set1", ...) {
    824.     scale_fill_brewer(palette = palname, ...)
    825. }
    826. x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
    827. x = sort(x)
    828. sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
    829. x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
    830. x = sort(x)
    831. sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))
    832. ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    833.   theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
    834. ```
    836. ## Group6 vs Group8
    838. ```{r, echo=TRUE, warning=FALSE}
    839. ps.ng.tax_sel <- ps.ng.tax
    840. otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("33","34","35","36","37","38","39","51",    "47","48","49","50","52","53","55")]
    841. diagdds = phyloseq_to_deseq2(ps.ng.tax_sel, ~Group)
    842. diagdds$Group <- relevel(diagdds$Group, "Group8")
    843. diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
    844. resultsNames(diagdds)
    846. res = results(diagdds, cooksCutoff = FALSE)
    847. alpha = 0.05
    848. sigtab = res[which(res$padj < alpha), ]
    849. sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(ps.ng.tax_sel)[rownames(sigtab), ], "matrix"))
    850. sigtab <- sigtab[rownames(sigtab) %in% rownames(tax_table(ps.ng.tax_most_copied)), ]
    851. kable(sigtab) %>%
    852.   kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    854. library("ggplot2")
    855. theme_set(theme_bw())
    856. scale_fill_discrete <- function(palname = "Set1", ...) {
    857.     scale_fill_brewer(palette = palname, ...)
    858. }
    859. x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
    860. x = sort(x)
    861. sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
    862. x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
    863. x = sort(x)
    864. sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))
    865. ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    866.   theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
    867. ```

  2. MicrobiotaProcess.R

    1. # -----------------------------------
    2. # ---- prepare the R environment ----
    3. #Rscript MicrobiotaProcess.R
    4. #NOTE: exit R script, then login again R-environment; rm -rf Phyloseq*_cache
    5. rmarkdown::render('Phyloseq.Rmd',output_file='Phyloseq.html')
    7. # -----------------------------
    8. # ---- 3.1. bridges other tools
    9. ##https://github.com/YuLab-SMU/MicrobiotaProcess
    10. ##https://www.bioconductor.org/packages/release/bioc/vignettes/MicrobiotaProcess/inst/doc/MicrobiotaProcess.html
    11. ##https://chiliubio.github.io/microeco_tutorial/intro.html#framework
    12. ##https://yiluheihei.github.io/microbiomeMarker/reference/plot_cladogram.html
    13. #BiocManager::install("MicrobiotaProcess")
    14. #install.packages("microeco")
    15. #install.packages("ggalluvial")
    16. #install.packages("ggh4x")
    18. library(MicrobiotaProcess)    
    19. library(microeco)
    20. library(ggalluvial)
    21. library(ggh4x)
    22. library(gghalves)
    24. ## Convert the phyloseq object to a MicrobiotaProcess object
    25. #mp <- as.MicrobiotaProcess(ps.ng.tax)
    27. #mt <- phyloseq2microeco(ps.ng.tax) #--> ERROR
    28. #abundance_table <- mt$abun_table
    29. #taxonomy_table <- mt$tax_table
    31. #ps.ng.tax_abund <- phyloseq::filter_taxa(ps.ng.tax, function(x) sum(x > total*0.01) > 0, TRUE)
    32. #ps.ng.tax_most = phyloseq::filter_taxa(ps.ng.tax_rel, function(x) mean(x) > 0.001, TRUE)
    34. ##OPTION1: take all samples, prepare mpse_abund!
    35. ##mpse <- ps.ng.tax %>% as.MPSE()
    36. #mpse_abund <- ps.ng.tax_abund %>% as.MPSE()
    38. ##OPTION2: take partial samples, prepare mpse_abund
    39. ps.ng.tax_sel <- ps.ng.tax  #IMPORTANT
    40. ##otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("1","2","5","6","7",  "15","16","17","18","19","20",  "29","30","31","32",  "40","41","42","43","44","46")]
    41. ##NOTE: Only choose Group2, Group4, Group6, Group8
    42. #> ps.ng.tax_sel
    43. #otu_table()   OTU Table:         [ 37465 taxa and 29 samples ]
    44. #sample_data() Sample Data:       [ 29 samples by 10 sample variables ]
    45. #tax_table()   Taxonomy Table:    [ 37465 taxa by 7 taxonomic ranks ]
    46. #phy_tree()    Phylogenetic Tree: [ 37465 tips and 37461 internal nodes ]
    47. otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax)[,c("8","9","10","12","13","14",  "21","22","23","24","25","26","27","28",  "33","34","35","36","37","38","39","51",  "47","48","49","50","52","53","55")]
    48. #For quick calculation
    49. #otu_table(ps.ng.tax_sel) <- otu_table(ps.ng.tax_abund)[,c("8","9","10","12","13","14",  "21","22","23","24","25","26","27","28",  "33","34","35","36","37","38","39","51",  "47","48","49","50","52","53","55")]
    50. mpse_abund <- ps.ng.tax_sel %>% as.MPSE()
    52. # -----------------------------------
    53. # ---- 3.2. alpha diversity analysis 
    54. # Rarefied species richness + RareAbundance
    55. mpse_abund %<>% mp_rrarefy()
    56. # 'chunks' represent the split number of each sample to calculate alpha
    57. # diversity, default is 400. e.g. If a sample has total 40000
    58. # reads, if chunks is 400, it will be split to 100 sub-samples
    59. # (100, 200, 300,..., 40000), then alpha diversity index was
    60. # calculated based on the sub-samples. 
    61. # '.abundance' the column name of abundance, if the '.abundance' is not be 
    62. # rarefied calculate rarecurve, user can specific 'force=TRUE'.
    64. # + RareAbundance
    65. mpse_abund %<>% 
    66.     mp_cal_rarecurve(
    67.         .abundance = RareAbundance,
    68.         chunks = 400
    69.     )
    70. # The RareAbundanceRarecurve column will be added the colData slot 
    71. # automatically (default action="add")
    72. mpse_abund %>% print(width=180, n=100)
    74. # default will display the confidence interval around smooth.
    75. # se=TRUE
    76. p1 <- mpse_abund %>% 
    77.       mp_plot_rarecurve(
    78.         .rare = RareAbundanceRarecurve, 
    79.         .alpha = Observe,
    80.       )
    82. p2 <- mpse_abund %>% 
    83.       mp_plot_rarecurve(
    84.         .rare = RareAbundanceRarecurve, 
    85.         .alpha = Observe, 
    86.         .group = pre_post_stroke
    87.       ) +
    88.       scale_color_manual(values=c("#00A087FF", "#3C5488FF")) +
    89.       scale_fill_manual(values=c("#00A087FF", "#3C5488FF"), guide="none")
    91. # combine the samples belong to the same groups if plot.group=TRUE
    92. p3 <- mpse_abund %>% 
    93.       mp_plot_rarecurve(
    94.         .rare = RareAbundanceRarecurve, 
    95.         .alpha = "Observe", 
    96.         .group = pre_post_stroke, 
    97.         plot.group = TRUE
    98.       ) +
    99.       scale_color_manual(values=c("#00A087FF", "#3C5488FF")) +
    100.       scale_fill_manual(values=c("#00A087FF", "#3C5488FF"),guide="none")
    102. png("rarefaction_of_samples_or_groups.png", width=1080, height=600)
    103. p1 + p2 + p3
    104. dev.off()
    106. # ------------------------------------------
    107. # 3.3. calculate alpha index and visualization
    108. library(ggplot2)
    109. library(MicrobiotaProcess)
    110. mpse_abund %<>% 
    111.     mp_cal_alpha(.abundance=RareAbundance)
    112. mpse_abund
    114. f1 <- mpse_abund %>% 
    115.       mp_plot_alpha(
    116.         .group=pre_post_stroke, 
    117.         .alpha=c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)
    118.       ) +
    119.       scale_fill_manual(values=c("#00A087FF", "#3C5488FF"), guide="none") +
    120.       scale_color_manual(values=c("#00A087FF", "#3C5488FF"), guide="none")
    122. f2 <- mpse_abund %>%
    123.       mp_plot_alpha(
    124.         .alpha=c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)
    125.       )
    126. png("alpha_diversity_comparison.png", width=1000, height=1000)
    127. f1 / f2
    128. dev.off()
    130. # -------------------------------------------
    131. # 3.4. The visualization of taxonomy abundance (Class)
    132. mpse_abund %<>%
    133.     mp_cal_abundance( # for each samples
    134.       .abundance = RareAbundance
    135.     ) %>%
    136.     mp_cal_abundance( # for each groups 
    137.       .abundance=RareAbundance,
    138.       .group=pre_post_stroke
    139.     )
    140. mpse_abund
    142. # visualize the relative abundance of top 20 phyla for each sample.
    143. p1 <- mpse_abund %>%
    144.         mp_plot_abundance(
    145.           .abundance=RareAbundance,
    146.           .group=time, 
    147.           taxa.class = Class, 
    148.           topn = 20,
    149.           relative = TRUE
    150.         )
    151. # visualize the abundance (rarefied) of top 20 phyla for each sample.
    152. p2 <- mpse_abund %>%
    153.           mp_plot_abundance(
    154.             .abundance=RareAbundance,
    155.             .group=time,
    156.             taxa.class = Class,
    157.             topn = 20,
    158.             relative = FALSE
    159.           )
    160. png("relative_abundance_and_abundance.png", width= 1200, height=600) #NOT PRODUCED!
    161. p1 / p2
    162. dev.off()
    164. h1 <- mpse_abund %>%
    165.         mp_plot_abundance(
    166.           .abundance = RareAbundance,
    167.           .group = pre_post_stroke,
    168.           taxa.class = Class,
    169.           relative = TRUE,
    170.           topn = 20,
    171.           geom = 'heatmap',
    172.           features.dist = 'euclidean',
    173.           features.hclust = 'average',
    174.           sample.dist = 'bray',
    175.           sample.hclust = 'average'
    176.         )
    178. h2 <- mpse_abund %>%
    179.           mp_plot_abundance(
    180.             .abundance = RareAbundance,
    181.             .group = pre_post_stroke,
    182.             taxa.class = Class,
    183.             relative = FALSE,
    184.             topn = 20,
    185.             geom = 'heatmap',
    186.             features.dist = 'euclidean',
    187.             features.hclust = 'average',
    188.             sample.dist = 'bray',
    189.             sample.hclust = 'average'
    190.           )
    191. # the character (scale or theme) of figure can be adjusted by set_scale_theme
    192. # refer to the mp_plot_dist
    193. png("relative_abundance_and_abundance_heatmap.png", width= 1200, height=600)
    194. aplot::plot_list(gglist=list(h1, h2), tag_levels="A")
    195. dev.off()
    197. # visualize the relative abundance of top 20 class for each .group (pre_post_stroke)
    198. p3 <- mpse_abund %>%
    199.         mp_plot_abundance(
    200.             .abundance=RareAbundance, 
    201.             .group=pre_post_stroke,
    202.             taxa.class = Class,
    203.             topn = 20,
    204.             plot.group = TRUE
    205.           )
    207. # visualize the abundance of top 20 phyla for each .group (time)
    208. p4 <- mpse_abund %>%
    209.           mp_plot_abundance(
    210.             .abundance=RareAbundance,
    211.             .group= pre_post_stroke,
    212.             taxa.class = Class,
    213.             topn = 20,
    214.             relative = FALSE,
    215.             plot.group = TRUE
    216.           )
    217. png("relative_abundance_and_abundance_groups.png", width= 1000, height=1000)
    218. p3 / p4
    219. dev.off()
    221. # ---------------------------
    222. # 3.5. Beta diversity analysis
    224. # ---------------------------------------------
    225. # 3.5.1 The distance between samples or groups
    226. # standardization
    227. # mp_decostand wraps the decostand of vegan, which provides
    228. # many standardization methods for community ecology.
    229. # default is hellinger, then the abundance processed will
    230. # be stored to the assays slot. 
    231. mpse_abund %<>% 
    232.     mp_decostand(.abundance=Abundance)
    233. mpse_abund
    235. # calculate the distance between the samples.
    236. # the distance will be generated a nested tibble and added to the
    237. # colData slot.
    238. mpse_abund %<>% mp_cal_dist(.abundance=hellinger, distmethod="bray")
    239. mpse_abund
    241. # mp_plot_dist provides there methods to visualize the distance between the samples or groups
    242. # when .group is not provided, the dot heatmap plot will be return
    243. p1 <- mpse_abund %>% mp_plot_dist(.distmethod = bray)
    244. png("distance_between_samples.png", width= 1000, height=1000)
    245. p1
    246. dev.off()
    248. # when .group is provided, the dot heatmap plot with group information will be return.
    249. p2 <- mpse_abund %>% mp_plot_dist(.distmethod = bray, .group = pre_post_stroke)
    250. # The scale or theme of dot heatmap plot can be adjusted using set_scale_theme function.
    251. p2 %>% set_scale_theme(
    252.           x = scale_fill_manual(
    253.                 values=c("orange", "deepskyblue"), 
    254.                 guide = guide_legend(
    255.                             keywidth = 1, 
    256.                             keyheight = 0.5, 
    257.                             title.theme = element_text(size=8),
    258.                             label.theme = element_text(size=6)
    259.                 )
    260.               ), 
    261.           aes_var = pre_post_stroke # specific the name of variable 
    262.       ) %>%
    263.       set_scale_theme(
    264.           x = scale_color_gradient(
    265.                 guide = guide_legend(keywidth = 0.5, keyheight = 0.5)
    266.               ),
    267.           aes_var = bray
    268.       ) %>%
    269.       set_scale_theme(
    270.           x = scale_size_continuous(
    271.                 range = c(0.1, 3),
    272.                 guide = guide_legend(keywidth = 0.5, keyheight = 0.5)
    273.               ),
    274.           aes_var = bray
    275.       )
    276. png("distance_between_samples_with_group_info.png", width= 1000, height=1000)
    277. p2
    278. dev.off()
    280. # when .group is provided and group.test is TRUE, the comparison of different groups will be returned
    281. p3 <- mpse_abund %>% mp_plot_dist(.distmethod = bray, .group = pre_post_stroke, group.test=TRUE, textsize=2)
    282. png("comparison_of_distance.png", width= 1000, height=1000)
    283. p3
    284. dev.off()
    286. # -----------------------
    287. # 3.5.2 The PCoA analysis
    289. #install.packages("corrr")
    290. library(corrr)
    291. #install.packages("ggside")
    292. library(ggside)
    293. mpse_abund %<>% 
    294.     mp_cal_pcoa(.abundance=hellinger, distmethod="bray")
    295. # The dimensions of ordination analysis will be added the colData slot (default).
    296. mpse_abund
    297. #> methods(class=class(mpse_abund))
    298. # [1] [                        [[<-                     [<-                     
    299. # [4] $                        $<-                      arrange                 
    300. # [7] as_tibble                as.data.frame            as.phyloseq             
    301. #[10] coerce                   coerce<-                 colData<-               
    302. #[13] distinct                 filter                   group_by                
    303. #[16] left_join                mp_adonis                mp_aggregate_clade      
    304. #[19] mp_aggregate             mp_anosim                mp_balance_clade        
    305. #[22] mp_cal_abundance         mp_cal_alpha             mp_cal_cca              
    306. #[25] mp_cal_clust             mp_cal_dca               mp_cal_dist             
    307. #[28] mp_cal_nmds              mp_cal_pca               mp_cal_pcoa             
    308. #[31] mp_cal_pd_metric         mp_cal_rarecurve         mp_cal_rda              
    309. #[34] mp_cal_upset             mp_cal_venn              mp_decostand            
    310. #[37] mp_diff_analysis         mp_diff_clade            mp_envfit               
    311. #[40] mp_extract_abundance     mp_extract_assays        mp_extract_dist         
    312. #[43] mp_extract_feature       mp_extract_internal_attr mp_extract_rarecurve    
    313. #[46] mp_extract_refseq        mp_extract_sample        mp_extract_taxonomy     
    314. #[49] mp_extract_tree          mp_filter_taxa           mp_mantel               
    315. #[52] mp_mrpp                  mp_plot_abundance        mp_plot_alpha           
    316. #[55] mp_plot_diff_boxplot     mp_plot_diff_res         mp_plot_dist            
    317. #[58] mp_plot_ord              mp_plot_rarecurve        mp_plot_upset           
    318. #[61] mp_plot_venn             mp_rrarefy               mp_select_as_tip        
    319. #[64] mp_stat_taxa             mutate                   otutree                 
    320. #[67] otutree<-                print                    pull                    
    321. #[70] refsequence              refsequence<-            rename                  
    322. #[73] rownames<-               select                   show                    
    323. # [ reached getOption("max.print") -- omitted 6 entries ]
    324. #see '?methods' for accessing help and source code
    326. # We also can perform adonis or anosim to check whether it is significant to the dissimilarities of groups.
    327. mpse_abund %<>%
    328.     mp_adonis(.abundance=hellinger, .formula=~Group, distmethod="bray", permutations=9999, action="add")
    329. mpse_abund %>% mp_extract_internal_attr(name=adonis)
    331. # ("1","2","5","6","7",  "15","16","17","18","19","20",  "29","30","31","32",  "40","41","42","43","44","46")
    332. #div.df2[div.df2 == "Group1"] <- "aged.post"
    333. #div.df2[div.df2 == "Group3"] <- "young.post"
    334. #div.df2[div.df2 == "Group5"] <- "aged.post"
    335. #div.df2[div.df2 == "Group7"] <- "young.post"
    337. # ("8","9","10","12","13","14",  "21","22","23","24","25","26","27","28",  "33","34","35","36","37","38","39","51",  "47","48","49","50","52","53","55")
    338. #div.df2[div.df2 == "Group2"] <- "aged.pre"
    339. #div.df2[div.df2 == "Group4"] <- "young.pre"
    340. #div.df2[div.df2 == "Group6"] <- "aged.pre"
    341. #div.df2[div.df2 == "Group8"] <- "young.pre"
    343. #Group1: f.aged and post
    344. #Group2: f.aged and pre
    345. #Group3: f.young and post
    346. #Group4: f.young and pre
    347. #Group5: m.aged and post
    348. #Group6: m.aged and pre
    349. #Group7: m.young and post
    350. #Group8: m.young and pre
    352. #[,c("1","2","5","6","7",                "8","9","10","12","13","14")]
    353. #[,c("15","16","17","18","19","20",      "21","22","23","24","25","26","27","28")]
    354. #[,c("29","30","31","32",                "33","34","35","36","37","38","39","51")]
    355. #[,c("40","41","42","43","44","46",      "47","48","49","50","52","53","55")]
    357. p1 <- mpse_abund %>%
    358.         mp_plot_ord(
    359.           .ord = pcoa, 
    360.           .group = Group, 
    361.           .color = Group, 
    362.           .size = 2.4,
    363.           .alpha = 1,
    364.           ellipse = TRUE,
    365.           show.legend = FALSE # don't display the legend of stat_ellipse
    366.         ) +
    367.         scale_fill_manual(
    368.           #values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#cab2d6", "#6a3d9a"), 
    369.           #values = c("#a6cee3", "#b2df8a", "#fb9a99", "#cab2d6"), 
    370.           values = c("#1f78b4", "#33a02c", "#e31a1c", "#6a3d9a"), 
    371.           guide = guide_legend(keywidth=1.6, keyheight=1.6, label.theme=element_text(size=12))
    372.         ) +
    373.         scale_color_manual(
    374.           #values=c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#cab2d6", "#6a3d9a"),
    375.           #values = c("#a6cee3", "#b2df8a", "#fb9a99", "#cab2d6"),
    376.           values = c("#1f78b4", "#33a02c", "#e31a1c", "#6a3d9a"), 
    377.           guide = guide_legend(keywidth=1.6, keyheight=1.6, label.theme=element_text(size=12))
    378.         )
    379.         #scale_fill_manual(values=c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#00FF00", "#FFFF00", "#00FFFF", "#FFA500")) +
    380.         #scale_color_manual(values=c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#00FF00", "#FFFF00", "#00FFFF", "#FFA500"))
    381.         #scale_fill_manual(values=c("#00A087FF", "#3C5488FF")) +
    382.         #scale_color_manual(values=c("#00A087FF", "#3C5488FF")) 
    383. #png("PCoA.png", width= 1000, height=1000)
    384. svg("PCoA.svg", width= 11, height=10)
    385. #svg("PCoA_.svg", width=10, height=10)
    386. #svg("PCoA.svg")
    387. pdf("PCoA.pdf")
    388. p1
    389. dev.off()
    391. # The size of point also can be mapped to other variables such as Observe, or Shannon 
    392. # Then the alpha diversity and beta diversity will be displayed simultaneously.
    393. p2 <- mpse_abund %>% 
    394.         mp_plot_ord(
    395.           .ord = pcoa, 
    396.           .group = Group, 
    397.           .color = Group, 
    398.           .size = Shannon,
    399.           .alpha = Observe,
    400.           ellipse = TRUE,
    401.           show.legend = FALSE # don't display the legend of stat_ellipse 
    402.         ) +
    403.         scale_fill_manual(
    404.           values = c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#00FF00", "#FFFF00", "#00FFFF", "#FFA500"), 
    405.           guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=8))
    406.         ) +
    407.         scale_color_manual(
    408.           values=c("#FF0000", "#000000", "#0000FF", "#C0C0C0", "#00FF00", "#FFFF00", "#00FFFF", "#FFA500"),
    409.           guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=8))
    410.         ) +
    411.         scale_size_continuous(
    412.           range=c(0.5, 3),
    413.           guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=8))
    414.         )
    416. # ------------------------------------------
    417. # 3.5.3 Hierarchical cluster (tree) analysis
    418. #input should contain hellinger!
    419. mpse_abund %<>%
    420.       mp_cal_clust(
    421.         .abundance = hellinger, 
    422.         distmethod = "bray",
    423.         hclustmethod = "average", # (UPGAE)
    424.         action = "add" # action is used to control which result will be returned
    425.       )
    426. mpse_abund
    428. # if action = 'add', the result of hierarchical cluster will be added to the MPSE object
    429. # mp_extract_internal_attr can extract it. It is a treedata object, so it can be visualized
    430. # by ggtree.
    431. sample.clust <- mpse_abund %>% mp_extract_internal_attr(name='SampleClust')
    432. sample.clust
    434. library(ggtree)
    435. p_cluster <- ggtree(sample.clust) + 
    436.       geom_tippoint(aes(color=pre_post_stroke)) +
    437.       geom_tiplab(as_ylab = TRUE) +
    438.       ggplot2::scale_x_continuous(expand=c(0, 0.01))
    439. png("hierarchical_cluster1.png", width= 1000, height=800)
    440. p_cluster
    441. dev.off()
    443. library(ggtreeExtra)
    444. library(ggplot2)
    445. # # Extract relative abundance of phyla
    446. # phyla.tb <- mouse.time.mpse %>% 
    447. #             mp_extract_abundance(taxa.class=Phylum, topn=30)
    448. # # The abundance of each samples is nested, it can be flatted using the unnest of tidyr.
    449. # phyla.tb %<>% tidyr::unnest(cols=RareAbundanceBySample) %>% dplyr::rename(Phyla="label")
    450. # phyla.tb
    451. #
    452. # p1 <- p + 
    453. #       geom_fruit(
    454. #          data=phyla.tb,
    455. #          geom=geom_col,
    456. #          mapping = aes(x = RelRareAbundanceBySample, 
    457. #                        y = Sample, 
    458. #                        fill = Phyla
    459. #                  ),
    460. #          orientation = "y",
    461. #          #offset = 0.4,
    462. #          pwidth = 3, 
    463. #          axis.params = list(axis = "x", 
    464. #                             title = "The relative abundance of phyla (%)",
    465. #                             title.size = 4,
    466. #                             text.size = 2, 
    467. #                             vjust = 1),
    468. #          grid.params = list()
    469. #       )
    470. # png("hierarchical_cluster2.png", width = 1000, height = 800)
    471. # p1
    472. # dev.off()
    474. # Extract relative abundance of classes
    475. mpse_abund %>% print(width=150)
    476. class.tb <- mpse_abund %>% 
    477.             mp_extract_abundance(taxa.class = Class, topn = 30)
    478. # Flatten and rename the columns
    479. class.tb %<>% tidyr::unnest(cols = RareAbundanceBySample) %>% dplyr::rename(Class = "label")
    480. # View the data frame
    481. class.tb
    482. # Create the plot
    483. p1 <- p + 
    484.       geom_fruit(
    485.         data = class.tb,
    486.         geom = geom_col,
    487.         mapping = aes(x = RelRareAbundanceBySample, 
    488.                       y = Sample, 
    489.                       fill = Class
    490.                 ),
    491.         orientation = "y",
    492.         pwidth = 3, 
    493.         axis.params = list(axis = "x", 
    494.                             title = "The relative abundance of classes (%)",
    495.                             title.size = 4,
    496.                             text.size = 2, 
    497.                             vjust = 1),
    498.         grid.params = list()
    499.       )
    501. # Save the plot to a file
    502. png("hierarchical_cluster2.png", width = 1000, height = 800)
    503. print(p1)
    504. dev.off()
    506. # -----------------------
    507. # 3.6 Biomarker discovery
    508. library(ggtree)
    509. library(ggtreeExtra)
    510. library(ggplot2)
    511. library(MicrobiotaProcess)
    512. library(tidytree)
    513. library(ggstar)
    514. library(forcats)
    515. mpse_abund %>% print(width=150)
    517. mpse_abund %<>%
    518.     mp_cal_abundance( # for each samples
    519.       .abundance = RareAbundance
    520.     ) %>%
    521.     mp_cal_abundance( # for each groups 
    522.       .abundance=RareAbundance,
    523.       .group=pre_post_stroke
    524.     )
    525. mpse_abund
    527. mpse_abund %<>%
    528.     mp_diff_analysis(
    529.       .abundance = RelRareAbundanceBySample,
    530.       .group = pre_post_stroke,
    531.       first.test.alpha = 0.01
    532.     )
    534. # The result is stored to the taxatree or otutree slot, you can use mp_extract_tree to extract the specific slot.
    535. taxa.tree <- mpse_abund %>% 
    536.               mp_extract_tree(type="taxatree")
    537. taxa.tree
    539. ## And the result tibble of different analysis can also be extracted 
    540. ## with tidytree (>=0.3.5)
    541. taxa.tree %>% select(label, nodeClass, LDAupper, LDAmean, LDAlower, Sign_pre_post_stroke, pvalue, fdr) %>% dplyr::filter(!is.na(fdr))
    542. taxa.tree %>% print(width=150, n=100)
    544. #.data, layout, tree.type, .taxa.class, tiplab.size, offset.abun, pwidth.abun, offset.effsize, pwidth.effsize, group.abun, tiplab.linetype
    545. p <- mpse_abund %>%
    546.       mp_plot_diff_res(
    547.         group.abun = TRUE,
    548.         pwidth.abun=0.05,
    549.         offset.abun=0.02,
    550.         pwidth.effsize=0.3,
    551.         offset.effsize=0.46,
    552.         tiplab.size = 4.9
    553.       ) +
    554.       scale_fill_manual(values=c("deepskyblue", "orange")) +
    555.       scale_fill_manual(
    556.         aesthetics = "fill_new", # The fill aes was renamed to "fill_new" for the abundance dotplot layer
    557.         values = c("deepskyblue", "orange")
    558.       ) +
    559.       scale_fill_manual(
    560.         aesthetics = "fill_new_new", # The fill aes for hight light layer of tree was renamed to 'fill_new_new'
    561.         values = c("#E41A1C", "#377EB8", "#4DAF4A",
    562.                     "#984EA3", "#FF7F00", "#FFFF33",
    563.                     "#A65628", "#F781BF", "#00FFFF", "#999999"
    564.                   )
    565.       ) + 
    566.       theme(
    567.     axis.title = element_text(size = 28),        # Font size for axis titles
    568.     axis.text = element_text(size = 28),          # Font size for axis text
    569.     plot.title = element_text(size = 28),         # Font size for plot title
    570.     legend.title = element_text(size = 16),       # Font size for legend title
    571.     legend.text = element_text(size = 14)         # Font size for legend text
    572.   )
    573. #p$layers[[2]]$geom <- geom_tiplab(fontsize = 22)  # Change 12 to the desired font size
    575. png("differently_expressed_otu.png", width=2000, height=2000)
    576. #svg("p7.svg",width=8, height=8)
    577. p
    578. dev.off()
    580. f <- mpse_abund %>%
    581.     mp_plot_diff_cladogram(
    582.       label.size = 2.5,
    583.       hilight.alpha = .3,
    584.       bg.tree.size = .5,
    585.       bg.point.size = 2,
    586.       bg.point.stroke = .25
    587.     ) +
    588.     scale_fill_diff_cladogram( # set the color of different group.
    589.       values = c('deepskyblue', 'orange')
    590.     ) +
    591.     scale_size_continuous(range = c(1, 4))
    592. #png("cladogram.png", width=1000, height=1000)
    593. svg("cladogram.svg", width=10, height=10)
    594. f
    595. dev.off()
    597. ## Extract the OTU table and taxonomy table from the phyloseq object
    598. #otu_table <- phyloseq::otu_table(ps.ng.tax_abund) %>% as.data.frame() %>% as.matrix()
    599. #tax_table <- phyloseq::tax_table(ps.ng.tax_abund) %>% as.data.frame() %>% as.matrix()
    600. #write.csv(otu_table, file="otu_table.csv")
    601. #write.csv(tax_table, file="tax_table.csv")
    602. #~/Tools/csv2xls-0.4/csv_to_xls.py otu_table.csv tax_table.csv -d',' -o otu_tax.xls

3.1. Environment Setup: It sets up a Conda environment named picrust2, using the conda create command and then activates this environment using conda activate picrust2.

  1.     #https://github.com/picrust/picrust2/wiki/PICRUSt2-Tutorial-(v2.2.0-beta)#minimum-requirements-to-run-full-tutorial
  2.     conda create -n picrust2 -c bioconda -c conda-forge picrust2=2.2.0_b
  3.     conda activate picrust2

3.2. Data Preparation: The script creates a new directory called picrust2_out, then enters it using mkdir and cd commands. It then identifies input files that are needed for the analysis: metadata.tsv, seqs.fna, table.biom. The biom commands are used to inspect and convert the BIOM format files.

  1.     mkdir picrust2_out
  2.     cd picrust2_out
  4.     # Identifying input data
  5.     # Note: Replace the paths and filenames with your actual data if different
  6.     # metadata.tsv == ../map_corrected.txt
  7.     # seqs.fna     == ../clustering/seqs.fna
  8.     # table.biom   == ../core_diversity_e42369/table_even42369.biom
  10.     # Inspect and convert the BIOM format files
  11.     biom head -i ../core_diversity_e42369/table_even42369.biom
  12.     biom summarize-table -i ../core_diversity_e42369/table_even42369.biom
  13.     biom convert -i ../core_diversity_e42369/table_even42369.biom -o table_even42369.tsv --to-tsv

3.3. Running PiCRUST2: The place_seqs.py command aligns the input sequences to a reference tree. The hsp.py commands generate hidden state prediction for multiple functional categories.

  1.     #insert reads into reference tree using EPA-NG
  2.     cp ../clustering/rep_set.fna ./
  3.     grep ">" rep_set.fna | wc -l  #44238
  4.     vim table_even42369.tsv       #40596-2
  6.     samtools faidx rep_set.fna
  7.     cut -f1-1 table_even42369.tsv > table_even42369.id
  8.     #manually modify table_even42369.id by replacing "\n" with " >> seqs.fna\nsamtools faidx rep_set.fna "
  9.     run table_even42369.id
  11.     rm -rf intermediate/
  12.     place_seqs.py -s seqs.fna -o out.tre -p 64 --intermediate intermediate/place_seqs
  14.     #castor: Efficient Phylogenetics on Large Trees
  15.     #https://github.com/picrust/picrust2/wiki/Hidden-state-prediction
  17.     hsp.py -i 16S -t out.tre -o 16S_predicted_and_nsti.tsv.gz -p 100 -n
  18.     hsp.py -i COG -t out.tre -o COG_predicted.tsv.gz -p 100
  19.     hsp.py -i PFAM -t out.tre -o PFAM_predicted.tsv.gz -p 100
  20.     hsp.py -i KO -t out.tre -o KO_predicted.tsv.gz -p 100
  21.     hsp.py -i EC -t out.tre -o EC_predicted.tsv.gz -p 100
  22.     hsp.py -i TIGRFAM -t out.tre -o TIGRFAM_predicted.tsv.gz -p 100
  23.     hsp.py -i PHENO -t out.tre -o PHENO_predicted.tsv.gz -p 100
  25. >In this table the predicted copy number of all Enzyme Classification (EC) numbers is shown for each ASV. The NSTI values per ASV are not in this table since we did not specify the -n option. EC numbers are a type of gene family defined based on the chemical reactions they catalyze. For instance, EC:1.1.1.1 corresponds to alcohol dehydrogenase. In this tutorial we are focusing on EC numbers since they can be used to infer MetaCyc pathway levels (see below).
  27.     zless -S EC_predicted.tsv.gz
  28.     sequence        EC:1.1.1.1      EC:1.1.1.10     EC:1.1.1.100    ...
  29.     20e568023c10eaac834f1c110aacea18        2       0       3    ...
  30.     23fe12a325dfefcdb23447f43b6b896e        0       0       1    ...
  31.     288c8176059111c4c7fdfb0cd5afce64        1       0       1    ...
  32.     ...
  34.     ##Why import the tsv file to MyData?
  35.     #MyData <- read.csv(file="./COG_predicted.tsv", header=TRUE, sep="\t", row.names=1)   #6806 4598  e.g. COG5665
  36.     #MyData <- read.csv(file="./PFAM_predicted.tsv", header=TRUE, sep="\t", row.names=1)  #6806 11089 e.g. PF17225
  37.     #MyData <- read.csv(file="./KO_predicted.tsv", header=TRUE, sep="\t", row.names=1)    #6806 10543 e.g. K19791
  38.     #MyData <- read.csv(file="./EC_predicted.tsv", header=TRUE, sep="\t", row.names=1)    #6806 2913  e.g. EC.6.6.1.2
  39.     #MyData <- read.csv(file="./16S_predicted.tsv", header=TRUE, sep="\t", row.names=1)   #6806    1     e.g. X16S_rRNA_Count
  40.     #MyData <- read.csv(file="./TIGRFAM_predicted.tsv", header=TRUE, sep="\t", row.names=1)  #6806 4287  e.g. TIGR04571
  41.     #MyData <- read.csv(file="./PHENO_predicted.tsv", header=TRUE, sep="\t", row.names=1)    #6806   41  e.g. Use_of_nitrate_as_electron_acceptor, Xylose_utilizing

3.4. The metagenome_pipeline.py commands perform metagenomic prediction for several functional categories. Predicted gene families weighted by the relative abundance of ASVs in their community. In other words, we are interested in inferring the metagenomes of the communities.

  1.     #Generate metagenome predictions using EC numbers https://en.wikipedia.org/wiki/List_of_enzymes#Category:EC_1.1_(act_on_the_CH-OH_group_of_donors)
  2.     metagenome_pipeline.py -i ../core_diversity_e42369/table_even42369.biom -m 16S_predicted_and_nsti.tsv.gz -f COG_predicted.tsv.gz -o COG_metagenome_out --strat_out
  3.     metagenome_pipeline.py -i ../core_diversity_e42369/table_even42369.biom -m 16S_predicted_and_nsti.tsv.gz -f EC_predicted.tsv.gz -o EC_metagenome_out --strat_out
  4.     metagenome_pipeline.py -i ../core_diversity_e42369/table_even42369.biom -m 16S_predicted_and_nsti.tsv.gz -f KO_predicted.tsv.gz -o KO_metagenome_out --strat_out
  5.     metagenome_pipeline.py -i ../core_diversity_e42369/table_even42369.biom -m 16S_predicted_and_nsti.tsv.gz -f PFAM_predicted.tsv.gz -o PFAM_metagenome_out --strat_out
  6.     metagenome_pipeline.py -i ../core_diversity_e42369/table_even42369.biom -m 16S_predicted_and_nsti.tsv.gz -f TIGRFAM_predicted.tsv.gz -o TIGRFAM_metagenome_out --strat_out

3.5. Pathway-level inference: By default this script infers MetaCyc pathway abundances based on EC number abundances, although different gene families and pathways can also be optionally specified. This script performs a number of steps by default, which are based on the approach implemented in HUMAnN2:

  • Regroups EC numbers to MetaCyc reactions.
  • Infers which MetaCyc pathways are present based on these reactions with MinPath.
  • Calculates and returns the abundance of pathways identified as present.
    1. pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_contrib.tsv.gz -o pathways_out -p 15
    3. #Note that the path of map files is under /home/jhuang/anaconda3/envs/picrust2/lib/python3.6/site-packages/picrust2/default_files/pathway_mapfiles
    4. pathway_pipeline.py -i COG_metagenome_out/pred_metagenome_contrib.tsv.gz -o KEGG_pathways_out -p 15 --no_regroup --map /home/jhuang/anaconda3/envs/picrust2/lib/python3.6/site-packages/picrust2/default_files/pathway_mapfiles/KEGG_pathways_to_KO.tsv
    6. #Mapping predicted KO abundances to legacy KEGG pathways (with stratified output that represents contributions to community-wide abundances):
    7. pathway_pipeline.py -i KO_metagenome_out/pred_metagenome_strat.tsv.gz -o KEGG_pathways_out --no_regroup --map /home/jhuang/anaconda3/envs/picrust2/lib/python3.6/site-packages/picrust2/default_files/pathway_mapfiles/KEGG_pathways_to_KO.tsv
    9. #Map EC numbers to MetaCyc pathways and get stratified output corresponding to contribution of predicted gene family abundances within each predicted genome:
    10. pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -o pathways_out_per_seq --per_sequence_contrib --per_sequence_abun EC_metagenome_out/seqtab_norm.tsv.gz --per_sequence_function EC_predicted.tsv.gz

3.6. Add functional descriptions: Finally, it can be useful to have a description of each functional id in the output abundance tables. The below commands will add these descriptions as new column in gene family and pathway abundance tables

  1.     #--6.1. Add descriptions in gene family tables
  2.     add_descriptions.py -i COG_metagenome_out/pred_metagenome_unstrat.tsv.gz -m COG -o COG_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
  3.     add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC -o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
  4.     add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO -o KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz   # EC and METACYC is a pair, EC for gene_annotation and METACYC for pathway_annotation
  5.     add_descriptions.py -i PFAM_metagenome_out/pred_metagenome_unstrat.tsv.gz -m PFAM -o PFAM_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
  6.     add_descriptions.py -i TIGRFAM_metagenome_out/pred_metagenome_unstrat.tsv.gz -m TIGRFAM -o TIGRFAM_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
  8.     #--6.2. Add descriptions in pathway abundance tables
  9.     add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC -o pathways_out/path_abun_unstrat_descrip.tsv.gz
  10.     gunzip path_abun_unstrat_descrip.tsv.gz
  12.     #Error - no rows remain after regrouping input table. The default pathway and regroup mapfiles are meant for EC numbers. Note that KEGG pathways are not supported since KEGG is a closed-source database, but you can input custom pathway mapfiles if you have access. If you are using a custom function database did you mean to set the --no-regroup flag and/or change the default pathways mapfile used?
  13.     #If ERROR --> USE the METACYC for downstream analyses!!!
  15.     add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -o KEGG_pathways_out/path_abun_unstrat_descrip.tsv.gz --custom_map_table /home/jhuang/anaconda3/envs/picrust2/lib/python3.6/site-packages/picrust2/default_files/description_mapfiles/KEGG_pathways_info.tsv.gz

3.7. Visualization

  1.     #7.1 STAMP
  2.     #https://github.com/picrust/picrust2/wiki/STAMP-example
  3.     conda activate stamp; #import path_abun_unstrat_descrip.tsv to STAMP
  4.     conda deactivate
  5.     conda install -c bioconda stamp
  7.     #conda install -c bioconda stamp
  8.     #sudo pip install pyqi
  9.     #sudo apt-get install libblas-dev liblapack-dev gfortran
  10.     #sudo apt-get install freetype* python-pip python-dev python-numpy python-scipy python-matplotlib
  11.     #sudo pip install STAMP
  12.     #conda install -c bioconda stamp
  14.     conda create -n stamp -c bioconda/label/cf201901 stamp
  15.     brew install pyqt
  17.     #DEBUG the environment
  18.     conda install pyqt=4
  19.     #conda install icu=56
  21.     e.g. path_abun_unstrat_descrip.tsv.gz and metadata.tsv from the tutorial)
  22.     cut -d$'\t' -f1 map_corrected.txt > 1
  23.     cut -d$'\t' -f5 map_corrected.txt > 5
  24.     cut -d$'\t' -f6 map_corrected.txt > 6
  25.     paste -d$'\t' 1 5 > 1_5
  26.     paste -d$'\t' 1_5 6 > metadata.tsv
  27.     #SampleID --> SampleID
  28.     SampleID    Facility    Genotype
  29.     100CHE6KO   PaloAlto    KO
  30.     101CHE6WT   PaloAlto    WT
  32.     #7.2. ALDEx2
  33.     https://bioconductor.org/packages/release/bioc/html/ALDEx2.html
  35.     #7.3. Convert png to svg and pdf
  36.     inkscape error_bar.png --export-plain-svg=error_bar.svg (embbed)
  38.     sudo apt update
  39.     sudo apt install autotrace
  41.     sudo apt-get install -y libpng-dev libtiff-dev imagemagick
  42.     git clone https://github.com/autotrace/autotrace.git
  44.     cd autotrace
  45.     #sudo apt install intltool
  46.     #sudo apt install gettext libglib2.0-dev
  47.     #sudo apt install libtool libtool-bin
  48.     #sudo apt install automake
  49.     sudo apt-get install libxml-parser-perl
  50.     ./autogen.sh
  51.     ./configure
  52.     make
  54.     autotrace -output-format svg -output-file error_bar.svg error_bar.png

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. pheatmap vs heatmap.2
  6. Setup conda environments
  7. Should the inputs for GSVA be normalized or raw?
  8. File format for single channel analysis of Agilent microarray data with Limma?
  9. PiCRUST2 Pipeline for Functional Prediction and Pathway Analysis in Metagenomics
  10. Kraken2 Installation and Usage Guide

最新文章

最多评论文章

推荐相似文章

QIIME + Phyloseq + MicrobiotaProcess (v1)

Plotting Alpha Diversities from 16S rRNA Sequencing Data

Processing for Data_Marius_16S_2025

MicrobiotaProcess Group2 vs Group6 (v1)

© 2023 XGenes.com Impressum