16S-rDNA-V3-V4 repository: https://github.com/ycl6/16S-rDNA-V3-V4

export2graphlan: GitHub, Documentation

GraPhlAn: GitHub, Documentation, Paper

Demo Dataset: PRJEB27564 from Gut Microbiota in Parkinson’s Disease: Temporal Stability and Relations to Disease Progression. EBioMedicine. 2019;44:691-707

License: GPL-3.0

Workflow (Continues from Part 2)

Note: export2graphlan and GraPhlAn can only run in Python 2.7.

Start `R`

cd /ngs/16S-Demo/PRJEB27564

R

Load package and set path

library("data.table")
library("phyloseq")
library("ggplot2")
library("grid")
library("tidyverse")
Sys.setlocale("LC_COLLATE", "C")

## [1] "C"

Set the full-path to additional scripts

source("/path-to-script/lefse.R", local = TRUE)

Load saved workspace from Part 2

load("2_phyloseq_tutorial.RData")

Create a folder lefse

if(!dir.exists("lefse")) { dir.create("lefse") }

Prepare input data

Subset samples

# Patient vs. control at baseline
ps1a = prune_samples(sample_names(ps1)[grep("B$", sample_names(ps1))], ps1)
ps1a = prune_taxa(taxa_sums(ps1a) > 0, ps1a)

# Patient followup vs. patient baseline
ps1b = prune_samples(sample_names(ps1)[grep("^P", sample_names(ps1))], ps1)
ps1b = prune_taxa(taxa_sums(ps1b) > 0, ps1b)

Prepare LEfSe input

We are now using new version of LEfSe. Since we are interested in testing the biomarker using the main class labels (control vs. patient), hence the subject line is removed when creating the LEfSe input table.

# Add 'subject = TRUE' to include subject in the data.frame
tax1 = lefse_1name_obj(ps1a, sample_data(ps1a)$Group)
lefse1 = lefse_obj(ps1a)
lefse1 = rbind(tax1, lefse1)

tax2 = lefse_1name_obj(ps1b, sample_data(ps1b)$Time)
lefse2 = lefse_obj(ps1b)
lefse2 = rbind(tax2, lefse2)

# Replace unsupported chars with underscore
lefse1$name = gsub(" ","_",lefse1$name)
lefse1$name = gsub("-","_",lefse1$name)
lefse1$name = gsub("/","_",lefse1$name)

lefse2$name = gsub(" ","_",lefse2$name)
lefse2$name = gsub("-","_",lefse2$name)
lefse2$name = gsub("/","_",lefse2$name)

For Silva v132 users, apply below to fix identifical phylum & class names (Actinobacteria and Deferribacteres)

lefse1 = fix_taxa_silva132(lefse1)
lefse2 = fix_taxa_silva132(lefse2)

Output to file

expr1: Patient vs. control at baseline
expr2: Patient at baseline vs. patient at followup

write.table(lefse1, file = "lefse/expr1.lefse_table.txt", sep = "\t", quote = F, 
        row.names = F, col.names = F)
write.table(lefse2, file = "lefse/expr2.lefse_table.txt", sep = "\t", quote = F, 
        row.names = F, col.names = F)

Perform LEfSe analysis

We use linear discriminant analysis (LDA) effect size (LEfSe) to determines the features (in this case the clades at each taxonomic rank) most likely to explain the differences between Parkinson’s patients and control subjects

Set the full-path to LEfSe folder

Set the full-paths to LEfSe’s python scripts if they are not in PATH

lefse-format_input = "/path-to-lefse-format_input.py"   # e.g. /home/ngs/lefse/lefse-format_input.py
run_lefse = "/path-to-run_lefse.py"         # e.g. /home/ngs/lefse/run_lefse.py

Run LEfSe

Run LEfSe in terminal/console

/path-to-lefse-format_input.py lefse/expr1.lefse_table.txt lefse/expr1.lefse_table.in -c 1 -o 1000000
/path-to-lefse-format_input.py lefse/expr2.lefse_table.txt lefse/expr2.lefse_table.in -c 1 -o 1000000

/path-to-run_lefse.py lefse/expr1.lefse_table.in lefse/expr1.lefse_table.res -b 100 -a 1 -l 1
/path-to-run_lefse.py lefse/expr2.lefse_table.in lefse/expr2.lefse_table.res -b 100 -a 1 -l 1

Run LEfSe in R

system2(lefse_format_input, 
    args = c("lefse/expr1.lefse_table.txt", "lefse/expr1.lefse_table.in", 
         "-c 1", "-o 1000000"))

system2(lefse_format_input, 
    args = c("lefse/expr2.lefse_table.txt", "lefse/expr2.lefse_table.in", 
         "-c 1", "-o 1000000"))

# set Kruskal-Wallis alpha (-a) to 1 to allow returning of all P-value to perform adjustment later
system2(run_lefse, 
    args = c("lefse/expr1.lefse_table.in", "lefse/expr1.lefse_table.res", 
         "-b 100", "-a 1", "-l 1"), stdout = TRUE)

## [1] "Number of significantly discriminative features: 36 ( 269 ) before internal wilcoxon"
## [2] "Number of discriminative features with abs LDA score > 1.0 : 36"

system2(run_lefse, 
    args = c("lefse/expr2.lefse_table.in", "lefse/expr2.lefse_table.res", 
         "-b 100", "-a 1", "-l 1"), stdout = TRUE)

## [1] "Number of significantly discriminative features: 11 ( 267 ) before internal wilcoxon"
## [2] "Number of discriminative features with abs LDA score > 1.0 : 9"

Multiple testing correction

Perform multiple testing correction on reported P-values

# set fdr threshold at 0.1
q = 0.1

expr1 = data.frame(fread("lefse/expr1.lefse_table.res"))
expr1 = pcorrection(expr1, q)

expr2 = data.frame(fread("lefse/expr2.lefse_table.res"))
expr2 = pcorrection(expr2, q)

# You can use table() to check number of taxa pass the threshold
table(expr1$V3 != "")

## 
## FALSE  TRUE 
##   234    36

table(expr2$V3 != "")

## 
## FALSE  TRUE 
##   261     9

# Write new result file with the P-value column replaced by the FDR
write.table(expr1[,c(1:4,6)], file = "lefse/expr1.lefse_table.res.padj", sep = "\t", quote = F, 
        row.names = F, col.names = F)
write.table(expr2[,c(1:4,6)], file = "lefse/expr2.lefse_table.res.padj", sep = "\t", quote = F, 
        row.names = F, col.names = F)

Plot cladogram

Note: export2graphlan and GraPhlAn can only run in Python 2.7. They do not support Python 3.

Note: Update the coloring lefse/expr1.colors and lefse/expr2.colors if necessary. The colors should be defined in HSV (hue, saturation, value) scale.

Regarding the color file format:

The spacing between the group name and HSV scales should be tab, not space
The spacing between the HSV scales should be space not tab
There should be no new line at the end of the file

A       150., 100., 53.
B       217., 100., 91.
C       4., 85., 84.
D       39., 100., 100.

*Activate conda environment

If you have used conda to create a Python 2.7-specific environment for export2graphlan and GraPhlAn, for example an environment called graphlan, you can activate it now to run both programs.

Note: Use conda deactivate to deactivate from a currently active conda environment if required.

conda activate graphlan

Run `export2graphlan` and `GraPhlAn` in terminal/console

Note: Remember to change the PATHs to each script accordingly.

Patient vs. control at baseline (expr1)

# Run export2graphlan
/path-to-export2graphlan.py -i lefse/expr1.lefse_table.txt -o lefse/expr1.lefse_table.res.padj \
-t lefse/expr1.graphlan_tree.txt -a lefse/expr1.graphlan_annot.txt --external_annotations 2,3,4,5,6 \
--fname_row 0 --biomarkers2colors lefse/expr1.colors

# Run graphlan
/path-to-graphlan_annotate.py --annot lefse/expr1.graphlan_annot.txt lefse/expr1.graphlan_tree.txt lefse/expr1.graphlan_outtree.txt

# Convert 'lsqb' and 'rsqb' back to square bracket symbols 
sed 's/lsqb/[/' lefse/expr1.graphlan_outtree.txt | sed 's/rsqb/]/' > lefse/expr1.graphlan_outtree_sqb.txt
sed 's/lsqb/[/' lefse/expr1.lefse_table.res.padj | sed 's/rsqb/]/' > lefse/expr1.lefse_table.res_sqb.padj

# Create cladogram
/path-to-graphlan.py --dpi 150 lefse/expr1.graphlan_outtree_sqb.txt lefse/expr1.graphlan.png --external_legends --size 8 --pad 0.2

# Convert to readble output
perl /path-to-lefse.pl lefse/expr1.lefse_table.res_sqb.padj lefse/expr1.graphlan_outtree_sqb.txt lefse/expr1.out

“lefse/expr1.graphlan.png”

Patient at baseline vs. patient at followup (expr2)

# Run export2graphlan
/path-to-export2graphlan.py -i lefse/expr2.lefse_table.txt -o lefse/expr2.lefse_table.res.padj \
-t lefse/expr2.graphlan_tree.txt -a lefse/expr2.graphlan_annot.txt --external_annotations 2,3,4,5,6 \
--fname_row 0 --biomarkers2colors lefse/expr2.colors

# Run graphlan
/path-to-graphlan_annotate.py --annot lefse/expr2.graphlan_annot.txt lefse/expr2.graphlan_tree.txt lefse/expr2.graphlan_outtree.txt

# Convert 'lsqb' and 'rsqb' back to square bracket symbols
sed 's/lsqb/[/' lefse/expr2.graphlan_outtree.txt | sed 's/rsqb/]/' > lefse/expr2.graphlan_outtree_sqb.txt
sed 's/lsqb/[/' lefse/expr2.lefse_table.res.padj | sed 's/rsqb/]/' > lefse/expr2.lefse_table.res_sqb.padj

# Create cladogram
/path-to-graphlan.py --dpi 150 lefse/expr2.graphlan_outtree_sqb.txt lefse/expr2.graphlan.png --external_legends --size 8 --pad 0.2

# Convert to readble output
perl /path-to-lefse.pl lefse/expr2.lefse_table.res_sqb.padj lefse/expr2.graphlan_outtree_sqb.txt lefse/expr2.out

“lefse/expr2.graphlan.png”

Note: Use conda deactivate to deactivate the graphlan environment if required.

Plot results

Load parsed results

Patient vs. control at baseline (expr1)

res1 = data.frame(fread("lefse/expr1.out"))
res1 = res1[order(res1$order),]
names(res1)[c(5:7)] = c("Group","LDA","FDR")
res1$taxon = paste0(res1$taxon, "(",res1$rank, ";", " ",res1$order, ")")
res1$taxon = as.factor(res1$taxon)
res1$taxon = factor(res1$taxon, levels = unique(res1$taxon[order(res1$order, decreasing = T)]))
res1$Group = as.factor(res1$Group)
levels(res1$taxon) = gsub("Escherichia_Shigella","Escherichia/Shigella",levels(res1$taxon))
levels(res1$taxon) = gsub("_UCG_","_UCG-",levels(res1$taxon))

Patient at baseline vs. patient at followup (expr2)

res2 = data.frame(fread("lefse/expr2.out"))
res2 = res2[order(res2$order),]
names(res2)[c(5:7)] = c("Group","LDA","FDR")
res2$taxon = paste0(res2$taxon, "(",res2$rank, ";", " ",res2$order, ")")
res2$taxon = as.factor(res2$taxon)
res2$taxon = factor(res2$taxon, levels = unique(res2$taxon[order(res2$order, decreasing = T)]))
res2$Group = as.factor(res2$Group)
levels(res2$taxon) = gsub("Escherichia_Shigella","Escherichia/Shigella",levels(res2$taxon))
levels(res2$taxon) = gsub("_UCG_","_UCG-",levels(res2$taxon))

Plot bar

Patient vs. control at baseline (expr1)

png("lefse/expr1.lefse_table.png", width = 6, height = 5, units = "in", res = 300)
ggplot(res1, aes(taxon, LDA, fill = Group)) + geom_bar(stat = "identity", width = 0.7, size = 0.5) + 
    coord_flip() + theme_bw() + facet_grid(rows = vars(Group), scales = "free_y", space = "free_y") +
    scale_fill_manual(values = c("control" = "blue", "patient" = "red")) +
    theme(legend.position = "none", 
          axis.text.y = element_text(face = "bold", size = 8), 
          strip.placement = "outside",
          strip.text.y = element_text(face = "bold", angle = 0)) +
    labs(title = "LEfSe Analysis (Baseline)", x = "Taxon", y = "LDA")
invisible(dev.off())

“lefse/expr1.lefse_table.png”

Patient at baseline vs. patient at followup (expr2)

png("lefse/expr2.lefse_table.png", width = 6, height = 3, units = "in", res = 300)
ggplot(res2, aes(taxon, LDA, fill = Group)) + geom_bar(stat = "identity", width = 0.7, size = 0.5) + 
    coord_flip() + theme_bw() + facet_grid(rows = vars(Group), scales = "free_y", space = "free_y") +
        scale_fill_manual(values = c("baseline" = "blue", "followup" = "red")) +
        theme(legend.position = "none", 
              axis.text.y = element_text(face = "bold", size = 8), 
          strip.placement = "outside",
              strip.text.y = element_text(face = "bold", angle = 0)) +
        labs(title = "LEfSe Analysis (followup)", x = "Taxon", y = "LDA")
invisible(dev.off())

“lefse/expr2.lefse_table.png”

Save current workspace

save.image(file = "3_lefse_tutorial.RData")

Session information

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /home/ihsuan/miniconda3_4.10.0/envs/r4-16S/lib/libopenblasp-r0.3.15.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8       
##  [4] LC_COLLATE=C               LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.7       purrr_0.3.4      
##  [5] readr_1.4.0       tidyr_1.1.3       tibble_3.1.2      tidyverse_1.3.1  
##  [9] ggplot2_3.3.4     phyloseq_1.36.0   data.table_1.14.0 knitr_1.33       
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-2            hwriter_1.3.2               ellipsis_0.3.2             
##   [4] XVector_0.32.0              GenomicRanges_1.44.0        fs_1.5.0                   
##   [7] rstudioapi_0.13             bit64_4.0.5                 AnnotationDbi_1.54.0       
##  [10] fansi_0.4.2                 lubridate_1.7.10            xml2_1.3.2                 
##  [13] codetools_0.2-18            splines_4.1.0               cachem_1.0.5               
##  [16] geneplotter_1.70.0          dada2_1.20.0                ade4_1.7-17                
##  [19] jsonlite_1.7.2              Rsamtools_2.8.0             annotate_1.70.0            
##  [22] broom_0.7.8                 cluster_2.1.2               dbplyr_2.1.1               
##  [25] png_0.1-7                   compiler_4.1.0              httr_1.4.2                 
##  [28] backports_1.2.1             fastmap_1.1.0               assertthat_0.2.1           
##  [31] Matrix_1.3-4                cli_2.5.0                   htmltools_0.5.1.1          
##  [34] tools_4.1.0                 igraph_1.2.6                gtable_0.3.0               
##  [37] glue_1.4.2                  GenomeInfoDbData_1.2.6      reshape2_1.4.4             
##  [40] ShortRead_1.50.0            Rcpp_1.0.6                  Biobase_2.52.0             
##  [43] cellranger_1.1.0            vctrs_0.3.8                 Biostrings_2.60.0          
##  [46] rhdf5filters_1.4.0          multtest_2.48.0             ape_5.5                    
##  [49] nlme_3.1-152                iterators_1.0.13            xfun_0.24                  
##  [52] ps_1.6.0                    rvest_1.0.0                 lifecycle_1.0.0            
##  [55] XML_3.99-0.6                zlibbioc_1.38.0             MASS_7.3-54                
##  [58] scales_1.1.1                hms_1.1.0                   MatrixGenerics_1.4.0       
##  [61] parallel_4.1.0              SummarizedExperiment_1.22.0 biomformat_1.20.0          
##  [64] rhdf5_2.36.0                RColorBrewer_1.1-2          yaml_2.2.1                 
##  [67] memoise_2.0.0               RSQLite_2.2.5               latticeExtra_0.6-29        
##  [70] stringi_1.6.2               highr_0.9                   genefilter_1.74.0          
##  [73] S4Vectors_0.30.0            foreach_1.5.1               permute_0.9-5              
##  [76] BiocGenerics_0.38.0         BiocParallel_1.26.0         GenomeInfoDb_1.28.0        
##  [79] rlang_0.4.11                pkgconfig_2.0.3             bitops_1.0-7               
##  [82] matrixStats_0.59.0          evaluate_0.14               lattice_0.20-44            
##  [85] Rhdf5lib_1.14.0             GenomicAlignments_1.28.0    bit_4.0.4                  
##  [88] tidyselect_1.1.1            DESeq2_1.32.0               plyr_1.8.6                 
##  [91] magrittr_2.0.1              R6_2.5.0                    IRanges_2.26.0             
##  [94] generics_0.1.0              DelayedArray_0.18.0         DBI_1.1.1                  
##  [97] pillar_1.6.1                haven_2.4.1                 withr_2.4.2                
## [100] mgcv_1.8-36                 KEGGREST_1.32.0             survival_3.2-11            
## [103] RCurl_1.98-1.3              modelr_0.1.8                crayon_1.4.1               
## [106] utf8_1.2.1                  rmarkdown_2.9               jpeg_0.1-8.1               
## [109] locfit_1.5-9.4              readxl_1.3.1                blob_1.2.1                 
## [112] vegan_2.5-7                 reprex_2.0.0                digest_0.6.27              
## [115] xtable_1.8-4                RcppParallel_5.1.4          stats4_4.1.0               
## [118] munsell_0.5.0

16S rDNA amplicon sequencing analysis using R (Part 3: LEfSe & GraPhlAn)

I-Hsuan Lin

July 19, 2021

Workflow (Continues from Part 2)

Start `R`

Load package and set path

Load saved workspace from Part 2

Prepare input data

Subset samples

Prepare LEfSe input

Output to file

Perform LEfSe analysis

Set the full-path to LEfSe folder

Run LEfSe

Run LEfSe in terminal/console

Run LEfSe in R

Multiple testing correction

Plot cladogram

*Activate conda environment

Run `export2graphlan` and `GraPhlAn` in terminal/console

Plot results

Load parsed results

Plot bar

Save current workspace

Session information

16S rDNA amplicon sequencing analysis using R (Part 3: LEfSe & GraPhlAn)

I-Hsuan Lin

July 19, 2021

Workflow (Continues from Part 2)

Start R

Load package and set path

Load saved workspace from Part 2

Prepare input data

Subset samples

Prepare LEfSe input

Output to file

Perform LEfSe analysis

Set the full-path to LEfSe folder

Run LEfSe

Run LEfSe in terminal/console

Run LEfSe in R

Multiple testing correction

Plot cladogram

*Activate conda environment

Run export2graphlan and GraPhlAn in terminal/console

Plot results

Load parsed results

Plot bar

Save current workspace

Session information

Start `R`

Run `export2graphlan` and `GraPhlAn` in terminal/console