Tessa Pierce Ward
    • Create new note
    • Create a note from template
      • Sharing URL Link copied
      • /edit
      • View mode
        • Edit mode
        • View mode
        • Book mode
        • Slide mode
        Edit mode View mode Book mode Slide mode
      • Customize slides
      • Note Permission
      • Read
        • Only me
        • Signed-in users
        • Everyone
        Only me Signed-in users Everyone
      • Write
        • Only me
        • Signed-in users
        • Everyone
        Only me Signed-in users Everyone
      • Engagement control Commenting, Suggest edit, Emoji Reply
    • Invite by email
      Invitee

      This note has no invitees

    • Publish Note

      Share your work with the world Congratulations! 🎉 Your note is out in the world Publish Note No publishing access yet

      Your note will be visible on your profile and discoverable by anyone.
      Your note is now live.
      This note is visible on your profile and discoverable online.
      Everyone on the web can find and read all notes of this public team.

      Your account was recently created. Publishing will be available soon, allowing you to share notes on your public page and in search results.

      Your team account was recently created. Publishing will be available soon, allowing you to share notes on your public page and in search results.

      Explore these features while you wait
      Complete general settings
      Bookmark and like published notes
      Write a few more notes
      Complete general settings
      Write a few more notes
      See published notes
      Unpublish note
      Please check the box to agree to the Community Guidelines.
      View profile
    • Commenting
      Permission
      Disabled Forbidden Owners Signed-in users Everyone
    • Enable
    • Permission
      • Forbidden
      • Owners
      • Signed-in users
      • Everyone
    • Suggest edit
      Permission
      Disabled Forbidden Owners Signed-in users Everyone
    • Enable
    • Permission
      • Forbidden
      • Owners
      • Signed-in users
    • Emoji Reply
    • Enable
    • Versions and GitHub Sync
    • Note settings
    • Note Insights New
    • Engagement control
    • Make a copy
    • Transfer ownership
    • Delete this note
    • Save as template
    • Insert from template
    • Import from
      • Dropbox
      • Google Drive
      • Gist
      • Clipboard
    • Export to
      • Dropbox
      • Google Drive
      • Gist
    • Download
      • Markdown
      • HTML
      • Raw HTML
Menu Note settings Note Insights Versions and GitHub Sync Sharing URL Create Help
Create Create new note Create a note from template
Menu
Options
Engagement control Make a copy Transfer ownership Delete this note
Import from
Dropbox Google Drive Gist Clipboard
Export to
Dropbox Google Drive Gist
Download
Markdown HTML Raw HTML
Back
Sharing URL Link copied
/edit
View mode
  • Edit mode
  • View mode
  • Book mode
  • Slide mode
Edit mode View mode Book mode Slide mode
Customize slides
Note Permission
Read
Only me
  • Only me
  • Signed-in users
  • Everyone
Only me Signed-in users Everyone
Write
Only me
  • Only me
  • Signed-in users
  • Everyone
Only me Signed-in users Everyone
Engagement control Commenting, Suggest edit, Emoji Reply
  • Invite by email
    Invitee

    This note has no invitees

  • Publish Note

    Share your work with the world Congratulations! 🎉 Your note is out in the world Publish Note No publishing access yet

    Your note will be visible on your profile and discoverable by anyone.
    Your note is now live.
    This note is visible on your profile and discoverable online.
    Everyone on the web can find and read all notes of this public team.

    Your account was recently created. Publishing will be available soon, allowing you to share notes on your public page and in search results.

    Your team account was recently created. Publishing will be available soon, allowing you to share notes on your public page and in search results.

    Explore these features while you wait
    Complete general settings
    Bookmark and like published notes
    Write a few more notes
    Complete general settings
    Write a few more notes
    See published notes
    Unpublish note
    Please check the box to agree to the Community Guidelines.
    View profile
    Engagement control
    Commenting
    Permission
    Disabled Forbidden Owners Signed-in users Everyone
    Enable
    Permission
    • Forbidden
    • Owners
    • Signed-in users
    • Everyone
    Suggest edit
    Permission
    Disabled Forbidden Owners Signed-in users Everyone
    Enable
    Permission
    • Forbidden
    • Owners
    • Signed-in users
    Emoji Reply
    Enable
    Import from Dropbox Google Drive Gist Clipboard
       Owned this note    Owned this note      
    Published Linked with GitHub
    • Any changes
      Be notified of any changes
    • Mention me
      Be notified of mention me
    • Unsubscribe
    # K-mer-based Taxonomic Profiling of an Environmental Metagenome [![hackmd-github-sync-badge](https://hackmd.io/0ErEGoopTHaSwyRVoDovdg/badge)](https://hackmd.io/0ErEGoopTHaSwyRVoDovdg) --- ## Setup In this session, we'll use sourmash to analyze the composition of a metagenome, both genomically and taxonomically. We'll also use sourmash to classify some MAGs and integrate them into our analysis. ### Install sourmash First, we need to install the software! We'll use conda/mamba to do this. The below command installs [sourmash](http://sourmash.readthedocs.io/) and [GNU parallel](https://www.gnu.org/software/parallel/). Run: ``` # create a new environment mamba create -n smash -y -c conda-forge -c bioconda sourmash parallel ``` to install the software, and then run ``` conda activate smash ``` to activate the conda environment so you can run the software. :::success Victory conditions: - your prompt should start with `(smash) ` - you should now be able to run `sourmash` and have it output usage information!! ::: ### Create a working subdirectory Make a directory named `kmers` and change into it. ``` mkdir -p ~/kmers cd ~/kmers ``` ### Get the reference and metagenome data :::info **We've prepared your remote computers with all data so we can avoid issues with download time.** If you're using this tutorial outside of STAMPS, there is download information at the bottom of the tutorial. ::: Copy the data to your `kmers` folder ``` cp -r /opt/shared/kmers/smash-data/* ./ ``` > This will copy two folders: > - `reference` > - `metagenome` > These contain the important files we'll be using for this tutorial. ### The Lemonade Metagenome Dataset We're going to start by doing a reference-based _compositional analysis_ of a metagenome sample sequenced by the microbial diversity course. Paper: [Microbial community dynamics and coexistence in a sulfide-driven phototrophic bloom](https://environmentalmicrobiome.biomedcentral.com/articles/10.1186/s40793-019-0348-0) ![](https://hackmd.io/_uploads/H1kvgHhq3.png) > The samples come from the Trunk River lagoon in Falmouth (just a short ways up the road!) which has large amounts of biomass and which experiences periodic bright yellow microbial blooms, typically associated with disturbances to the ecosystem (source). These metagenomes were sampled from human-made disturbances that caused blooms. The accession numbers are SRR8859675 and SRR8859678. We'll use sample SRR8859675 for today, and you can view sample info [here](https://www.ebi.ac.uk/ena/browser/view/SRR8859675?show=reads) on the ENA. ### Prepare sample reads First, we're going to prepare the metagenome for use with sourmash by converting it into a _signature file_ containing _sketches_ of the k-mers in the metagenome. This is the step that "shreds" all of the reads into k-mers of size 31, and then does further data reduction by [sketching](https://en.wikipedia.org/wiki/Streaming_algorithm) the resulting k-mers. To build a signature file, we run `sourmash sketch dna` like so: ``` sourmash sketch dna -p k=31,abund metagenome/SRR8859675*.gz \ -o SRR8859675.sig.gz --name SRR8859675 ``` Here we're telling sourmash to sketch at k=31, and to track k-mer multiplicity (with 'abund'). We sketch _both_ metagenome files together into a single signature named `SRR8859675` and stored in the file `SRR8859675.sig.gz`. When we run this, we should see: >`calculated 1 signature for 3452142 sequences taken from 2 files` which tells you how many reads there are in these two files! If you look at the resulting files, ``` ls -lh SRR8859675* ``` you'll see that the signature file is _much_ smaller (2.5mb) than the metagenome files (~600mb). This is because of the way sourmash uses a reduced representation of the data, and it's what makes sourmash fast. Please see the paper above for more info! ### Look at Reference Database For this purpose, we're need a database of known genomes. We'll use the GTDB genomic representatives database, containing ~85,000 genomes - that's because it's smaller than the full GTDB database (~403,000) or Genbank (~1.3m), and hence faster. But you can download and use those on your own, if you like! You can see all sourmash prepared reference databases on [the sourmash prepared databases page](https://sourmash.readthedocs.io/en/latest/databases.html). This file is 2.2 GB: ``` ls -lh reference/gtdb-rs214-reps.k31.zip ``` and you can examine the contents with sourmash `sig summarize`: ``` sourmash sig summarize reference/gtdb-rs214-reps.k31.zip ``` which will show you: ``` ** loading from 'reference/gtdb-rs214-reps.k31.zip' path filetype: ZipFileLinearIndex location: /home/stamps/kmers/smash-data/reference/gtdb-rs214-reps.k31.zip is database? yes has manifest? yes num signatures: 85205 ** examining manifest... total hashes: 271122208 summary of sketches: 85205 sketches with DNA, k=31, scaled=1000 271122208 total hashes ``` There's a lot of things to digest in this output but the two main ones are: * there are 85,205 genome sketches in this database, for a k-mer size of 31 * this database represents 271 *billion* k-mers (multiply number of hashes by the scaled number) Note that the GTDB prepared database we downloaded above was built using the same `sourmash sketch dna` command we used above, but applied to 85,000 genomes and stored in a zip file. The associated taxonomy spreadsheet contains information connecting Genbank Assembly dataset identifiers (e.g. GCA_000008085.1) to the GTDB taxonomy - take a look: ``` head -2 reference/gtdb-rs214.lineages.csv ``` will show you: ``` ident,gtdb_representative,superkingdom,phylum,class,order,family,genus,species GCA_000008085.1,t,d__Archaea,p__Nanoarchaeota,c__Nanoarchaeia,o__Nanoarchaeales,f__Nanoarchaeaceae,g__Nanoarchaeum,s__Nanoarchaeum equitans ``` Let's prepare the taxonomy database for faster access: ``` sourmash tax prepare --taxonomy reference/gtdb-rs214.lineages.csv \ -o gtdb-rs214.taxonomy.sqldb -F sql ``` This creates a file `gtdb-rs214.taxonomy.sqldb` that contains all the information in the CSV file, but which is faster to load than the CSV file. ## Find matching reference genomes with `sourmash gather` At last, we have the ingredients we need to analyze the metagenome against GTDB! * the software is installed * the GTDB database is downloaded * the metagenome is downloaded and sketched Now, we'll run the [sourmash gather](https://sourmash.readthedocs.io/en/latest/command-line.html#sourmash-gather-find-metagenome-members) command to find matching genomes. Run gather - this will take ~6 minutes: ``` sourmash gather SRR8859675.sig.gz reference/gtdb-rs214-reps.k31.zip --save-matches matches.zip ``` Here we are saving the matching genome sketches to `matches.zip` so we can rerun the analysis if we like. The results will look like this: ``` Prefetch found 38 signatures with overlap >= 50.0 kbp. Doing gather to generate minimum metagenome cover. overlap p_query p_match avg_abund --------- ------- ------- --------- 2.0 Mbp 0.4% 31.8% 1.3 GCF_004138165.1 Candidatus Chloroploca sp. Khr17 strain=Khr17, ASM413816v1 1.9 Mbp 0.5% 66.9% 2.1 GCF_900101955.1 Desulfuromonas thiophila strain=DSM 8987, IMG-taxon 2602042025 annotated assembly 0.6 Mbp 0.3% 23.3% 3.2 GCA_016938795.1 Chromatiaceae bacterium, ASM1693879v1 0.6 Mbp 0.5% 27.3% 6.6 GCA_016931495.1 Chlorobiaceae bacterium, ASM1693149v1 found 24 matches total; the recovered matches hit 5.3% of the abundance-weighted query. the recovered matches hit 2.4% of the query k-mers (unweighted). ``` In this output: * the last column is the name of the matching GTDB genome * the first column is the estimated overlap between the metagenome and that genome, in base pairs (estimated from shared k-mers) * the second column, `p_query` is the percentage of metagenome k-mers (weighted by multiplicity) that match to the genome; this will approximate the percentage of _metagenome reads_ that will map to this genome, if you map. * the third column, `p_match`, is the percentage of the genome k-mers that are matched by the metagenome; this will approximate the percentage of _genome bases_ that will be covered by mapped reads; * the fourth column is the estimated mean abundance of this genome in the metagenome. The other interesting number is here: >`the recovered matches hit 5.3% of the abundance-weighted query` which tells you that you should expect about 5.3% of the metagenome reads to map to these 24 reference genomes. :::info ## Weighted vs Unweighted Let's look at the summary lines at the bottom of the `gather` results: ``` found 24 matches total; the recovered matches hit 5.3% of the abundance-weighted query. the recovered matches hit 2.4% of the query k-mers (unweighted). ``` The `5.3% (abundance-weighted)` is the percentage of the metagenome k-mers that are matched by all reference genomes, meaning ~95% of the dataset does not match any reference. The `2.4% (unweighted)` represents the proportion of _unique_ kmers in the metagenome that are found in any genome. This will approximate the percentage of _genome bases_ that will be covered by mapped reads The `unweighted` percent is (approximately) the following: * suppose you assembled the entire metagenome perfectly into perfect contigs (**note, this is impossible, although you can get close with "unitigs"**); * and then matched all the genomes to the contigs; * approximately 2.5% of the bases in the contigs would have genomes that match to them. Interestingly, this is the _only_ number in this entire tutorial that is essentially impossible to estimate any way other than with k-mers. This number is also a big underestimate of the "true" number for the metagenome - we'll discuss more later :) ::: ## Build a taxonomic profile of the metagenome We can use these matching genomes to build a taxonomic profile of the metagenome using [sourmash tax metagenome](https://sourmash.readthedocs.io/en/latest/command-line.html#sourmash-tax-subcommands-for-integrating-taxonomic-information-into-gather-results). This method has been shown to perform quite well in benchmarking using mock metagenome datasets. Read more about it for short and long reads in [Portik et. al 2022](https://link.springer.com/article/10.1186/s12859-022-05103-0). Let's re-run gather so we can save CSV output. This run will be much faster, because instead of using the entire database, we'll used our pre-saved 'matches.zip' which contains all 24 genomes with matches to the metagenome. ``` # run gather, save the results to a CSV sourmash gather SRR8859675.sig.gz matches.zip -o SRR8859675.x.gtdb.csv # use tax metagenome to classify the metagenome sourmash tax metagenome --gather-csv SRR8859675.x.gtdb.csv \ --taxonomy gtdb-rs214.taxonomy.sqldb \ --output-format human --rank order ``` this shows you the rank, taxonomic lineage, and weighted fraction of the metagenome at the 'order' rank. ## Interlude: why reference-based analyses are problematic for environmental metagenomes Reference-based metagenome classification is highly dependent on the organisms present in our reference databases. For well-studied environments, such as human-associated microbiomes, your classification percentage is likely to be quite high. In contrast, this is an environmental metagenome, and you can see that we're estimating only 5.3% of it will map to GTDB reference genomes! Wow, that's **terrible**! Our taxonomic and/or functional analysis will be based on only 1/20th of the data! What could we do to improve that?? There are two basic options - (1) Use a more complete reference database, like the entire GTDB, or Genbank. This will only get you so far, unfortunately. (See exercises at end.) (2) Assemble and bin the metagenome to produce new reference genomes! There are other things you could think about doing here, too, but these are probably the "easiest" options. Last year at STAMPS we assembled these data into Metagenome-Assembled Genomes (MAGs) with ATLAS as part of [Taylor Reiter's STAMPS 2022 tutorial on assembly and binning](https://github.com/mblstamps/stamps2022/blob/main/assembly_and_binning/tutorial_assembly_and_binning.md). Let's add the MAGs generated there to our analysis. We'll need to: * sketch the MAGs with k=31; * re-run sourmash gather with both GTDB _and_ the MAGs. Let's do it!! ## Use MAGs as additional "reference" genomes Three MAGs were assembled and binned from the metagenomes. Let's see if including data from these metagenome-assembled genomes improves the classifable fraction of the metagenome. Let's sketch the MAGs: ``` sourmash sketch dna metagenome/MAG*.fasta --name-from-first ``` here, `--name-from-first` is a convenient way to give them distinguishing names based on the name of the first contig in the FASTA file; you can see the names of the signatures by doing: ``` sourmash sig describe MAG1.fasta.sig ``` Now, let's re-do the metagenome classification with the MAGs: ``` sourmash gather SRR8859675.sig.gz MAG*.sig matches.zip -o SRR8859675.x.gtdb+MAGS.csv ``` and look, we classify a lot more! ``` overlap p_query p_match avg_abund --------- ------- ------- --------- 2.3 Mbp 12.1% 99.9% 39.4 MAG2_1 2.2 Mbp 26.5% 99.9% 92.4 MAG3_1 2.0 Mbp 0.4% 31.8% 1.3 GCF_004138165.1 Candidatus Chloroploca sp. Khr17 strain=Khr17, ASM413816v1 1.9 Mbp 0.5% 66.9% 2.1 GCF_900101955.1 Desulfuromonas thiophila strain=DSM 8987, IMG-taxon 2602042025 annotated assembly 1.0 Mbp 2.7% 100.0% 20.3 MAG1_1 0.6 Mbp 0.3% 23.2% 3.1 GCA_016938795.1 Chromatiaceae bacterium, ASM1693879v1 0.6 Mbp 0.1% 24.5% 2.1 GCA_016931495.1 Chlorobiaceae bacterium, ASM1693149v1 ... found 26 matches total; the recovered matches hit 43.5% of the abundance-weighted query. the recovered matches hit 4.0% of the query k-mers (unweighted). ``` Here we see a few interesting things - (1) The three MAG matches are all ~100% present in the metagenome. (2) They are all at high abundance in the metagenome, because assembly needs genomes to be ~5x or more in abundance in order to work! (3) Because they're at high abundance and 100% present, they account for _a lot_ of the metagenome! What's the remaining 50%? There are several answers - (1) most of the genomes actually the metagenome aren't in our reference database; (2) not everything in the metagenome is high enough coverage to bin into MAGs; (3) not everything in the metagenome is bacterial or archaeal, and we didn't do viral or eukaryotic binning; (4) some of what's in the metagenome k-mers may simply be erroneous (although with abundance weighting, this is likely to be only a small chunk of things) ## Add taxonomic information from these MAGs. We can infer taxonomic information for the MAGs in a similar way as we did for the metagenome, above. ### Classify the MAGs First, classify the genomes using GTDB; this will use trace overlaps between contigs in the MAGs and GTDB genomes to find the reference genome that is most similar to the entire "MAG" (assembled contig bin). ``` for i in MAG*.fasta.sig do # get 'MAG' prefix. => NAME NAME=$(basename $i .fasta.sig) # search against GTDB echo sourmash gather $i reference/gtdb-rs214-reps.k31.zip \ --threshold-bp=5000 \ -o ${NAME}.x.gtdb.csv done | parallel ``` (This will take about a minute.) Here, we're using a for loop and [GNU parallel](https://www.gnu.org/software/parallel/) to classify the three genomes in parallel. If you scan the results quickly, you'll see that one MAG has matches in genus Prosthecochloris, another MAG has matches to Chlorobaculum, and one has matches to Candidatus Moranbacteria. Let's classify them "officially" using sourmash and an average nucleotide identity threshold of 0.8 - ``` sourmash tax genome --gather-csv MAG*.x.gtdb.csv \ --taxonomy gtdb-rs214.taxonomy.sqldb \ --output-format human \ --ani 0.8 ``` This is an extremely liberal ANI threshold, incidentally; in reality you'd probably want to do something more stringent, as at least one of these is probably a new species. You should see: ``` >sample name proportion lineage >----------- ---------- ------- >MAG3_1 5.3% d__Bacteria;p__Bacteroidota;c__Chlorobia;o__Chlorobiales;f__Chlorobiaceae;g__Prosthecochloris;s__Prosthecochloris vibrioformis >MAG2_1 5.0% d__Bacteria;p__Bacteroidota;c__Chlorobia;o__Chlorobiales;f__Chlorobiaceae;g__Chlorobaculum;s__Chlorobaculum parvum_B >MAG1_1 1.1% d__Bacteria;p__Patescibacteria;c__Paceibacteria;o__Moranbacterales;f__UBA1568;g__JAAXTX01;s__JAAXTX01 sp013334245 ``` The proportion here is the fraction of k-mers in the MAG that are annotated. Now let's turn this into a lineage spreadsheet: ``` sourmash tax genome --gather-csv MAG*.x.gtdb.csv \ --taxonomy gtdb-rs214.taxonomy.sqldb --output-format lineage_csv \ --ani 0.8 -o MAGs ``` This will produce a file `MAGs.lineage.csv`; let's take a look: ``` cat MAGs.lineage.csv ``` You should see: ``` >ident,superkingdom,phylum,class,order,family,genus,species >MAG1_1,d__Bacteria,p__Patescibacteria,c__Paceibacteria,o__Moranbacterales,f__UBA1 568,g__JAAXTX01,s__JAAXTX01 sp013334245 >MAG2_1,d__Bacteria,p__Bacteroidota,c__Chlorobia,o__Chlorobiales,f__Chlorobiaceae, g__Chlorobaculum,s__Chlorobaculum parvum_B >MAG3_1,d__Bacteria,p__Bacteroidota,c__Chlorobia,o__Chlorobiales,f__Chlorobiaceae, g__Prosthecochloris,s__Prosthecochloris vibrioformis ``` ## Now let's re-classify the metagenome If we re-classify the metagenome using the combined "reference" of GTDB + MAGs, we see: ``` sourmash tax metagenome --gather-csv SRR8859675.x.gtdb+MAGS.csv \ --taxonomy gtdb-rs214.taxonomy.sqldb MAGs.lineage.csv \ --output-format human --rank order ``` Now only 56.5% remains unclassified, which is much better than before! For the classified fraction, we can visualize the distribution of taxa using [sourmashconsumr](https://github.com/Arcadia-Science/sourmashconsumr/): ![](https://hackmd.io/_uploads/ryx1O4T53.png) :::spoiler :::warning To build this plot, run the following in your terminal: ``` sourmash tax annotate --gather-csv SRR8859675.x.gtdb+MAGS.csv --taxonomy gtdb-rs214.taxonomy.sqldb MAGs.lineage.csv ``` This adds a 'lineages' column to the gather results without doing any taxonomic summarization. **Then go into your R Console.** Install sourmashconsumr (warning, takes several minutes) ``` install.packages("remotes") remotes::install_github("Arcadia-Science/sourmashconsumr") ``` Load packages ``` library(sourmashconsumr) library(dplyr) library(ggplot2) library(RColorBrewer) ``` Load the data ``` gather_csv <- '~/kmers/SRR8859675.x.gtdb+MAGS.with-lineages.csv' inf = sourmashconsumr::read_taxonomy_annotate(gather_csv) %>% dplyr::mutate(query_name = basename(query_filename)) ``` plot ``` plot_taxonomy_annotate_sankey(taxonomy_annotate_df = inf, tax_glom_level = "order") ``` ::: ## Summary of Taxonomic Profiling To recap, we've done the following: * analyzed a metagenome's composition against 85,000 GTDB genomes, using 31-mers; * found that a disappointingly small fraction of the metagenome can be identified this way. * incorporated MAGs built from the metagenome into this analysis, bumping up the classification rate to ~45%; * added taxonomic output to both sets of analyses. ATLAS only bins bacterial and archaeal genomes, so we wouldn't expect much in the way of viral or eukaryotic genomes to be binned. :::info ## Optional Exercise: Evaluating how much of the metagenome assembled How much of our metagenome even _assembles_? First, sketch the assembled contigs into a sourmash signature: ``` sourmash sketch dna metagenome/SRR8859675_contigs.fasta --name-from-first ``` Let's pick a few of the matching genomes out from GTDB and evaluate how many of the k-mers from that genome match to the unassembled metagenome, and then how many of them match to the assembled contigs. Now, extract one of the top gather matches to use as a query; this is "Chromatiaceae bacterium": ``` sourmash sig cat matches.zip --include GCA_016938795.1 -o GCA_016938795.sig ``` ### Evaluate containment of known genomes in reads vs assembly Now do a containment search of this genome against both the unassembled metagenome and the assembled (but unbinned) contigs - ``` sourmash search --containment GCA_016938795.sig \ SRR8859675*.sig* --threshold=0 --ignore-abund ``` We see: ``` similarity match ---------- ----- 23.3% SRR8859675 4.7% SRR8859675_0 ``` where the first match (at 23.3% containment) is to the metagenome. (You'll note this matches the % in the gather output, too.) The second match is to the assembled contigs, and it's 4.7%. That means ~19% of the k-mers that match to this GTDB genome are present in the unassembled metagenome, but are lost during the assembly process. Why? Some thoughts and answers It _could_ be that the GTDB genome is full of errors, and those errors are shared with the metagenome, and assembly is squashing those errors. Yay! But this is extremely unlikely... This GTDB genome was built and validated entirely independently from this sample... It's more likely that one of two things is happening: (1) this sample is at low abundance in the metagenome, and assembly can only recover parts of it. (2) this sample contains _several_ strain variants of this genome, and assembly is squashing the strain variation, because that's what assembly does. Note, you can try the above with another one of the top gather matches and you'll see it's *entirely* lost in the process of assembly - ``` sourmash sig cat matches.zip --include GCF_004138165.1 -o GCF_004138165.sig sourmash search --containment GCF_004138165.sig \ SRR8859675*.sig* \ --ignore-abund --threshold=0 ``` ::: ## Concluding thoughts Above, we demonstrated a _reference-based_ analysis of shotgun metagenome data using sourmash. We then _updated our references_ using the MAGs produced from assembly and binning tutorial, which increased our classification rate substantially. Last but not least, we (optionally) looked at the loss of k-mer information due to metagenome assembly. All of these results were based on 31-mer overlap and containment with k-mers! A few points: * We would have gotten slightly different results using k=21 or k=51; more of the metagenome would have been classified with k=21, while the classification results would have been more closely specific to genomes with k=51; * sourmash is a nice multitool for doing this, but you could have gotten similar results with other methods. * Next steps here could include mapping reads to the genomes we found, and/or doing functional analysis on the matching genes and genomes. * You can do downstream analyses (e.g. diff abundance) with sourmash gather -> taxonomy results. Here's a tutorial for building a **phyloseq** object: [From raw metagenome reads to phyloseq taxonomy table using sourmash gather and sourmash taxonomy](https://taylorreiter.github.io/2022-07-28-From-raw-metagenome-reads-to-phyloseq-taxonomy-table-using-sourmash-gather-and-sourmash-taxonomy/). ## References and Further Links - If you want to read more about how sourmash sketching and gather, please see [Lightweight compositional analysis of metagenomes with FracMinHash and minimum metagenome covers](https://www.biorxiv.org/content/10.1101/2022.01.11.475838v2), Irber et al., 2022. - For more on sourmash taxonomic profiling, see our paper [here](https://link.springer.com/article/10.1186/s12859-022-05103-0). - We have a prototype frontend for running `sourmash gather` against the GTDB representatives database, [https://greyhound.sourmash.bio/](https://greyhound.sourmash.bio/). Try it out! It may not work for very large files :). ## After STAMPS: downloading this test data :::spoiler ### Reference ``` mkdir -p ~/kmers/smash-data/reference cd ~/kmers/smash-data/reference ``` Download the reference database: ``` curl -JLO https://farm.cse.ucdavis.edu/~ctbrown/sourmash-db/gtdb-rs214/gtdb-rs214-reps.k31.zip ``` Download and uncompress the taxonomy spreadsheet: ``` curl -JLO https://farm.cse.ucdavis.edu/~ctbrown/sourmash-db/gtdb-rs214/gtdb-rs214.lineages.csv.gz gunzip gtdb-rs214.lineages.csv.gz ``` ``` mkdir -p ~/kmers/smash-data/metagenome cd ~/kmers/smash-data/metagenome ``` ### Metagenome Data To download the metagenome from the ENA, run: ``` wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR885/005/SRR8859675/SRR8859675_1.fastq.gz wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR885/005/SRR8859675/SRR8859675_2.fastq.gz ``` Download 3 MAGs generated by ATLAS ``` curl -JLO https://osf.io/fejps/download curl -JLO https://osf.io/jf65t/download curl -JLO https://osf.io/2a4nk/download ``` We'll use sample SRR8859675 for today, and you can view sample info [here](https://www.ebi.ac.uk/ena/browser/view/SRR8859675?show=reads) on the ENA. If you want to just start at the Optional Assembly Loss exercise, you can download the files needed for the below sourmash searches from [this link](https://github.com/mblstamps/stamps2022/raw/main/kmers_and_sourmash/assembly-loss-files.zip). :::

    Import from clipboard

    Paste your markdown or webpage here...

    Advanced permission required

    Your current role can only read. Ask the system administrator to acquire write and comment permission.

    This team is disabled

    Sorry, this team is disabled. You can't edit this note.

    This note is locked

    Sorry, only owner can edit this note.

    Reach the limit

    Sorry, you've reached the max length this note can be.
    Please reduce the content or divide it to more notes, thank you!

    Import from Gist

    Import from Snippet

    or

    Export to Snippet

    Are you sure?

    Do you really want to delete this note?
    All users will lose their connection.

    Create a note from template

    Create a note from template

    Oops...
    This template has been removed or transferred.
    Upgrade
    All
    • All
    • Team
    No template.

    Create a template

    Upgrade

    Delete template

    Do you really want to delete this template?
    Turn this template into a regular note and keep its content, versions, and comments.

    This page need refresh

    You have an incompatible client version.
    Refresh to update.
    New version available!
    See releases notes here
    Refresh to enjoy new features.
    Your user state has changed.
    Refresh to load new user state.

    Sign in

    Forgot password
    or
    Sign in via Facebook Sign in via X(Twitter) Sign in via GitHub Sign in via Dropbox Sign in with Wallet
    Wallet ( )
    Connect another wallet

    New to HackMD? Sign up

    By signing in, you agree to our terms of service.

    Help

    • English
    • 中文
    • Français
    • Deutsch
    • 日本語
    • Español
    • Català
    • Ελληνικά
    • Português
    • italiano
    • Türkçe
    • Русский
    • Nederlands
    • hrvatski jezik
    • język polski
    • Українська
    • हिन्दी
    • svenska
    • Esperanto
    • dansk

    Documents

    Help & Tutorial

    How to use Book mode

    Slide Example

    API Docs

    Edit in VSCode

    Install browser extension

    Contacts

    Feedback

    Discord

    Send us email

    Resources

    Releases

    Pricing

    Blog

    Policy

    Terms

    Privacy

    Cheatsheet

    Syntax Example Reference
    # Header Header 基本排版
    - Unordered List
    • Unordered List
    1. Ordered List
    1. Ordered List
    - [ ] Todo List
    • Todo List
    > Blockquote
    Blockquote
    **Bold font** Bold font
    *Italics font* Italics font
    ~~Strikethrough~~ Strikethrough
    19^th^ 19th
    H~2~O H2O
    ++Inserted text++ Inserted text
    ==Marked text== Marked text
    [link text](https:// "title") Link
    ![image alt](https:// "title") Image
    `Code` Code 在筆記中貼入程式碼
    ```javascript
    var i = 0;
    ```
    var i = 0;
    :smile: :smile: Emoji list
    {%youtube youtube_id %} Externals
    $L^aT_eX$ LaTeX
    :::info
    This is a alert area.
    :::

    This is a alert area.

    Versions and GitHub Sync
    Get Full History Access

    • Edit version name
    • Delete

    revision author avatar     named on  

    More Less

    Note content is identical to the latest version.
    Compare
      Choose a version
      No search result
      Version not found
    Sign in to link this note to GitHub
    Learn more
    This note is not linked with GitHub
     

    Feedback

    Submission failed, please try again

    Thanks for your support.

    On a scale of 0-10, how likely is it that you would recommend HackMD to your friends, family or business associates?

    Please give us some advice and help us improve HackMD.

     

    Thanks for your feedback

    Remove version name

    Do you want to remove this version name and description?

    Transfer ownership

    Transfer to
      Warning: is a public team. If you transfer note to this team, everyone on the web can find and read this note.

        Link with GitHub

        Please authorize HackMD on GitHub
        • Please sign in to GitHub and install the HackMD app on your GitHub repo.
        • HackMD links with GitHub through a GitHub App. You can choose which repo to install our App.
        Learn more  Sign in to GitHub

        Push the note to GitHub Push to GitHub Pull a file from GitHub

          Authorize again
         

        Choose which file to push to

        Select repo
        Refresh Authorize more repos
        Select branch
        Select file
        Select branch
        Choose version(s) to push
        • Save a new version and push
        • Choose from existing versions
        Include title and tags
        Available push count

        Pull from GitHub

         
        File from GitHub
        File from HackMD

        GitHub Link Settings

        File linked

        Linked by
        File path
        Last synced branch
        Available push count

        Danger Zone

        Unlink
        You will no longer receive notification when GitHub file changes after unlink.

        Syncing

        Push failed

        Push successfully