ycl6 / 16s-rdna-v3-v4

16S rDNA V3-V4 amplicon sequencing analysis using dada2, phyloseq, LEfSe, picrust2 and other tools. Demo: https://ycl6.github.io/16S-Demo/

Hi there 👋

Welcome to my GitHub profile!

R Packages

• scRUtils: Single-Cell Sequencing Analysis Utilities [Vignette]

NGS Tutorials

• RNA-Seq: RNA-seq Data Analysis
  • Differential gene expression (DGE) analysis
  • Differential transcript expression (DTE) analysis
  • Differential transcript usage (DTU) analysis
• Microbiome: 16S rDNA V3-V4 amplicon sequencing analysis
• Tertiary analysis
  • GO enrichment analysis with topGO, clusterProfiler & Enrichr/enrichR
  • Visualize enrichment of genomic signals on targets or regions using EnrichedHeatmap

Misc.

• Create a GitHub Pages site for your repository
• Create stacked violin plots from Seurat, SingleCellExperiment and data.frame objects
• Gallery showcasing Pandoc syntax highlight styles before and after version 2.11
• Short examples of using R in data analysis and visualisation
  • National survey
    • Deaths registered monthly in England and Wales
    • COVID-19 Vaccination Statistics in England
  • Gene expression
    • Visualise differential expression analysis results with MA plot
    • Normalisation of bulk RNA-seq counts
    • Plot ligand-receptor expression on UMAP and t-SNE projection
  • Others
    • Visualise interaction between 2 variables with parallel sets diagrams
    • Testing of the Spatial Image Analysis of Tissues (SPIAT) R package

16s-rdna-v3-v4's Issues

downstream picrust2 analysis query

Hello I-Hsuan @ycl6
Many thanks for this tutorial (https://ycl6.github.io/16S-Demo/4_picrust2_tutorial.html).
I was wondering if there is a possibility to draw a figure like this:

Your help is much appreciated.

about dada2 filterandTrim parameter

Can i get how to determine dada2 filterAndTrim() paramer?
(truncLen=c(240,230), minLen=200, maxN=0, truncQ=2, maxEE=c(2,5) ..)

For example:
Illuminia MiSeq - (2*250bp reads)
321F: TAATTCGTCCTACGGRAGGCAG
800R: AGCCGGACTACCRGGGTAT

I am curious if there is any formula.

How to best root the tree

The current practice is to randomly pick a taxa (with sample) and assign it as outgroup. There have been suggestions to do this differently. Need to re-visit this when free.

phy_tree(ps) <- root(phy_tree(ps), sample(taxa_names(ps), 1), resolve.root = TRUE)

a

He

error in evaluating argument 'x'

Hi Hsuan Lin,
Another question arose according to the script.

When I follow the line:

> new_taxtab = taxtab
> rownames(new_taxtab$tax) = df[match(rownames(new_taxtab$tax), df$sequence),]$id

> Error in h(simpleError(msg, call)) :
  error in evaluating argument 'x' when selecting method for function 'match': 'error in evaluating argument 'x' when selecting method for function 'rownames': '$ operator is invalid for atomic vectors''

How can I solve this problem?

Relative abundance

Dear,
I'm trying to reproduce the script https://ycl6.github.io/16S-Demo/1_dada2_tutorial.html and I ran into a first problem.
How can I solve it, I know there will be others but in this first one the relative abundance is only counted as 1 or 0.
Did I do any previous procedures wrong?

taxcumsum = tdt[, .N, by = TotalCounts] 
setkey(taxcumsum, TotalCounts)
taxcumsum[, CumSum := cumsum(N)]

pCumSum = ggplot(taxcumsum, aes(TotalCounts, CumSum)) + geom_point() + theme_bw() + 
  xlab("Filtering Threshold") + ylab("ASV Filtered")
gridExtra::grid.arrange(pCumSum, pCumSum + xlim(0, 500), 
                        pCumSum + xlim(0, 100), pCumSum + xlim(0, 50), nrow = 2, 
                        top = "ASVs that would be filtered vs. minimum taxa counts threshold")

mdt = fast_melt(ps)
mdt = mdt[count > 0] [!is.na(count)]
mdt[, RelativeAbundance := count / sum(count), by = taxaID] 
mdt

       Kingdom         Phylum               Class                 Order               Family
   1: Bacteria Proteobacteria Gammaproteobacteria       Methylococcales      Methylomonaceae
   2: Bacteria Proteobacteria Gammaproteobacteria Betaproteobacteriales              SC-I-84
   3: Bacteria Proteobacteria Deltaproteobacteria   Syntrophobacterales Syntrophobacteraceae
   4: Bacteria           <NA>                <NA>                  <NA>                 <NA>
   5: Bacteria Proteobacteria Gammaproteobacteria    Steroidobacterales  Steroidobacteraceae
  ---                                                                                       
1043: Bacteria Proteobacteria                <NA>                  <NA>                 <NA>
1044: Bacteria Actinobacteria      Actinobacteria                  <NA>                 <NA>
1045: Bacteria Proteobacteria Deltaproteobacteria                  <NA>                 <NA>
1046: Bacteria Proteobacteria Gammaproteobacteria    Steroidobacterales  Steroidobacteraceae
1047: Bacteria  Bacteroidetes         Bacteroidia         Bacteroidales   Prolixibacteraceae

                Genus  taxaID SampleID count RelativeAbundance
   1:      Crenothrix    OTU1      sa1   332                 1
   2:            <NA>   OTU10      sa1   147                 1
   3: Syntrophobacter  OTU10  sa1    49                 1
   4:            <NA> OTU1000      sa1     1                 1
   5:            <NA> OTU1001      sa1     1                 1
  ---                                                         
1043:            <NA>  OTU995      sa1     1                 1
1044:            <NA>  OTU996      sa1     1                 1
1045:            <NA>  OTU997      sa1     1                 1
1046:            <NA>  OTU998      sa1     1                 1
1047:            <NA>  OTU999      sa1     1                 1

Thank you very much in advance

How to work with several different paired fastq files together?

Dear Hsuan Lin,
I am using your script to process some 16S rRNA sequencing samples.
When I use it for a single sequenced sample, in a paired sample, it works.
But I would like to work on 4 paired samples together. How do I do?
I concatenated the files but it didn't work, what should I do?

The following message appeared to me:
Error in filterAndTrim(RC_fnFs, RC_filtFs, RC_fnRs, RC_filtRs, truncLen = c(250, :
Every input file must have a corresponding output file.