HackMD
  • New!
    New!  Markdown linter lands.
    Write with consistency! Enable markdown linter that hints better markdown styling you can take. Click the lightbulb on the editor status bar 💡
    Got it
      • Create new note
      • Create a note from template
    • New!  Markdown linter lands.
      New!  Markdown linter lands.
      Write with consistency! Enable markdown linter that hints better markdown styling you can take. Click the lightbulb on the editor status bar 💡
      Got it
      • Options
      • Versions and GitHub Sync
      • Transfer ownership
      • Delete this note
      • Template
      • Save as template
      • Insert from template
      • Export
      • Dropbox
      • Google Drive
      • Gist
      • Import
      • Dropbox
      • Google Drive
      • Gist
      • Clipboard
      • Download
      • Markdown
      • HTML
      • Raw HTML
      • ODF (Beta)
      • Sharing Link copied
      • /edit
      • View mode
        • Edit mode
        • View mode
        • Book mode
        • Slide mode
        Edit mode View mode Book mode Slide mode
      • Note Permission
      • Read
        • Owners
        • Signed-in users
        • Everyone
        Owners Signed-in users Everyone
      • Write
        • Owners
        • Signed-in users
        • Everyone
        Owners Signed-in users Everyone
      • More (Comment, Invitee)
      • Publishing
        Everyone on the web can find and read all notes of this public team.
        After the note is published, everyone on the web can find and read this note.
        See all published notes on profile page.
      • Commenting Enable
        Disabled Forbidden Owners Signed-in users Everyone
      • Permission
        • Forbidden
        • Owners
        • Signed-in users
        • Everyone
      • Invitee
      • No invitee
    Menu Sharing Create Help
    Create Create new note Create a note from template
    Menu
    Options
    Versions and GitHub Sync Transfer ownership Delete this note
    Export
    Dropbox Google Drive Gist
    Import
    Dropbox Google Drive Gist Clipboard
    Download
    Markdown HTML Raw HTML ODF (Beta)
    Back
    Sharing
    Sharing Link copied
    /edit
    View mode
    • Edit mode
    • View mode
    • Book mode
    • Slide mode
    Edit mode View mode Book mode Slide mode
    Note Permission
    Read
    Owners
    • Owners
    • Signed-in users
    • Everyone
    Owners Signed-in users Everyone
    Write
    Owners
    • Owners
    • Signed-in users
    • Everyone
    Owners Signed-in users Everyone
    More (Comment, Invitee)
    Publishing
    Everyone on the web can find and read all notes of this public team.
    After the note is published, everyone on the web can find and read this note.
    See all published notes on profile page.
    More (Comment, Invitee)
    Commenting Enable
    Disabled Forbidden Owners Signed-in users Everyone
    Permission
    Owners
    • Forbidden
    • Owners
    • Signed-in users
    • Everyone
    Invitee
    No invitee
       owned this note    owned this note      
    Published Linked with GitHub
    Like BookmarkBookmarked
    Subscribed
    • Any changes
      Be notified of any changes
    • Mention me
      Be notified of mention me
    • Unsubscribe
    Subscribe
    # notes from snakemake live stream C. Titus Brown Jan 20, 2021 ctbrown@ucdavis.edu [toc] ## Recording A [recording of the livestream is available.](https://video.ucdavis.edu/media/1_ygjgu1ma) ## Getting started w/RStudio in a binder I'll be using [a binder](mybinder.org/) to do everything on an anonymous remote cloud computer. [![Binder](https://binder.pangeo.io/badge_logo.svg)](https://binder.pangeo.io/v2/gh/binder-examples/r-conda/master?urlpath=rstudio) To configure my terminal in the RStudio terminal window on the binder, I'll run: ``` conda init echo "PS1='\w $ '" >> .bashrc ``` and then restart the terminal. Then I'll add two channels to conda so that I can install necessary software, ``` conda config --add channels conda-forge conda config --add channels bioconda ``` ## The commands for variant calling I'm going to start with the following minimap2-based variant pipeline that ends by producing a VCF file. ``` # download some data wget https://osf.io/4rdza/download -O SRR2584857_1.fastq.gz # download a reference genome wget https://osf.io/8sm92/download -O ecoli-rel606.fa.gz # run minimap2 minimap2 -ax sr ecoli-rel606.fa.gz SRR2584857_1.fastq.gz > SRR2584857_1.ecoli-rel606.sam # convert SAM to BAM samtools view -b -F 4 SRR2584857_1.ecoli-rel606.sam > SRR2584857_1.ecoli-rel606.bam # sort the BAM by position samtools sort SRR2584857_1.ecoli-rel606.bam > SRR2584857_1.ecoli-rel606.bam.sorted # unzip genome for bcftools gunzip -c ecoli-rel606.fa.gz > ecoli-rel606.fa # run mpileup bcftools mpileup -Ou -f ecoli-rel606.fa SRR2584857_1.ecoli-rel606.bam.sorted > SRR2584857_1.ecoli-rel606.pileup # call variants bcftools call -mv -Ob SRR2584857_1.ecoli-rel606.pileup -o SRR2584857_1.ecoli-rel606.bcf # turn the BCF into a VCF bcftools view SRR2584857_1.ecoli-rel606.bcf > SRR2584857_1.ecoli-rel606.vcf # and done! ``` ## First stopping point Snakefile env-variant-calling.yml: ```yaml channels: - bioconda - conda-forge - defaults - https://conda.anaconda.org/conda-forge dependencies: - bcftools=1.9=h68d8f2e_9 - minimap2=2.17=hed695b0_3 - samtools=1.9=h10a08f8_12 ``` Snakefile: ```python # default rule - first one in the file rule all: input: "SRR2584857_1.ecoli-rel606.vcf" rule download_data: output: "SRR2584857_1.fastq.gz" shell: """ wget https://osf.io/4rdza/download -O {output} """ rule download_reference: output: "ecoli-rel606.fa.gz" shell: """ wget https://osf.io/8sm92/download -O {output} """ rule run_minimap: input: reads="SRR2584857_1.fastq.gz", genome="ecoli-rel606.fa.gz" output: "SRR2584857_1.ecoli-rel606.sam" conda: "env-variant-calling.yml" shell: """ minimap2 -ax sr {input.genome} {input.reads} > {output} """ rule sam_to_bam: input: rules.run_minimap.output #"SRR2584857_1.ecoli-rel606.sam" output: "SRR2584857_1.ecoli-rel606.bam" conda: "env-variant-calling.yml" shell: """ samtools view -b -F 4 {input} > {output} """ rule sort_bam: output: "SRR2584857_1.ecoli-rel606.bam.sorted" input: "SRR2584857_1.ecoli-rel606.bam" conda: "env-variant-calling.yml" shell: """ samtools sort {input} > {output} """ rule call_variants: input: genome = "ecoli-rel606.fa.gz", sorted_bam="SRR2584857_1.ecoli-rel606.bam.sorted" output: ungz_genome = "ecoli-rel606.fa", pileup = "SRR2584857_1.ecoli-rel606.pileup", bcf = "SRR2584857_1.ecoli-rel606.bcf", vcf = "SRR2584857_1.ecoli-rel606.vcf" conda: "env-variant-calling.yml" shell: """ gunzip -c {input.genome} > {output.ungz_genome} bcftools mpileup -Ou -f {output.ungz_genome} {input.sorted_bam} > {output.pileup} bcftools call -mv -Ob {output.pileup} -o {output.bcf} bcftools view {output.bcf} > {output.vcf} """ ``` ## Stopping point 2: workflow with full wildcards and sample dictionary ```python samples_dict = { 'SRR2584857_1': 'https://osf.io/aksmc/download', 'SRR2584403_1': 'https://osf.io/6jx7z/download', 'SRR2584404_1': 'https://osf.io/s24ky/download', 'SRR2584405_1': 'https://osf.io/7qek6/download', } print('rule all contains', expand("{s}.ecoli-rel606.vcf", s=list(samples_dict))) # default rule - first one in the file rule all: input: expand("{s}.ecoli-rel606.vcf", s=list(samples_dict)) ## all of the rules below here are now generic for ANY sample, ## as long as {sample}.fastq.gz exists, you can produce any of the ## rule outputs below. def lookup_url(wildcards): return samples_dict[wildcards.sample] rule download_data: output: "{sample}.fastq.gz" params: # use params instead of input for this, because # params do not need to be file names that exist ;) download_link = lookup_url shell: """ wget {params.download_link} -O {output} """ rule download_reference: output: "ecoli-rel606.fa.gz" shell: """ wget https://osf.io/8sm92/download -O {output} """ rule run_minimap: input: reads="{sample}.fastq.gz", genome="ecoli-rel606.fa.gz" output: "{sample}.ecoli-rel606.sam" conda: "env-variant-calling.yml" shell: """ minimap2 -ax sr {input.genome} {input.reads} > {output} """ rule sam_to_bam: input: rules.run_minimap.output output: "{sample}.ecoli-rel606.bam" conda: "env-variant-calling.yml" shell: """ samtools view -b -F 4 {input} > {output} """ rule sort_bam: output: "{sample}.ecoli-rel606.bam.sorted" input: "{sample}.ecoli-rel606.bam" conda: "env-variant-calling.yml" shell: """ samtools sort {input} > {output} """ rule ungz_genome: input: "ecoli-rel606.fa.gz", output: "ecoli-rel606.fa" shell: "gunzip -c {input} > {output}" rule call_variants: input: sorted_bam="{sample}.ecoli-rel606.bam.sorted", ungz_genome = "ecoli-rel606.fa", output: pileup = "{sample}.ecoli-rel606.pileup", bcf = "{sample}.ecoli-rel606.bcf", vcf = "{sample}.ecoli-rel606.vcf" conda: "env-variant-calling.yml" shell: """ bcftools mpileup -Ou -f {input.ungz_genome} {input.sorted_bam} > {output.pileup} bcftools call -mv -Ob {output.pileup} -o {output.bcf} bcftools view {output.bcf} > {output.vcf} """ ``` ## Stopping point 3: loading samples from a CSV file instead Put this at the top of the Snakefile: ```python import csv # loading the sample names and URLs from 'samples.csv' with open('samples.csv', 'rt') as fp: r = csv.DictReader(fp) samples_dict = {} for row in r: name = row['sample'] url = row['url'] samples_dict[name] = url ``` and write a `samples.csv` file like so: ```csvpreview= sample,url SRR2584857_1,https://osf.io/aksmc/download SRR2584403_1,https://osf.io/6jx7z/download SRR2584404_1,https://osf.io/s24ky/download SRR2584405_1,https://osf.io/7qek6/download ``` ## Other tutorials and links http://ivory.idyll.org/blog/2020-snakemake-hacks-collections-files.html https://training.nih-cfde.org/en/latest/Bioinformatics-Skills/Snakemake/

    Import from clipboard

    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 lost their connection.

    Create a note from template

    Create a note from template

    Oops...
    This template is not available.
    All
    • All
    • Team
    No template found.

    Create a template

    Delete template

    Do you really want to delete this template?

    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

    By clicking below, you agree to our terms of service.

    Sign in via Facebook Sign in via Twitter Sign in via GitHub Sign in via Dropbox Sign in via Google

    New to HackMD? Sign up

    Help

    Documents

    Tutorials
    YAML Metadata
    Slide Example
    Book Example

    Contacts

    Talk to us
    Report an issue
    Send us email

    Cheatsheet

    Example Syntax
    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~~
    19th 19^th^
    H2O H~2~O
    Inserted text ++Inserted text++
    Marked text ==Marked text==
    Link [link text](https:// "title")
    Image ![image alt](https:// "title")
    Code `Code`
    var i = 0;
    ```javascript
    var i = 0;
    ```
    :smile: :smile:
    Externals {%youtube youtube_id %}
    LaTeX $L^aT_eX$

    This is a alert area.

    :::info
    This is a alert area.
    :::

    Versions

    Versions and GitHub Sync

    Sign in to link this note to GitHub Learn more
    This note is not linked with GitHub Learn more
     
    Add badge Pull Push GitHub Link Settings

    Version named by    

    More Less
    • Edit
    • Delete

    Note content is identical to the latest version.
    Compare with
      Choose a version
      No search result
      Version not found

    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. Learn more

         Sign in to GitHub

        HackMD links with GitHub through a GitHub App. You can choose which repo to install our App.

        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

        Pull from GitHub

         
        File from GitHub
        File from HackMD

        GitHub Link Settings

        File linked

        Linked by
        File path
        Last synced branch

        Danger Zone

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

        Syncing

        Push failed

        Push successfully