HPCwiki - User contributions [en]

Whole genome alignment pipeline

2022-07-18T13:56:19Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/whole-genome-alignment here]

= Whole genome alignment pipeline =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/whole-genome-alignment


== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This pipeline aligns one or more genomes to a specified genome and plots the alignment.


==== Tools used: ====

* [https://github.com/lh3/minimap2 minimap2]
* R
** [https://github.com/tpoorten/dotPlotly/blob/master/pafCoordsDotPlotly.R pafCoordsDotPlotly.R]

{|
!align="center"| [[File:whole-genome-alignment-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files ===

<syntaxhighlight lang="yaml"># genome alignment parameters:
GENOME: /path/to/genome #genome fasta to be compared
COMPARISON_GENOME:
<genome1>: path/to/genome1.fasta
<genome2>: path/to/genome2.fasta
<genome3>: path/to/genome3.fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000

PREFIX: <prefix>

OUTDIR: /path/to/outdir</syntaxhighlight>
* GENOME: path to the genome fasta file (can be compressed). This is the genome that you want to be compared to all the others
* COMPARISON_GENOME: genome fasta (can be compressed) for whole genome comparison. Add your species name and the path to the fasta file. ex: chicken: /path/to/chicken.fna.gz. You can add several genomes, one on each line.
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there’s an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.
* PREFIX: name of your species (ex: turkey)
* OUTDIR: directory where snakemake will run and where the results will be written to

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>

== ADDITIONAL SET UP ==


=== Installing R packages ===

First load R: <code>module load R/4.0.2</code>

Enter the R environment by writing <code>R</code> and clicking enter. Install the packages:

<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>
If you get an error like this:

<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.


== RESULTS ==

The most important files and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''genome_alignment/{prefix}''vs''{species}.paf''' paf format file w

Whole genome alignment pipeline

2022-07-18T13:47:58Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/whole-genome-alignment here]

= Single Cell preprocessing pipeline =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/whole-genome-alignment


== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This pipeline aligns one or more genomes to a specified genome and plots the alignment.


==== Tools used: ====

* [https://github.com/lh3/minimap2 minimap2]
* R
** [https://github.com/tpoorten/dotPlotly/blob/master/pafCoordsDotPlotly.R pafCoordsDotPlotly.R]

{|
!align="center"| [[File:whole-genome-alignment-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files ===

<syntaxhighlight lang="yaml"># genome alignment parameters:
GENOME: /path/to/genome #genome fasta to be compared
COMPARISON_GENOME:
<genome1>: path/to/genome1.fasta
<genome2>: path/to/genome2.fasta
<genome3>: path/to/genome3.fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000

PREFIX: <prefix>

OUTDIR: /path/to/outdir</syntaxhighlight>
* GENOME: path to the genome fasta file (can be compressed). This is the genome that you want to be compared to all the others
* COMPARISON_GENOME: genome fasta (can be compressed) for whole genome comparison. Add your species name and the path to the fasta file. ex: chicken: /path/to/chicken.fna.gz. You can add several genomes, one on each line.
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there’s an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.
* PREFIX: name of your species (ex: turkey)
* OUTDIR: directory where snakemake will run and where the results will be written to

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>

== ADDITIONAL SET UP ==


=== Installing R packages ===

First load R: <code>module load R/4.0.2</code>

Enter the R environment by writing <code>R</code> and clicking enter. Install the packages:

<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>
If you get an error like this:

<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.


== RESULTS ==

The most important files and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''genome_alignment/{prefix}''vs''{species}.paf''' paf format file w

File:Whole-genome-alignment-workflow.png

2022-07-18T13:46:53Z

Moiti001:

Whole genome alignment pipeline

2022-07-18T13:46:02Z

Moiti001: added documentation

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/whole-genome-alignment here]

= Single Cell preprocessing pipeline =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/whole-genome-alignment


== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This pipeline aligns one or more genomes to a specified genome and plots the alignment.


==== Tools used: ====

* [https://github.com/lh3/minimap2 minimap2]
* R
** [https://github.com/tpoorten/dotPlotly/blob/master/pafCoordsDotPlotly.R pafCoordsDotPlotly.R]

{|
!align="center"| [[File:whole-genome-alignment-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files ===

<syntaxhighlight lang="yaml"># genome alignment parameters:
GENOME: /path/to/genome #genome fasta to be compared
COMPARISON_GENOME:
<genome1>: path/to/genome1.fasta
<genome2>: path/to/genome2.fasta
<genome3>: path/to/genome3.fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000

PREFIX: <prefix>

OUTDIR: /path/to/outdir</syntaxhighlight>
* GENOME: path to the genome fasta file (can be compressed). This is the genome that you want to be compared to all the others
* COMPARISON_GENOME: genome for whole genome comparison. Add your species name and the path to the fasta file. ex: chicken: /path/to/chicken.fna.gz. You can add several genomes, one on each line.
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there’s an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.
* PREFIX: name of your species (ex: turkey)
* OUTDIR: directory where snakemake will run and where the results will be written to

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>

== ADDITIONAL SET UP ==


=== Installing R packages ===

First load R: <code>module load R/4.0.2</code>

Enter the R environment by writing <code>R</code> and clicking enter. Install the packages:

<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>
If you get an error like this:

<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.


== RESULTS ==

The most important files and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''genome_alignment/{prefix}''vs''{species}.paf''' paf format file w

Bioinformatics tips tricks workflows

2022-07-18T13:41:09Z

Moiti001: added whole genome alignment pipeline page

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling| Nanopore assembly and variant calling pipeline]]
* [[Population variant calling pipeline | Population variant calling pipeline]]
* [[Single Cell preprocessing pipeline| Single Cell preprocessing pipeline]]
* [[Whole genome alignment pipeline | Whole genome alignment pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

File:Single-cell-processing-workflow.png

2022-07-18T09:03:43Z

Moiti001:

Single Cell preprocessing pipeline

2022-07-18T09:02:40Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/single-cell-data-processing here]

= Single Cell preprocessing pipeline =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/single-cell-data-processing


== First follow the instructions here ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This pipeline includes the first steps in the analysis of Single-cell data. 
The first step is getting the reference package for your species. This will be used for read alignment and gene expression quantification. If you're working with human or mouse, you can download the reference from the Cellranger website, if not, the pipeline can create the reference for you. (details below)

Once you have the reference package, the pipeline starts by running [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count Cellranger count]. Cellranger count performs alignment, filtering, barcode counting, and UMI counting. It uses the Chromium cellular barcodes to generate feature-barcode matrices, determine clusters, and perform gene expression analysis. 
If the fastq files are not named in the format accepted by Cellranger count: <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>, you can specify in the config file that these need to be renamed: option <code>RENAME: y</code> or you can rename them yourself to follow this naming convention.

The metrics from Cellranger count for all samples are combined into one file <code>cellranger_count_metrics_allsamples.tsv</code>. This will have information such as "estimated number of cells", and "mean reads per cell".

After the Cellranger count step, it's important to remove the ambient RNA, which is RNA that has been released from degraded or dying cells and is now in the cell suspension. The R package [https://github.com/constantAmateur/SoupX SoupX] is used to correct for ambient RNA. In addition to the output files with the corrected data, one html document is created per sample processed (<code>2_ambient_RNA_correction/Ambient_RNA_correction_<sample>.html</code>). This html file shows the code used to perform the ambient RNA correction, as well as a few plots that illustrate this process - for the 5 most affected genes and for 5 random genes:

* Plot 1: in which cells the gene is expressed
* Plot 2: ratio of observed to expected counts
* Plot 3: change in expression due to correction

Once the data has been corrected for ambient RNA, it's time for quality control filtering. This is a step that depends on the cell type, library preparation method used, etc, so you should always check if the default parameters make sense, use your own, or even run several times with different ones.

QC is run for every sample separately. First [https://scanpy.readthedocs.io/en/stable/ Scanpy] calculates some general QC metrics for genes and cells. It will also calculate the proportion of counts for mitochondrial genes. Several plots will be created to help assess the quality of the data: 
Before filtering:

* Violin plots showing:
** n_genes_by_counts: number of genes with positive counts in a cell
** total_counts: total number of counts for a cell
** pct_counts_mt: proportion of mitochondrial counts for a cell
* Scatter plot showing :
** total_counts vs pct_counts_mt
** total counts vs n_genes_by_counts

After filtering:

* Percentage of counts per gene for the top 20 genes after filtering
* Violin plots showing:
** n_genes_by_counts: number of genes with positive counts in a cell
** total_counts: total number of counts for a cell
** pct_counts_mt: proportion of mitochondrial counts for a cell

The final preprocessing step is doublet removal with [https://github.com/swolock/scrublet Scrublet]. This step may be run more than once to determine the ideal doublet score threshold. The histogram shown in <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code> should show a biomodal distribution, and the threshold shown in the "simulated doublets" plot should be at the minium between the two modes. The first run should be with parameter <code>SCRUB_TRESHOLD:</code>. Once the first run has finished, you should look at the histograms of all samples and see if you need to change the treshold. If you do, set it in the config file as explained further below and run the <code>remove_doublets</code> step again.

There will be an <code>h5ad</code> object containing the preprocessed data for each sample - ambient RNA removed, QC filtered and doublets removed - in the <code>4_Doublets</code> directory. 
This file can be used for further analysis.


==== Tools used ====

* Cellranger:
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkgtf mkgtf] - filter GTF. [default: off]
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkref mkref] - create reference
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count count] - create feature counts
* R
** combine Cellranger count sample metrics
** [https://github.com/constantAmateur/SoupX SoupX] - remove ambient RNA
* Python
** [https://scanpy.readthedocs.io/en/stable/index.html Scanpy] - QC filtering
** [https://github.com/swolock/scrublet Scrublet] - Doublet removal

{|
!align="center"| [[File:single-cell-processing-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files and set parameters ===

<syntaxhighlight lang="yaml">DATA: /path/to/data/dir
OUTDIR: /path/to/outdir

# mkref options
MKREF: <y/n>
FASTA: /path/to/fasta
GTF: /path/to/gtf # if creating own reference
REF_VERSION:
- "--ref-version=<version>"
CR_MKREF_EXTRA: ""

# Filter GTF
FILTER_GTF: y
# # see here for available biotypes https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkgtf
ATTRIBUTES:
- "--attribute=<biotype>"

PREFIX: <species>

# rename fastq files
RENAME: <y/n>

# Cell ranger count options
# https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count#cr-count
CR_COUNT_extra: ""

# QC parameters
MITO_PERCENTAGE: 10 # keep cells with less than X% mitochondrial read fraction
NUMBER_GENES_PER_CELL: 500 # keep cells with more than X genes
NUMBER_UMI_PER_CELL: 1000 # keep cells with more than X UMIs
ENSEMBLE_BIOMART_SPECIES: "<species>" # ensembl biomart species used to get the mitochondrial genes for that species

# threshold doublet score (should be at the minimum between two modes of the simulated doublet histogram)
SCRUB_THRESHOLD:
<sample1>: <value>
<sample2>: <empty></syntaxhighlight>
* '''DATA''' - path to directory containing fastq files. Preferrably, the files should be named in the format accepted by Cellranger Count <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>. If they are, set <code>RENAME: n</code>. If not, they should be in the format <code><sample>_R1.fastq.gz</code>. In this case, you should set <code>RENAME: y</code> so that the pipeline will rename the files according to the necessary format for Cellranger Count. The path can be to a directory that contains a subdirectory per sample. For example:
<pre>```text
DATA
├── SAMPLE_1
│ ├── sample_1_R1.fastq.gz
│ └── sample_1_R2.fastq.gz
└── SAMPLE_2
├── sample_2_R1.fastq.gz
└── sample_2_R2.fastq.gz
```

Or to a directory with fastqs for all samples:

```text
DATA
├── sample_1_R1.fastq.gz
├── sample_1_R2.fastq.gz
├── sample_2_R1.fastq.gz
└── sample_2_R2.fastq.gz
```</pre>
* '''OUTDIR''' - directory where snakemake will run and where the results will be written to. 
If you don't want the results to be written to a new directory, open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* '''MKREF''' <code>y</code> if Cell Ranger doesn't provide a reference package for your species (currently, the species available are Human and Mouse). <code>n</code> if you're using an existing reference package. In this case, you should create a directory in the pipeline directory named <code><prefix>_genome</code>, and download the reference from [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest here].

-----

'''You only need to set the following if <code>MKREF: y</code>:'''

* '''FASTA''' - path to fasta file
* '''GTF''': path to gft file
* '''PREFIX''': name of your species. Used to name the directory containing the reference package
* '''REF_VERSION''': Reference version string to include with reference
* '''CR_MKREF_EXTRA''': any other options for cellranger mkref. [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#singl see here]

-----

<ul>
<li>'''PREFIX''' - The name of your organism. The reference package used for cellranger count will be in the <code><prefix>_genome</code> directory</li>
<li>'''RENAME''' - <code>y</code> if your input fastqs are not named in this format <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>. Use <code>n</code> if they are.</li>
<li>Options for Cellrange count:
<ul>
<li>'''CR_COUNT_extra''' - any other options for cellranger count. [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count#cr-count Find other options here]. [Default: ""]</li></ul>
</li>
<li>QC parameters
<ul>
<li>'''MITO_PERCENTAGE''' - Keep cells with less than X% mitochondrial read fraction. [Default: 10]</li>
<li>'''NUMBER_GENES_PER_CELL''' - keep cells with more than X genes. [Default: 500]</li>
<li>'''NUMBER_UMI_PER_CELL''' - keep cells with more than X UMIs. [Default: 1000]</li>
<li>'''ENSEMBLE_BIOMART_SPECIES''' - ensembl biomart species used to get the mitochondrial genes</li></ul>
</li>
<li>Doublet removal score
<ul>
<li>'''SCRUB_THRESHOLD''' - threshold doublet score. It should be at the minimum between two modes of the simulated doublet histogram. In the first run it should be run as <code>SCRUB_TRESHOLD:</code> (with no parameters). After that is done, for each sample you should then look at the <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code> plot and see if the vertical line on the "simulated doublets" plot is at the minimum between the two modes. If it's not, you should manually set it in the config file as:
<pre>```yaml
SCRUB_THRESHOLD:
<sample 1>: <value>
<sample 2>: <empty>
```

In `SCRUB_THRESHOLD` there should be a line for each sample, even if you don't need to set the threshold for that sample. If you need to change the treshold, set `<sample>: <value>`, if not, set `<sample>:` (without value).</pre>

== Additional set up ==


=== Install Cellranger ===

Follow the instructions [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation here]

# First download the package from the [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest downloads page]. For example:

<syntaxhighlight lang="sh">wget -O cellranger-6.1.2.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-6.1.2.tar.gz?Expires=1648081234&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci02LjEuMi50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2NDgwODEyMzR9fX1dfQ__&Signature=HqCwx6eBEj~Lyw7C7UvsMAHzUH9aiPSM5yFcyflZiL2JRIwqzY2VWz1COtDQHNoJ48Ve41LZ5Q3eGv1yaAEf88SGhtxRUb2wJhFvvixBoR550bQ2wK7qfL6buLL9~u7MPw4q0-c1adXaSCm6otd6Xn0x2FIpZimOGJMYI9QEvNStN1Hi6MH4ZUOHGFFRBAvxlRxHmYBk-Vr~6qdc7nFXJW0C8OBWTn2g~XSKZRD50B5G5StMis0lLmgXZbRS0htQu8LPuUp8ZxqxQv20m9-HV9jEDVYEUP1sNJzAHGhAtq1FajN572Lptq0cWES8fheMexht1l-wRbQA-yOKAp7Bzg__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"</syntaxhighlight>
Attention: the <code>cellranger-6.1.2.tar.gz</code> name is an example, when you use this pipeline there might be a more recent version available.

<ol start="2" style="list-style-type: decimal;">
<li>Unpack Cellranger:</li></ol>

<syntaxhighlight lang="sh">tar -xzvf cellranger-6.1.2.tar.gz</syntaxhighlight>
This will create a new directory, <code>cellranger-6.1.2</code>, that contains cellranger and its dependencies.

<ol start="3" style="list-style-type: decimal;">
<li>Place the path to the <code>cellranger-6.1.2</code> (or the version you installed) in the config.yaml file. It should look like this</li></ol>

<syntaxhighlight lang="yaml">CELLRANGER_PATH: /path/to/cellranger-6.1.2</syntaxhighlight>
<blockquote>don't add the backslash ("\") after the directory name
</blockquote>

=== Reference package ===

If you're working with human or mouse data, download the reference from [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest here] and place it in a folder in the pipeline directory called <code><prefix>_genome</code>. 
If you're working with another organism, download the fasta file and gtf file for your organism and place them in a directory called <code><prefix>_genome</code> directory (should be in the pipeline main directory - where the Snakefile is). You can download these from [https://www.ensembl.org/index.html Ensembl].


=== How to run ===

<ol style="list-style-type: decimal;">
<li>Dry run to check if everything was correctly set up and if the pipeline is ready to run
<syntaxhighlight lang="sh">snakemake -np</syntaxhighlight></li>
<li>If all looks good, run the pipeline with
<syntaxhighlight lang="sh">snakemake --profile <name of hpc profile></syntaxhighlight></li>
<li>Once you have the <code>remove_doublets</code> results for each sample, you should look at the histogram, either in the jupyter notebook <code>4_Doublets/processed_notebook_<sample>.ipynb</code> or in the saved plot <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code>. The vertical bar (threshold) in the "simulated doublets" plot should be at the lowest point between the two modes. If not, you'll need to set it in the config file. Only change the threshold for the samples that need it. For some, the automatically set threshold may be good. To do this, edit the config file. Add your sample name (the same as the directory names in <code>4_Doublets</code>) and the new threshold or, in case you don't need to add a threshold, add the sample name and nothing in front of it.
<syntaxhighlight lang="yaml">SCRUB_THRESHOLD:
<sample 1>: <value>
<sample 2>:
<sample 3>: <value></syntaxhighlight>
In this case, the scrublet step for <code>sample 1</code> and <code>sample 3</code> will run with the user defined threshold and <code>sample 2</code> will run with the automatically defined treshold.</li>
<li>After changing the thresholds, you'll need to run the <code>remove_doublets</code> step again. Before you do this, you need to copy the previous results (<code>4_Doublets</code>) to another directory, or you need to delete those results. Example to move the results to another directory:
<syntaxhighlight lang="sh">mkdir -p 4_Doublets_first_run
mv -r 4_Doublets/* 4_Doublets_first_run</syntaxhighlight></li></ol>
<ol start="5" style="list-style-type: decimal;">
<li>Once that's done you can rerun the <code>remove_doublets</code> step with
<syntaxhighlight lang="sh">snakemake -np --forcerun remove_doublets
snakemake --profile <profile name> --forcerun remove_doublets</syntaxhighlight></li></ol>


== RESULTS ==

You will have results for each step of the pipeline.

<ul>
<li>1_renamed: directory with renamed fastq files softlinked</li>
<li>2_ambient_RNA_correction: directory containing results from ambient RNA correction
<ul>
<li>Ambient_RNA_correction_<sample>.html: shows the code used for the ambient RNA correction, as well as a few plots that illustrate this process - for the 5 most affected genes and for 5 random genes:
Plot 1: in which cells the gene is expressed 
Plot 2: ratio of observed to expected counts 
Plot 3: change in expression due to correction</li></ul>
</li>
<li>2_ambient_RNA_correction_data: inside, the ambient RNA corrected data for each sample is its corresponding directory.</li>
<li>3_QC: directory containing results from QC:
<ul>
<li>processed_notebook_<sample>.ipynp - Jupyter notebooks used to calculate QC for each sample. These are interactive and can be used to do further QC.</li>
<li><sample>.h5ad - Filtered data.</li>
<li><sample> - Directory with QC plots

===== Description of the QC plots =====

Before filtering:
<ul>
<li>Violin plots showing:
<ul>
<li>n_genes_by_counts: number of genes with positive counts in a cell</li>
<li>total_counts: total number of counts for a cell</li>
<li>pct_counts_mt: proportion of mitochondrial counts for a cell</li></ul>
</li>
<li>Scatter plot showing :
<ul>
<li>total_counts vs pct_counts_mt</li>
<li>total counts vs n_genes_by_counts</li></ul>
</li></ul>
After filtering:

<pre> - Percentage of counts per gene for the top 20 genes after filtering
- Violin plots showing:
- n_genes_by_counts: number of genes with positive counts in a cell
- total_counts: total number of counts for a cell
- pct_counts_mt: proportion of mitochondrial counts for a cell
These jupyter notebooks are interactive and can be used to do further QC control.</pre>
* 4_Doublets: directory containing results from doublet removal
** processed_notebook_<sample>.ipynb - Jupyter notebooks used to do the doublet removal step for each sample. These are interactive and can be used to test different scrublet tresholds
** <sample> - directory with saved plots
** <sample>_doublets.h5ad - filtered data
* <sample> - directory containing output from Cellranger count. See [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/overview here] for more information
** outs
*** web_summary.html - summary metrics and automated secondary analysis results. If an issue was detected during the pipeline run, an alert appears on this page.
*** metrics_summary.csv - metrics like "estimated number of cells"
*** possorted_genome_bam.bam - BAM file containing position-sorted reads aligned to the genome and transcriptome, as well as unaligned reads. Each read in this BAM file has Chromium cellular and molecular barcode information attached.
*** raw_feature_bc_matrix.h5 - Contains every barcode from the fixed list of known-good barcode sequences that has at least one read. This includes background and cell associated barcodes
*** filtered_feature_bc_matrix.h5 - Contains only detected cell-associated barcodes. For Targeted Gene Expression samples, non-targeted genes are removed from the filtered matrix.
*** analysis - directory containing secondary analysis results: clustering, differential expression analysis, PCA, t-SNE, UMAP
*** molecule_info.h5 - contains per-molecule information for all molecules that contain a valid barcode, a valid UMI, and were assigned with high confidence to a gene or Feature Barcode.
*** cloupe.cloup - file to be used with [https://support.10xgenomics.com/single-cell-gene-expression/software/visualization/latest/what-is-loupe-cell-browser Loupe Browser]

== Common issues ==


=== Only see <code>create_file_log</code> and <code>combine_cellranger_counter_metrics</code> when running <code>snakemake -np</code> ===

This can be because your fastq files are not named correctly and the <code>rename</code> option is not set to <code>y</code>.

Single Cell preprocessing pipeline

2022-07-18T08:54:06Z

Moiti001: added documentation

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/single-cell-data-processing here]

= Single Cell preprocessing pipeline =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/single-cell-data-processing


== First follow the instructions here ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This pipeline includes the first steps in the analysis of Single-cell data. 
The first step is getting the reference package for your species. This will be used for read alignment and gene expression quantification. If you're working with human or mouse, you can download the reference from the Cellranger website, if not, the pipeline can create the reference for you. (details below)

Once you have the reference package, the pipeline starts by running [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count Cellranger count]. Cellranger count performs alignment, filtering, barcode counting, and UMI counting. It uses the Chromium cellular barcodes to generate feature-barcode matrices, determine clusters, and perform gene expression analysis. 
If the fastq files are not named in the format accepted by Cellranger count: <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>, you can specify in the config file that these need to be renamed: option <code>RENAME: y</code> or you can rename them yourself to follow this naming convention.

The metrics from Cellranger count for all samples are combined into one file <code>cellranger_count_metrics_allsamples.tsv</code>. This will have information such as "estimated number of cells", and "mean reads per cell".

After the Cellranger count step, it's important to remove the ambient RNA, which is RNA that has been released from degraded or dying cells and is now in the cell suspension. The R package [https://github.com/constantAmateur/SoupX SoupX] is used to correct for ambient RNA. In addition to the output files with the corrected data, one html document is created per sample processed (<code>2_ambient_RNA_correction/Ambient_RNA_correction_<sample>.html</code>). This html file shows the code used to perform the ambient RNA correction, as well as a few plots that illustrate this process - for the 5 most affected genes and for 5 random genes:

* Plot 1: in which cells the gene is expressed
* Plot 2: ratio of observed to expected counts
* Plot 3: change in expression due to correction

Once the data has been corrected for ambient RNA, it's time for quality control filtering. This is a step that depends on the cell type, library preparation method used, etc, so you should always check if the default parameters make sense, use your own, or even run several times with different ones.

QC is run for every sample separately. First [https://scanpy.readthedocs.io/en/stable/ Scanpy] calculates some general QC metrics for genes and cells. It will also calculate the proportion of counts for mitochondrial genes. Several plots will be created to help assess the quality of the data: 
Before filtering:

* Violin plots showing:
** n_genes_by_counts: number of genes with positive counts in a cell
** total_counts: total number of counts for a cell
** pct_counts_mt: proportion of mitochondrial counts for a cell
* Scatter plot showing :
** total_counts vs pct_counts_mt
** total counts vs n_genes_by_counts

After filtering:

* Percentage of counts per gene for the top 20 genes after filtering
* Violin plots showing:
** n_genes_by_counts: number of genes with positive counts in a cell
** total_counts: total number of counts for a cell
** pct_counts_mt: proportion of mitochondrial counts for a cell

The final preprocessing step is doublet removal with [https://github.com/swolock/scrublet Scrublet]. This step may be run more than once to determine the ideal doublet score threshold. The histogram shown in <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code> should show a biomodal distribution, and the threshold shown in the "simulated doublets" plot should be at the minium between the two modes. The first run should be with parameter <code>SCRUB_TRESHOLD:</code>. Once the first run has finished, you should look at the histograms of all samples and see if you need to change the treshold. If you do, set it in the config file as explained further below and run the <code>remove_doublets</code> step again.

There will be an <code>h5ad</code> object containing the preprocessed data for each sample - ambient RNA removed, QC filtered and doublets removed - in the <code>4_Doublets</code> directory. 
This file can be used for further analysis.


==== Tools used ====

* Cellranger:
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkgtf mkgtf] - filter GTF. [default: off]
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkref mkref] - create reference
** [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count count] - create feature counts
* R
** combine Cellranger count sample metrics
** [https://github.com/constantAmateur/SoupX SoupX] - remove ambient RNA
* Python
** [https://scanpy.readthedocs.io/en/stable/index.html Scanpy] - QC filtering
** [https://github.com/swolock/scrublet Scrublet] - Doublet removal

{|
!align="center"| [[File:https://github.com/CarolinaPB/single-cell-data-processing/blob/master/workflow.png|DAG]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files and set parameters ===

<syntaxhighlight lang="yaml">DATA: /path/to/data/dir
OUTDIR: /path/to/outdir

# mkref options
MKREF: <y/n>
FASTA: /path/to/fasta
GTF: /path/to/gtf # if creating own reference
REF_VERSION:
- "--ref-version=<version>"
CR_MKREF_EXTRA: ""

# Filter GTF
FILTER_GTF: y
# # see here for available biotypes https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkgtf
ATTRIBUTES:
- "--attribute=<biotype>"

PREFIX: <species>

# rename fastq files
RENAME: <y/n>

# Cell ranger count options
# https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count#cr-count
CR_COUNT_extra: ""

# QC parameters
MITO_PERCENTAGE: 10 # keep cells with less than X% mitochondrial read fraction
NUMBER_GENES_PER_CELL: 500 # keep cells with more than X genes
NUMBER_UMI_PER_CELL: 1000 # keep cells with more than X UMIs
ENSEMBLE_BIOMART_SPECIES: "<species>" # ensembl biomart species used to get the mitochondrial genes for that species

# threshold doublet score (should be at the minimum between two modes of the simulated doublet histogram)
SCRUB_THRESHOLD:
<sample1>: <value>
<sample2>: <empty></syntaxhighlight>
* '''DATA''' - path to directory containing fastq files. Preferrably, the files should be named in the format accepted by Cellranger Count <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>. If they are, set <code>RENAME: n</code>. If not, they should be in the format <code><sample>_R1.fastq.gz</code>. In this case, you should set <code>RENAME: y</code> so that the pipeline will rename the files according to the necessary format for Cellranger Count. The path can be to a directory that contains a subdirectory per sample. For example:
<pre>```text
DATA
├── SAMPLE_1
│ ├── sample_1_R1.fastq.gz
│ └── sample_1_R2.fastq.gz
└── SAMPLE_2
├── sample_2_R1.fastq.gz
└── sample_2_R2.fastq.gz
```

Or to a directory with fastqs for all samples:

```text
DATA
├── sample_1_R1.fastq.gz
├── sample_1_R2.fastq.gz
├── sample_2_R1.fastq.gz
└── sample_2_R2.fastq.gz
```</pre>
* '''OUTDIR''' - directory where snakemake will run and where the results will be written to. 
If you don't want the results to be written to a new directory, open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* '''MKREF''' <code>y</code> if Cell Ranger doesn't provide a reference package for your species (currently, the species available are Human and Mouse). <code>n</code> if you're using an existing reference package. In this case, you should create a directory in the pipeline directory named <code><prefix>_genome</code>, and download the reference from [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest here].

-----

'''You only need to set the following if <code>MKREF: y</code>:'''

* '''FASTA''' - path to fasta file
* '''GTF''': path to gft file
* '''PREFIX''': name of your species. Used to name the directory containing the reference package
* '''REF_VERSION''': Reference version string to include with reference
* '''CR_MKREF_EXTRA''': any other options for cellranger mkref. [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#singl see here]

-----

<ul>
<li>'''PREFIX''' - The name of your organism. The reference package used for cellranger count will be in the <code><prefix>_genome</code> directory</li>
<li>'''RENAME''' - <code>y</code> if your input fastqs are not named in this format <code>[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz</code>. Use <code>n</code> if they are.</li>
<li>Options for Cellrange count:
<ul>
<li>'''CR_COUNT_extra''' - any other options for cellranger count. [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/count#cr-count Find other options here]. [Default: ""]</li></ul>
</li>
<li>QC parameters
<ul>
<li>'''MITO_PERCENTAGE''' - Keep cells with less than X% mitochondrial read fraction. [Default: 10]</li>
<li>'''NUMBER_GENES_PER_CELL''' - keep cells with more than X genes. [Default: 500]</li>
<li>'''NUMBER_UMI_PER_CELL''' - keep cells with more than X UMIs. [Default: 1000]</li>
<li>'''ENSEMBLE_BIOMART_SPECIES''' - ensembl biomart species used to get the mitochondrial genes</li></ul>
</li>
<li>Doublet removal score
<ul>
<li>'''SCRUB_THRESHOLD''' - threshold doublet score. It should be at the minimum between two modes of the simulated doublet histogram. In the first run it should be run as <code>SCRUB_TRESHOLD:</code> (with no parameters). After that is done, for each sample you should then look at the <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code> plot and see if the vertical line on the "simulated doublets" plot is at the minimum between the two modes. If it's not, you should manually set it in the config file as:
<pre>```yaml
SCRUB_THRESHOLD:
<sample 1>: <value>
<sample 2>: <empty>
```

In `SCRUB_THRESHOLD` there should be a line for each sample, even if you don't need to set the threshold for that sample. If you need to change the treshold, set `<sample>: <value>`, if not, set `<sample>:` (without value).</pre>

== Additional set up ==


=== Install Cellranger ===

Follow the instructions [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation here]

# First download the package from the [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest downloads page]. For example:

<syntaxhighlight lang="sh">wget -O cellranger-6.1.2.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-6.1.2.tar.gz?Expires=1648081234&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci02LjEuMi50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2NDgwODEyMzR9fX1dfQ__&Signature=HqCwx6eBEj~Lyw7C7UvsMAHzUH9aiPSM5yFcyflZiL2JRIwqzY2VWz1COtDQHNoJ48Ve41LZ5Q3eGv1yaAEf88SGhtxRUb2wJhFvvixBoR550bQ2wK7qfL6buLL9~u7MPw4q0-c1adXaSCm6otd6Xn0x2FIpZimOGJMYI9QEvNStN1Hi6MH4ZUOHGFFRBAvxlRxHmYBk-Vr~6qdc7nFXJW0C8OBWTn2g~XSKZRD50B5G5StMis0lLmgXZbRS0htQu8LPuUp8ZxqxQv20m9-HV9jEDVYEUP1sNJzAHGhAtq1FajN572Lptq0cWES8fheMexht1l-wRbQA-yOKAp7Bzg__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"</syntaxhighlight>
Attention: the <code>cellranger-6.1.2.tar.gz</code> name is an example, when you use this pipeline there might be a more recent version available.

<ol start="2" style="list-style-type: decimal;">
<li>Unpack Cellranger:</li></ol>

<syntaxhighlight lang="sh">tar -xzvf cellranger-6.1.2.tar.gz</syntaxhighlight>
This will create a new directory, <code>cellranger-6.1.2</code>, that contains cellranger and its dependencies.

<ol start="3" style="list-style-type: decimal;">
<li>Place the path to the <code>cellranger-6.1.2</code> (or the version you installed) in the config.yaml file. It should look like this</li></ol>

<syntaxhighlight lang="yaml">CELLRANGER_PATH: /path/to/cellranger-6.1.2</syntaxhighlight>
<blockquote>don't add the backslash ("\") after the directory name
</blockquote>

=== Reference package ===

If you're working with human or mouse data, download the reference from [https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest here] and place it in a folder in the pipeline directory called <code><prefix>_genome</code>. 
If you're working with another organism, download the fasta file and gtf file for your organism and place them in a directory called <code><prefix>_genome</code> directory (should be in the pipeline main directory - where the Snakefile is). You can download these from [https://www.ensembl.org/index.html Ensembl].


=== How to run ===

<ol style="list-style-type: decimal;">
<li>Dry run to check if everything was correctly set up and if the pipeline is ready to run
<syntaxhighlight lang="sh">snakemake -np</syntaxhighlight></li>
<li>If all looks good, run the pipeline with
<syntaxhighlight lang="sh">snakemake --profile <name of hpc profile></syntaxhighlight></li>
<li>Once you have the <code>remove_doublets</code> results for each sample, you should look at the histogram, either in the jupyter notebook <code>4_Doublets/processed_notebook_<sample>.ipynb</code> or in the saved plot <code>4_Doublets/<sample>/histogram_<sample>_doublets.pdf</code>. The vertical bar (threshold) in the "simulated doublets" plot should be at the lowest point between the two modes. If not, you'll need to set it in the config file. Only change the threshold for the samples that need it. For some, the automatically set threshold may be good. To do this, edit the config file. Add your sample name (the same as the directory names in <code>4_Doublets</code>) and the new threshold or, in case you don't need to add a threshold, add the sample name and nothing in front of it.
<syntaxhighlight lang="yaml">SCRUB_THRESHOLD:
<sample 1>: <value>
<sample 2>:
<sample 3>: <value></syntaxhighlight>
In this case, the scrublet step for <code>sample 1</code> and <code>sample 3</code> will run with the user defined threshold and <code>sample 2</code> will run with the automatically defined treshold.</li>
<li>After changing the thresholds, you'll need to run the <code>remove_doublets</code> step again. Before you do this, you need to copy the previous results (<code>4_Doublets</code>) to another directory, or you need to delete those results. Example to move the results to another directory:
<syntaxhighlight lang="sh">mkdir -p 4_Doublets_first_run
mv -r 4_Doublets/* 4_Doublets_first_run</syntaxhighlight></li></ol>
<ol start="5" style="list-style-type: decimal;">
<li>Once that's done you can rerun the <code>remove_doublets</code> step with
<syntaxhighlight lang="sh">snakemake -np --forcerun remove_doublets
snakemake --profile <profile name> --forcerun remove_doublets</syntaxhighlight></li></ol>


== RESULTS ==

You will have results for each step of the pipeline.

<ul>
<li>1_renamed: directory with renamed fastq files softlinked</li>
<li>2_ambient_RNA_correction: directory containing results from ambient RNA correction
<ul>
<li>Ambient_RNA_correction_<sample>.html: shows the code used for the ambient RNA correction, as well as a few plots that illustrate this process - for the 5 most affected genes and for 5 random genes:
Plot 1: in which cells the gene is expressed 
Plot 2: ratio of observed to expected counts 
Plot 3: change in expression due to correction</li></ul>
</li>
<li>2_ambient_RNA_correction_data: inside, the ambient RNA corrected data for each sample is its corresponding directory.</li>
<li>3_QC: directory containing results from QC:
<ul>
<li>processed_notebook_<sample>.ipynp - Jupyter notebooks used to calculate QC for each sample. These are interactive and can be used to do further QC.</li>
<li><sample>.h5ad - Filtered data.</li>
<li><sample> - Directory with QC plots

===== Description of the QC plots =====

Before filtering:
<ul>
<li>Violin plots showing:
<ul>
<li>n_genes_by_counts: number of genes with positive counts in a cell</li>
<li>total_counts: total number of counts for a cell</li>
<li>pct_counts_mt: proportion of mitochondrial counts for a cell</li></ul>
</li>
<li>Scatter plot showing :
<ul>
<li>total_counts vs pct_counts_mt</li>
<li>total counts vs n_genes_by_counts</li></ul>
</li></ul>
After filtering:

<pre> - Percentage of counts per gene for the top 20 genes after filtering
- Violin plots showing:
- n_genes_by_counts: number of genes with positive counts in a cell
- total_counts: total number of counts for a cell
- pct_counts_mt: proportion of mitochondrial counts for a cell
These jupyter notebooks are interactive and can be used to do further QC control.</pre>
* 4_Doublets: directory containing results from doublet removal
** processed_notebook_<sample>.ipynb - Jupyter notebooks used to do the doublet removal step for each sample. These are interactive and can be used to test different scrublet tresholds
** <sample> - directory with saved plots
** <sample>_doublets.h5ad - filtered data
* <sample> - directory containing output from Cellranger count. See [https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/overview here] for more information
** outs
*** web_summary.html - summary metrics and automated secondary analysis results. If an issue was detected during the pipeline run, an alert appears on this page.
*** metrics_summary.csv - metrics like "estimated number of cells"
*** possorted_genome_bam.bam - BAM file containing position-sorted reads aligned to the genome and transcriptome, as well as unaligned reads. Each read in this BAM file has Chromium cellular and molecular barcode information attached.
*** raw_feature_bc_matrix.h5 - Contains every barcode from the fixed list of known-good barcode sequences that has at least one read. This includes background and cell associated barcodes
*** filtered_feature_bc_matrix.h5 - Contains only detected cell-associated barcodes. For Targeted Gene Expression samples, non-targeted genes are removed from the filtered matrix.
*** analysis - directory containing secondary analysis results: clustering, differential expression analysis, PCA, t-SNE, UMAP
*** molecule_info.h5 - contains per-molecule information for all molecules that contain a valid barcode, a valid UMI, and were assigned with high confidence to a gene or Feature Barcode.
*** cloupe.cloup - file to be used with [https://support.10xgenomics.com/single-cell-gene-expression/software/visualization/latest/what-is-loupe-cell-browser Loupe Browser]

== Common issues ==


=== Only see <code>create_file_log</code> and <code>combine_cellranger_counter_metrics</code> when running <code>snakemake -np</code> ===

This can be because your fastq files are not named correctly and the <code>rename</code> option is not set to <code>y</code>.

Nanopore assembly and variant calling

2022-06-13T14:36:22Z

Moiti001: added link to github repo

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/nanopore-assembly here]
= Assemble nanopore reads and do variant calling with short and long reads =
Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/nanopore-assembly

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline that uses <code>Flye</code> to create a nanopore assembly. It also does variant calling with long and short reads. 
The pipeline starts by using <code>porechop</code> to trim the adaptors, then it uses <code>Flye</code> to create the assembly. After that, <code>ntLink-arks</code> from <code>Lonstitch</code> is used to scaffold the assembly using the nanopore reads. The scaffolded assembly is polished with <code>polca</code>. <code>Bwa-mem2</code> is used to map the short reads to the assembly and <code>Freebayes</code> to do variant calling using these reads. <code>Minimap2</code> is used to map the long reads to the assembly, and <code>longshot</code> for variant calling using these. In the end, in addition to your assembly and variant calling results, you'll also get assembly statistics and busco scores before and after the polishing.

==== Tools used: ====

* [https://github.com/rrwick/Porechop Porechop] - trim adaptors
* [https://github.com/fenderglass/Flye Flye] - assembly
* [https://github.com/lh3/seqtk Seqtk] - convert fasta to one line fasta
* [https://github.com/bcgsc/longstitch LongStitch (ntLink-arks)] - scaffolding with nanopore reads
* [https://busco.ezlab.org/ BUSCO] - assess assembly completeness
* [https://github.com/alekseyzimin/masurca MaSuRCA (polca)] - polish assembly
* Python - get assembly stats
* [https://github.com/lh3/minimap2 Minimap2] - map long reads to reference. Genome alignment
* [http://www.htslib.org/ Samtools] - sort and index mapped reads and vcf files
* [https://github.com/pjedge/longshot Longshot] - variant calling with nanopore reads
* [https://github.com/bwa-mem2/bwa-mem2 Bwa-mem2] - map short reads to reference
* [https://github.com/freebayes/freebayes Freebayes] - variant calling using short reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* R - [https://github.com/tpoorten/dotPlotly pafCoordsDotPlotly] - plot genome alignment

{|
!align="center"| [[File:nanopore-assembly-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>LONGREADS: <nanopore_reads.fq.gz>
SHORTREADS:
- /path/to/short/reads_1.fq.gz
- /path/to/short/reads_2.fq.gz
GENOME_SIZE: <approximate genome size>
PREFIX: <prefix>
OUTDIR: /path/to/outdir
BUSCO_LINEAGE:
- <lineage>

# genome alignment parameters:
COMPARISON_GENOME:
<species>: /path/to/genome/fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000</pre>
* LONGREADS - name of file with long reads. This file should be in the working directory (where this config and the Snakefile are)
* SHORTREADS - paths to short reads fq.gz
* GENOME_SIZE - approximate genome size <code>haploid genome size (bp)(e.g. '3e9' for human genome)</code> from [https://github.com/bcgsc/longstitch#full-help-page longstitch]
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* BUSCO_LINEAGE - lineage used for busco. Can be one or more (one per line). To see available lineages run <code>busco --list-datasets</code>
* COMPARISON_GENOME - genome for whole genome comparison. Add your species name and the path to the fasta file. ex: <code>chicken: /path/to/chicken.fna.gz</code>. You can add several genomes, one on each line.
** If you don't want to run the genome alignment step, comment out

<pre>COMPARISON_GENOME:
<species>: /path/to/genome/fasta</pre>
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there's an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.

-----
If you have your long reads in several fastq files and need to create one file compressed file with all the reads:

# In your pipeline directory create one file with all the reads

<pre>cat /path/to/fastq/directory/*.fastq > <name of file>.fq</pre>
<ol start="2" style="list-style-type: decimal;">
<li>Compress the file you just created:</li></ol>

<pre>gzip <name of file>.fq</pre>

=== After installing the conda environment (first step of this guide) you'll need to edit the polca.sh file. ===

First go to the directory where miniconda3 is installed (usually your home directory). Go to <code>/<home>/miniconda/envs/<env_name>/bin</code> and open the file <code>polca.sh</code>. In my case the path looks like this: <code>/home/WUR/<username>/miniconda3/envs/<env_name>/bin/</code>. In your editor open <code>polca.sh</code> and replace this line:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) $BASM.alignSorted 2>>samtools.err && \</pre>
With this:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) -o $BASM.alignSorted.bam 2>>samtools.err && \</pre>
== RESULTS ==

The most important files are and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''{prefix}_oneline.k32.w100.ntLink-arks.longstitch-scaffolds.fa.PolcaCorrected.fa''' final assembly
** assembly_stats_<prefix>.txt file with assembly statistics for the final assembly
** '''variant_calling''' directory with variant calling VCF files with long and short reads, as well as VCF stats
*** {prefix}_shortreads.vcf.gz
*** {prefix}_shortreads.vcf.gz.stats
*** {prefix}_longreads.vcf.gz
*** {prefix}_longreads.vcf.gz.stats
Both the short reads and long reads variant calling VCF is filtered for <code>QUAL > 20 </code>.
Freebayes (short read var calling) is ran with parameters <code>--use-best-n-alleles 4 --min-base-quality 10 --min-alternate-fraction 0.2 --haplotype-length 0 --ploidy 2 --min-alternate-count 2 </code>. For more details check the Snakefile.
** '''genome_alignment''' directory with results and figure from whole genome alignment
*** {prefix}_{species}.png
* '''mapped''' directory that contains the bam file with long reads mapped to the new assembly
** {prefix}_longreads.mapped.sorted.bam
* '''busco_{prefix}_before_polish_''' and '''busco_{prefix}_after_polish''' directories - contain busco results before and after polishing
** short_summary.specific.{lineage}.{prefix}_before_polish.txt
** short_summary.specific.{lineage}.{prefix}_after_polish.txt"
* '''other_files''' - directory containing other files created during the pipeline
* '''assembly''' - directory containing files created during the assembly step

Population variant calling pipeline

2022-06-13T14:34:13Z

Moiti001: added author info

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

For up-to-date documentation see [https://github.com/CarolinaPB/population-variant-calling here]

= Population level variant calling =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/population-variant-calling


== First follow the instructions here ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This is a pipeline that takes short reads aligned to a genome (in <code>.bam</code> format) and performs population level variant calling with <code>Freebayes</code>. It uses VEP to annotate the resulting VCF, calculates statistics, and calculates and plots a PCA.

It was developed to work with the results of [https://github.com/CarolinaPB/population-mapping this population mapping pipeline]. There are a few <code>Freebayes</code> requirements that you need to take into account if you don't use the mapping pipeline mentioned above to map your reads. You should make sure that:

* Alignments have read groups
* Alignments are sorted
* Duplicates are marked

See [https://github.com/freebayes/freebayes#calling-variants-from-fastq-to-vcf here] for more details.


==== Tools used ====

* [https://github.com/freebayes/freebayes Freebayes] - variant calling using short reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* [https://www.cog-genomics.org/plink/ Plink] - compute PCA
* R - Plot PCA

{|
!align="center"| [[File:population-var-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files ===

<syntaxhighlight lang="yaml">ASSEMBLY: /path/to/fasta
MAPPING_DIR: /path/to/bams/dir
PREFIX: <prefix>
OUTDIR: /path/to/outdir
SPECIES: <species>
NUM_CHRS: <number of chromosomes></syntaxhighlight>
* ASSEMBLY - path to genome fasta file
* MAPPING_DIR - path to directory with bam files to be used
** the pipeline will use all bam files in the directory, if you want to use a subset of those, create a file named <code>bam_list.txt</code> that contains the paths to the bam files you want to use. One path per line.

<syntaxhighlight lang="text">/path/to/file.bam
/path/to/file2.bam</syntaxhighlight>
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* SPECIES - species name to be used for VEP
* NUM_CHRS - number of chromosomes for your species (necessary for plink). ex: 38


== ADDITIONAL SET UP ==


=== Configuring VEP ===

This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:


==== For people using WUR's Anunna: ====

Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn't, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where "assembly name" is the name of the assembly you want to use.


==== For those not from WUR: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where "assembly name" is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>


=== Installing R packages ===

First load R: <code>module load R/3.6.2</code>

Enter the R environment by writing <code>R</code> and clicking enter. Install the packages:

<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>
If you get an error like this:

<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.


== RESULTS ==

The most important files and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''final_VCF''' directory with variant calling VCF files, as well as VCF stats
*** {prefix}.vep.vcf.gz - final VCF file
*** {prefix}.vep.vcf.gz.stats
** '''PCA''' PCA results and plot
*** {prefix}.eigenvec and {prefix}.eigenval - file with PCA eigenvectors and eigenvalues, respectively
*** {prefix}.pdf - PCA plot

The VCF file has been filtered for <code>QUAL > 20</code>. Freebayes is ran with parameters <code>--use-best-n-alleles 4 --min-base-quality 10 --min-alternate-fraction 0.2 --haplotype-length 0 --ploidy 2 --min-alternate-count 2</code>. These parameters can be changed in the Snakefile.

Bioinformatics tips tricks workflows

2022-06-13T14:25:46Z

Moiti001:

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling| Nanopore assembly and variant calling pipeline]]
* [[Population variant calling pipeline | Population variant calling pipeline]]
* [[Single Cell preprocessing pipeline| Single Cell preprocessing pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Bioinformatics tips tricks workflows

2022-06-13T14:21:46Z

Moiti001: added single cell page

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling| Nanopore assembly and variant calling pipeline]]
* [[Population variant calling pipeline | Population variant calling pipeline]]
* [[Single Cell preprocessing | Single Cell preprocessing]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Population variant calling pipeline

2022-03-02T14:18:57Z

Moiti001: added info on how to configure VEP and install R packages

= Population level variant calling =

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/population-variant-calling


== First follow the instructions here ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake


== ABOUT ==

This is a pipeline that takes short reads aligned to a genome (in <code>.bam</code> format) and performs population level variant calling with <code>Freebayes</code>. It uses VEP to annotate the resulting VCF, calculates statistics, and calculates and plots a PCA.

It was developed to work with the results of [https://github.com/CarolinaPB/population-mapping this population mapping pipeline]. There are a few <code>Freebayes</code> requirements that you need to take into account if you don't use the mapping pipeline mentioned above to map your reads. You should make sure that:

* Alignments have read groups
* Alignments are sorted
* Duplicates are marked

See [https://github.com/freebayes/freebayes#calling-variants-from-fastq-to-vcf here] for more details.


==== Tools used ====

* [https://github.com/freebayes/freebayes Freebayes] - variant calling using short reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* [https://www.cog-genomics.org/plink/ Plink] - compute PCA
* R - Plot PCA

{|
!align="center"| [[File:population-var-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}


=== Edit config.yaml with the paths to your files ===

<syntaxhighlight lang="yaml">ASSEMBLY: /path/to/fasta
MAPPING_DIR: /path/to/bams/dir
PREFIX: <prefix>
OUTDIR: /path/to/outdir
SPECIES: <species>
NUM_CHRS: <number of chromosomes></syntaxhighlight>
* ASSEMBLY - path to genome fasta file
* MAPPING_DIR - path to directory with bam files to be used
** the pipeline will use all bam files in the directory, if you want to use a subset of those, create a file named <code>bam_list.txt</code> that contains the paths to the bam files you want to use. One path per line.

<syntaxhighlight lang="text">/path/to/file.bam
/path/to/file2.bam</syntaxhighlight>
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* SPECIES - species name to be used for VEP
* NUM_CHRS - number of chromosomes for your species (necessary for plink). ex: 38


== ADDITIONAL SET UP ==


=== Configuring VEP ===

This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:


==== For people using WUR's Anunna: ====

Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn't, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where "assembly name" is the name of the assembly you want to use.


==== For those not from WUR: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where "assembly name" is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>


=== Installing R packages ===

First load R: <code>module load R/3.6.2</code>

Enter the R environment by writing <code>R</code> and clicking enter. Install the packages:

<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>
If you get an error like this:

<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.


== RESULTS ==

The most important files and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''final_VCF''' directory with variant calling VCF files, as well as VCF stats
*** {prefix}.vep.vcf.gz - final VCF file
*** {prefix}.vep.vcf.gz.stats
** '''PCA''' PCA results and plot
*** {prefix}.eigenvec and {prefix}.eigenval - file with PCA eigenvectors and eigenvalues, respectively
*** {prefix}.pdf - PCA plot

The VCF file has been filtered for <code>QUAL > 20</code>. Freebayes is ran with parameters <code>--use-best-n-alleles 4 --min-base-quality 10 --min-alternate-fraction 0.2 --haplotype-length 0 --ploidy 2 --min-alternate-count 2</code>. These parameters can be changed in the Snakefile.

File:Population-var-calling-workflow.png

2022-03-02T13:45:46Z

Moiti001:

Population variant calling pipeline

2022-03-02T13:45:28Z

Moiti001: added workflow image

Population variant calling pipeline

2022-03-02T13:43:18Z

Moiti001: added info on pipeline

Bioinformatics tips tricks workflows

2022-03-02T13:36:56Z

Moiti001: fixed nanopore assembly pipeline page link

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling| Nanopore assembly and variant calling pipeline]]
* [[Population variant calling pipeline | Population variant calling pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Bioinformatics tips tricks workflows

2022-03-02T13:33:36Z

Moiti001: added page for population variant calling pipeline

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling pipeline| Nanopore assembly and variant calling pipeline]]
* [[Population variant calling pipeline | Population variant calling pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Running Snakemake pipelines

2022-01-05T11:25:11Z

Moiti001: added info on how to install miniconda

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path 
In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp -r <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it
''Update 05/01/2022:'' 
Here I show how to install miniconda in a linux system 
[https://docs.conda.io/en/latest/miniconda.html Download installer] 
[https://conda.io/projects/conda/en/latest/user-guide/install/index.html Installation instructions]

# Download the installer to your home directory. Choose the version according to your operating system. You can right click the link, copy and download with

<pre>wget <link></pre>
At the time of writing this update, for me it would be:

<pre>wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh</pre>
To install miniconda, run:

<pre>bash <installer name></pre>
installer name could be <code>Miniconda3-latest-Linux-x86_64.sh</code>

Set up the conda channels in this order:

<pre>conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge</pre>

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>

You'll need to activate the conda environment again
<pre>conda activate <name-of-pipeline></pre>

Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see if there’s any mistake.

If all looks ok, you can now run your pipeline

<pre>snakemake --profile <name-of-pipeline></pre>
If everything was set up correctly, the jobs should be submitted and you should be able to see the progress of the pipeline in your terminal.

Mapping and variant calling pipeline

2021-11-22T11:56:38Z

Moiti001: info on how to concatenate fq files

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/mapping-variant-calling

[https://github.com/CarolinaPB/WUR_mapping-variant-calling Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads to a reference assembly. It outputs the mapped reads, a qualimap report and does variant calling.

==== Tools used: ====

* Bwa-mem2 - mapping
* Samtools - processing
* Qualimap - mapping summary
* Freebayes - variant calling
* Bcftools - VCF statistics

{|
!align="center"| [[File:mapping-variant-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* ASSEMBLY - path to the assembly file
* PREFIX - prefix for the final mapped reads file

If you want the results to be written to this directory (not to a new directory), comment out <pre>OUTDIR: /path/to/output</pre>

For the mapping step you should have one _1 fastq file and one _2 fastq file in <code>READS_DIR</code>. If you have several _1 and _2 fastq files from the same sample, you can combine them so you have one file for all _1 reads and one for all the _2 reads. This can be done by concatenating them using <code>cat</code>, if the original files are not compressed (<code>fastq</code> or <code>fq</code> extension), or <code>zcat</code> if the original files are compressed (<code>fastq.gz</code> or <code>fq.gz</code> extension). Example where your files are in the same directory and are compressed:

<pre>zcat *_1.fastq.gz > <new file name>_1.fastq.gz
zcat *_2.fastq.gz > <new file name>_2.fastq.gz</pre>

== RESULTS ==

* dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''sorted_reads''' directory with the file containing the mapped reads
* '''results''' directory containing the qualimap results
* '''variant_calling''' directory containing the variant calling VCF file and a file with with VCF statistics

Nanopore assembly and variant calling

2021-11-04T08:44:32Z

Moiti001: added info about variant calling filtering

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

= Assemble nanopore reads and do variant calling with short and long reads =
Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/nanopore-assembly

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline that uses <code>Flye</code> to create a nanopore assembly. It also does variant calling with long and short reads. 
The pipeline starts by using <code>porechop</code> to trim the adaptors, then it uses <code>Flye</code> to create the assembly. After that, <code>ntLink-arks</code> from <code>Lonstitch</code> is used to scaffold the assembly using the nanopore reads. The scaffolded assembly is polished with <code>polca</code>. <code>Bwa-mem2</code> is used to map the short reads to the assembly and <code>Freebayes</code> to do variant calling using these reads. <code>Minimap2</code> is used to map the long reads to the assembly, and <code>longshot</code> for variant calling using these. In the end, in addition to your assembly and variant calling results, you'll also get assembly statistics and busco scores before and after the polishing.

==== Tools used: ====

* [https://github.com/rrwick/Porechop Porechop] - trim adaptors
* [https://github.com/fenderglass/Flye Flye] - assembly
* [https://github.com/lh3/seqtk Seqtk] - convert fasta to one line fasta
* [https://github.com/bcgsc/longstitch LongStitch (ntLink-arks)] - scaffolding with nanopore reads
* [https://busco.ezlab.org/ BUSCO] - assess assembly completeness
* [https://github.com/alekseyzimin/masurca MaSuRCA (polca)] - polish assembly
* Python - get assembly stats
* [https://github.com/lh3/minimap2 Minimap2] - map long reads to reference. Genome alignment
* [http://www.htslib.org/ Samtools] - sort and index mapped reads and vcf files
* [https://github.com/pjedge/longshot Longshot] - variant calling with nanopore reads
* [https://github.com/bwa-mem2/bwa-mem2 Bwa-mem2] - map short reads to reference
* [https://github.com/freebayes/freebayes Freebayes] - variant calling using short reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* R - [https://github.com/tpoorten/dotPlotly pafCoordsDotPlotly] - plot genome alignment

{|
!align="center"| [[File:nanopore-assembly-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>LONGREADS: <nanopore_reads.fq.gz>
SHORTREADS:
- /path/to/short/reads_1.fq.gz
- /path/to/short/reads_2.fq.gz
GENOME_SIZE: <approximate genome size>
PREFIX: <prefix>
OUTDIR: /path/to/outdir
BUSCO_LINEAGE:
- <lineage>

# genome alignment parameters:
COMPARISON_GENOME:
<species>: /path/to/genome/fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000</pre>
* LONGREADS - name of file with long reads. This file should be in the working directory (where this config and the Snakefile are)
* SHORTREADS - paths to short reads fq.gz
* GENOME_SIZE - approximate genome size <code>haploid genome size (bp)(e.g. '3e9' for human genome)</code> from [https://github.com/bcgsc/longstitch#full-help-page longstitch]
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* BUSCO_LINEAGE - lineage used for busco. Can be one or more (one per line). To see available lineages run <code>busco --list-datasets</code>
* COMPARISON_GENOME - genome for whole genome comparison. Add your species name and the path to the fasta file. ex: <code>chicken: /path/to/chicken.fna.gz</code>. You can add several genomes, one on each line.
** If you don't want to run the genome alignment step, comment out

<pre>COMPARISON_GENOME:
<species>: /path/to/genome/fasta</pre>
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there's an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.

-----
If you have your long reads in several fastq files and need to create one file compressed file with all the reads:

# In your pipeline directory create one file with all the reads

<pre>cat /path/to/fastq/directory/*.fastq > <name of file>.fq</pre>
<ol start="2" style="list-style-type: decimal;">
<li>Compress the file you just created:</li></ol>

<pre>gzip <name of file>.fq</pre>

=== After installing the conda environment (first step of this guide) you'll need to edit the polca.sh file. ===

First go to the directory where miniconda3 is installed (usually your home directory). Go to <code>/<home>/miniconda/envs/<env_name>/bin</code> and open the file <code>polca.sh</code>. In my case the path looks like this: <code>/home/WUR/<username>/miniconda3/envs/<env_name>/bin/</code>. In your editor open <code>polca.sh</code> and replace this line:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) $BASM.alignSorted 2>>samtools.err && \</pre>
With this:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) -o $BASM.alignSorted.bam 2>>samtools.err && \</pre>
== RESULTS ==

The most important files are and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''{prefix}_oneline.k32.w100.ntLink-arks.longstitch-scaffolds.fa.PolcaCorrected.fa''' final assembly
** assembly_stats_<prefix>.txt file with assembly statistics for the final assembly
** '''variant_calling''' directory with variant calling VCF files with long and short reads, as well as VCF stats
*** {prefix}_shortreads.vcf.gz
*** {prefix}_shortreads.vcf.gz.stats
*** {prefix}_longreads.vcf.gz
*** {prefix}_longreads.vcf.gz.stats
Both the short reads and long reads variant calling VCF is filtered for <code>QUAL > 20 </code>.
Freebayes (short read var calling) is ran with parameters <code>--use-best-n-alleles 4 --min-base-quality 10 --min-alternate-fraction 0.2 --haplotype-length 0 --ploidy 2 --min-alternate-count 2 </code>. For more details check the Snakefile.
** '''genome_alignment''' directory with results and figure from whole genome alignment
*** {prefix}_{species}.png
* '''mapped''' directory that contains the bam file with long reads mapped to the new assembly
** {prefix}_longreads.mapped.sorted.bam
* '''busco_{prefix}_before_polish_''' and '''busco_{prefix}_after_polish''' directories - contain busco results before and after polishing
** short_summary.specific.{lineage}.{prefix}_before_polish.txt
** short_summary.specific.{lineage}.{prefix}_after_polish.txt"
* '''other_files''' - directory containing other files created during the pipeline
* '''assembly''' - directory containing files created during the assembly step

File:Nanopore-assembly-workflow.png

2021-11-01T10:03:49Z

Moiti001: Moiti001 uploaded a new version of File:Nanopore-assembly-workflow.png

Nanopore assembly and variant calling

2021-11-01T10:03:12Z

Moiti001: added info on var calling with short reads

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

= Assemble nanopore reads and do variant calling with short and long reads =
Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/nanopore-assembly

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline that uses <code>Flye</code> to create a nanopore assembly. It also does variant calling with long and short reads. 
The pipeline starts by using <code>porechop</code> to trim the adaptors, then it uses <code>Flye</code> to create the assembly. After that, <code>ntLink-arks</code> from <code>Lonstitch</code> is used to scaffold the assembly using the nanopore reads. The scaffolded assembly is polished with <code>polca</code>. <code>Bwa-mem2</code> is used to map the short reads to the assembly and <code>Freebayes</code> to do variant calling using these reads. <code>Minimap2</code> is used to map the long reads to the assembly, and <code>longshot</code> for variant calling using these. In the end, in addition to your assembly and variant calling results, you'll also get assembly statistics and busco scores before and after the polishing.

==== Tools used: ====

* [https://github.com/rrwick/Porechop Porechop] - trim adaptors
* [https://github.com/fenderglass/Flye Flye] - assembly
* [https://github.com/lh3/seqtk Seqtk] - convert fasta to one line fasta
* [https://github.com/bcgsc/longstitch LongStitch (ntLink-arks)] - scaffolding with nanopore reads
* [https://busco.ezlab.org/ BUSCO] - assess assembly completeness
* [https://github.com/alekseyzimin/masurca MaSuRCA (polca)] - polish assembly
* Python - get assembly stats
* [https://github.com/lh3/minimap2 Minimap2] - map long reads to reference. Genome alignment
* [http://www.htslib.org/ Samtools] - sort and index mapped reads and vcf files
* [https://github.com/pjedge/longshot Longshot] - variant calling with nanopore reads
* [https://github.com/bwa-mem2/bwa-mem2 Bwa-mem2] - map short reads to reference
* [https://github.com/freebayes/freebayes Freebayes] - variant calling using short reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* R - [https://github.com/tpoorten/dotPlotly pafCoordsDotPlotly] - plot genome alignment

{|
!align="center"| [[File:nanopore-assembly-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>LONGREADS: <nanopore_reads.fq.gz>
SHORTREADS:
- /path/to/short/reads_1.fq.gz
- /path/to/short/reads_2.fq.gz
GENOME_SIZE: <approximate genome size>
PREFIX: <prefix>
OUTDIR: /path/to/outdir
BUSCO_LINEAGE:
- <lineage>

# genome alignment parameters:
COMPARISON_GENOME:
<species>: /path/to/genome/fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000</pre>
* LONGREADS - name of file with long reads. This file should be in the working directory (where this config and the Snakefile are)
* SHORTREADS - paths to short reads fq.gz
* GENOME_SIZE - approximate genome size <code>haploid genome size (bp)(e.g. '3e9' for human genome)</code> from [https://github.com/bcgsc/longstitch#full-help-page longstitch]
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* BUSCO_LINEAGE - lineage used for busco. Can be one or more (one per line). To see available lineages run <code>busco --list-datasets</code>
* COMPARISON_GENOME - genome for whole genome comparison. Add your species name and the path to the fasta file. ex: <code>chicken: /path/to/chicken.fna.gz</code>. You can add several genomes, one on each line.
** If you don't want to run the genome alignment step, comment out

<pre>COMPARISON_GENOME:
<species>: /path/to/genome/fasta</pre>
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there's an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.

-----
If you have your long reads in several fastq files and need to create one file compressed file with all the reads:

# In your pipeline directory create one file with all the reads

<pre>cat /path/to/fastq/directory/*.fastq > <name of file>.fq</pre>
<ol start="2" style="list-style-type: decimal;">
<li>Compress the file you just created:</li></ol>

<pre>gzip <name of file>.fq</pre>

=== After installing the conda environment (first step of this guide) you'll need to edit the polca.sh file. ===

First go to the directory where miniconda3 is installed (usually your home directory). Go to <code>/<home>/miniconda/envs/<env_name>/bin</code> and open the file <code>polca.sh</code>. In my case the path looks like this: <code>/home/WUR/<username>/miniconda3/envs/<env_name>/bin/</code>. In your editor open <code>polca.sh</code> and replace this line:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) $BASM.alignSorted 2>>samtools.err && \</pre>
With this:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) -o $BASM.alignSorted.bam 2>>samtools.err && \</pre>
== RESULTS ==

The most important files are and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''{prefix}_oneline.k32.w100.ntLink-arks.longstitch-scaffolds.fa.PolcaCorrected.fa''' final assembly
** assembly_stats_<prefix>.txt file with assembly statistics for the final assembly
** '''variant_calling''' directory with variant calling VCF files with long and short reads, as well as VCF stats
*** {prefix}_shortreads.vcf.gz
*** {prefix}_shortreads.vcf.gz.stats
*** {prefix}_longreads.vcf.gz
*** {prefix}_longreads.vcf.gz.stats
** '''genome_alignment''' directory with results and figure from whole genome alignment
*** {prefix}_{species}.png
* '''mapped''' directory that contains the bam file with long reads mapped to the new assembly
** {prefix}_longreads.mapped.sorted.bam
* '''busco_{prefix}_before_polish_''' and '''busco_{prefix}_after_polish''' directories - contain busco results before and after polishing
** short_summary.specific.{lineage}.{prefix}_before_polish.txt
** short_summary.specific.{lineage}.{prefix}_after_polish.txt"
* '''other_files''' - directory containing other files created during the pipeline
* '''assembly''' - directory containing files created during the assembly step

Mapping and variant calling pipeline

2021-11-01T09:38:22Z

Moiti001: changed mapping tool

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/mapping-variant-calling

[https://github.com/CarolinaPB/WUR_mapping-variant-calling Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads to a reference assembly. It outputs the mapped reads, a qualimap report and does variant calling.

==== Tools used: ====

* Bwa-mem2 - mapping
* Samtools - processing
* Qualimap - mapping summary
* Freebayes - variant calling
* Bcftools - VCF statistics

{|
!align="center"| [[File:mapping-variant-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* ASSEMBLY - path to the assembly file
* PREFIX - prefix for the final mapped reads file

If you want the results to be written to this directory (not to a new directory), comment out <pre>READS_DIR: /path/to/reads/</pre>
== RESULTS ==

* dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''sorted_reads''' directory with the file containing the mapped reads
* '''results''' directory containing the qualimap results
* '''variant_calling''' directory containing the variant calling VCF file and a file with with VCF statistics

File:Nanopore-assembly-workflow.png

2021-10-29T08:38:21Z

Moiti001:

Nanopore assembly and variant calling

2021-10-29T08:37:36Z

Moiti001: added info about nanopore pipeline

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

= Assemble nanopore reads and do variant calling with short and long reads =
Path to pipeline: /lustre/nobackup/WUR/ABGC/shared/PIPELINES/nanopore-assembly

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline that uses <code>Flye</code> to create a nanopore assembly. It also does variant calling with long and short reads. 
The pipeline starts by using <code>porechop</code> to trim the adaptors, then it uses <code>Flye</code> to create the assembly. After that, <code>ntLink-arks</code> from <code>Lonstitch</code> is used to scaffold the assembly using the nanopore reads. The scaffolded assembly is polished with <code>polca</code>. <code>Polca</code> also does variant calling with the short reads, while <code>longshot</code> does variant calling with the nanopore reads. To run <code>longshot</code>, first the long reads are aligned to the assembly with <code>minimap2</code>. 
In the end, in addition to your assembly and variant calling results, you'll also get assembly statistics and busco scores before and after the polishing.

==== Tools used: ====

* [https://github.com/rrwick/Porechop Porechop] - trim adaptors
* [https://github.com/fenderglass/Flye Flye] - assembly
* [https://github.com/lh3/seqtk Seqtk] - convert fasta to one line fasta
* [https://github.com/bcgsc/longstitch LongStitch (ntLink-arks)] - scaffolding with nanopore reads
* [https://busco.ezlab.org/ BUSCO] - assess assembly completeness
* [https://github.com/alekseyzimin/masurca MaSuRCA (polca)] - polish assembly and do variant calling with short reads
* Python - get assembly stats
* [https://github.com/lh3/minimap2 Minimap2] - map long reads to reference. Genome alignment
* [http://www.htslib.org/ Samtools] - sort and index mapped reads and vcf files
* [https://github.com/pjedge/longshot Longshot] - variant calling with nanopore reads
* [https://samtools.github.io/bcftools/bcftools.html bcftools] - vcf statistics
* R - [https://github.com/tpoorten/dotPlotly pafCoordsDotPlotly] - plot genome alignment

{|
!align="center"| [[File:nanopore-assembly-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>LONGREADS: <nanopore_reads.fq.gz>
SHORTREADS:
- /path/to/short/reads_1.fq.gz
- /path/to/short/reads_2.fq.gz
GENOME_SIZE: <approximate genome size>
PREFIX: <prefix>
OUTDIR: /path/to/outdir
BUSCO_LINEAGE:
- <lineage>

# genome alignment parameters:
COMPARISON_GENOME:
<species>: /path/to/genome/fasta

# filter alignments less than cutoff X bp
MIN_ALIGNMENT_LENGTH: 10000
MIN_QUERY_LENGTH: 50000</pre>
* LONGREADS - name of file with long reads. This file should be in the working directory (where this config and the Snakefile are)
* SHORTREADS - paths to short reads fq.gz
* GENOME_SIZE - approximate genome size <code>haploid genome size (bp)(e.g. '3e9' for human genome)</code> from [https://github.com/bcgsc/longstitch#full-help-page longstitch]
* PREFIX - prefix for the created files
* OUTDIR - directory where snakemake will run and where the results will be written to 
If you want the results to be written to this directory (not to a new directory), open config.yaml and comment out <code>OUTDIR: /path/to/outdir</code>
* BUSCO_LINEAGE - lineage used for busco. Can be one or more (one per line). To see available lineages run <code>busco --list-datasets</code>
* COMPARISON_GENOME - genome for whole genome comparison. Add your species name and the path to the fasta file. ex: <code>chicken: /path/to/chicken.fna.gz</code>. You can add several genomes, one on each line.
** If you don't want to run the genome alignment step, comment out

<pre>COMPARISON_GENOME:
<species>: /path/to/genome/fasta</pre>
* MIN_ALIGNMENT_LENGTH and MIN_QUERY_LENGTH - parameters for plotting. If your plot is coming out blank or if there's an error with the plotting step, try lowering these thresholds. This happens because the alignments are not large enough.

-----
If you have your long reads in several fastq files and need to create one file compressed file with all the reads:

# In your pipeline directory create one file with all the reads

<pre>cat /path/to/fastq/directory/*.fastq > <name of file>.fq</pre>
<ol start="2" style="list-style-type: decimal;">
<li>Compress the file you just created:</li></ol>

<pre>gzip <name of file>.fq</pre>

=== After installing the conda environment (first step of this guide) you'll need to edit the polca.sh file. ===

First go to the directory where miniconda3 is installed (usually your home directory). Go to <code>/<home>/miniconda/envs/<env_name>/bin</code> and open the file <code>polca.sh</code>. In my case the path looks like this: <code>/home/WUR/<username>/miniconda3/envs/<env_name>/bin/</code>. In your editor open <code>polca.sh</code> and replace this line:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) $BASM.alignSorted 2>>samtools.err && \</pre>
With this:

<pre>$SAMTOOLS sort -m $MEM -@ $NUM_THREADS <(samtools view -uhS $BASM.unSorted.sam) -o $BASM.alignSorted.bam 2>>samtools.err && \</pre>
== RESULTS ==

The most important files are and directories are:

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''results''' directory that contains
** '''{prefix}_oneline.k32.w100.ntLink-arks.longstitch-scaffolds.fa.PolcaCorrected.fa''' final assembly
** assembly_stats_<prefix>.txt file with assembly statistics for the final assembly
** '''variant_calling''' directory with variant calling VCF files with long and short reads, as well as VCF stats
*** {prefix}_shortreads.vcf.gz
*** {prefix}_shortreads.vcf.gz.stats
*** {prefix}_longreads.vcf.gz
*** {prefix}_longreads.vcf.gz.stats
** '''genome_alignment''' directory with results and figure from whole genome alignment
*** {prefix}_{species}.png
* '''mapped''' directory that contains the bam file with long reads mapped to the new assembly
** {prefix}_longreads.mapped.sorted.bam
* '''busco_{prefix}_before_polish_''' and '''busco_{prefix}_after_polish''' directories - contain busco results before and after polishing
** short_summary.specific.{lineage}.{prefix}_before_polish.txt
** short_summary.specific.{lineage}.{prefix}_after_polish.txt"
* '''other_files''' - directory containing other files created during the pipeline
* '''assembly''' - directory containing files created during the assembly step

Population structural variant calling pipeline

2021-10-28T08:50:29Z

Moiti001: added info about split read support option

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python
** PyVcf - add depth to vcf and create final table
** bamgroupreads.py + samblaster - create bam files with split and discordant reads

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name>
NUM_CHRS: <number of chromosomes>
BWA_MEM_M: <Y/N></pre>

* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example:
<blockquote>sample1,sample1.bam,Pop1 
sample2,sample2.bam,Pop1 
sample3,sample3.bam,Pop2 
sample4,sample4.bam,Pop2 </blockquote>
Tip: use the name of the bam file without the .bam extension as the sample name. Ex: from sample1.bam to sample1
* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files
* NUM_CHRS - number of chromosomes for your species (necessary for plink). ex: 38
* BWA_MEM_M - if you mapped your reads with `bwa mem` using the `-M` parameter and you want split read support in your VCF you need to run an extra step. For this write `Y`.
For a more detailed explanation see here [https://carolinapb.github.io/2021-10-28-smoove-SR-support/ Smoove SR support]

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==
=== Configuring VEP ===
This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

=== Installing R packages ===
First load R:
<pre>module load R/3.6.2</pre>
Enter the R environment by writing <pre>R</pre> and clicking enter.
Install the packages:
<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>

If you get an error like this:
<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF - with VEP annotation, not filtered for quality
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information. The variants have been filtered with Minimum Quality score = 30
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions. It has been filtered with Minimum Quality score = 30 and the DEPTH* field was added
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

-----

*The '''DEPTH''' field in the vcf has six fields, corresponding to the average depth across all samples.

<pre>DEPTH=(DHBFC_1/1, DHBFC_0/1, DHBFC_0/0, DHFFC_1/1, DHFFC_0/1, DHFFC_0/0)</pre>
Depth fold change as described in [https://github.com/brentp/duphold duphold]:
<blockquote>DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
DHFFC: fold-change for the variant depth relative to Flanking regions.</blockquote>

These fields are also in the <code>{prefix}_DUP_DEL_INV_table.tsv</code> file

File:Pop-sv-calling-workflow.png

2021-10-28T08:43:27Z

Moiti001: Moiti001 uploaded a new version of File:Pop-sv-calling-workflow.png

population structural variation calling pipeline workflow

File:Mapping-variant-calling-workflow.png

2021-10-22T11:48:00Z

Moiti001: workflow for mapping and variant calling pipeline

workflow for mapping and variant calling pipeline

Mapping and variant calling pipeline

2021-10-22T11:46:42Z

Moiti001: added pipeline workflow, small details

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/WUR_mapping-variant-calling Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads to a reference assembly. It outputs the mapped reads, a qualimap report and does variant calling.

==== Tools used: ====

* Bwa - mapping
* Samtools - processing
* Qualimap - mapping summary
* Freebayes - variant calling
* Bcftools - VCF statistics

{|
!align="center"| [[File:mapping-variant-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* ASSEMBLY - path to the assembly file
* PREFIX - prefix for the final mapped reads file

If you want the results to be written to this directory (not to a new directory), comment out <pre>READS_DIR: /path/to/reads/</pre>
== RESULTS ==

* dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''sorted_reads''' directory with the file containing the mapped reads
* '''results''' directory containing the qualimap results
* '''variant_calling''' directory containing the variant calling VCF file and a file with with VCF statistics

Population structural variant calling pipeline

2021-10-15T11:49:16Z

Moiti001: added new config parameter - num of chrs

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python - add depth to vcf and create final table
** PyVcf

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name>
NUM_CHRS: <number of chromosomes></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example:
<blockquote>sample1,sample1.bam,Pop1 
sample2,sample2.bam,Pop1 
sample3,sample3.bam,Pop2 
sample4,sample4.bam,Pop2 </blockquote>
Tip: use the name of the bam file without the .bam extension as the sample name. Ex: from sample1.bam to sample1
* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files
* NUM_CHRS - number of chromosomes for your species (necessary for plink). ex: 38

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==
=== Configuring VEP ===
This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

=== Installing R packages ===
First load R:
<pre>module load R/3.6.2</pre>
Enter the R environment by writing <pre>R</pre> and clicking enter.
Install the packages:
<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>

If you get an error like this:
<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF - with VEP annotation, not filtered for quality
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information. The variants have been filtered with Minimum Quality score = 30
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions. It has been filtered with Minimum Quality score = 30 and the DEPTH* field was added
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

-----

*The '''DEPTH''' field in the vcf has six fields, corresponding to the average depth across all samples.

<pre>DEPTH=(DHBFC_1/1, DHBFC_0/1, DHBFC_0/0, DHFFC_1/1, DHFFC_0/1, DHFFC_0/0)</pre>
Depth fold change as described in [https://github.com/brentp/duphold duphold]:
<blockquote>DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
DHFFC: fold-change for the variant depth relative to Flanking regions.</blockquote>

These fields are also in the <code>{prefix}_DUP_DEL_INV_table.tsv</code> file

Bioinformatics tips tricks workflows

2021-10-15T09:15:25Z

Moiti001: added nanopore assembly pipeline page

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]
* [[Nanopore assembly and variant calling | Nanopore assembly and variant calling]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Running Snakemake pipelines

2021-10-08T12:57:02Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path 
In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp -r <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>

You'll need to activate the conda environment again
<pre>conda activate <name-of-pipeline></pre>

Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see if there’s any mistake.

If all looks ok, you can now run your pipeline

<pre>snakemake --profile <name-of-pipeline></pre>
If everything was set up correctly, the jobs should be submitted and you should be able to see the progress of the pipeline in your terminal.

Population structural variant calling pipeline

2021-10-08T12:52:38Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python - add depth to vcf and create final table
** PyVcf

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example:
<blockquote>sample1,sample1.bam,Pop1 
sample2,sample2.bam,Pop1 
sample3,sample3.bam,Pop2 
sample4,sample4.bam,Pop2 </blockquote>
Tip: use the name of the bam file without the .bam extension as the sample name. Ex: from sample1.bam to sample1
* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==
=== Configuring VEP ===
This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

=== Installing R packages ===
First load R:
<pre>module load R/3.6.2</pre>
Enter the R environment by writing <pre>R</pre> and clicking enter.
Install the packages:
<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>

If you get an error like this:
<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here] and try to install the packages again.

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF - with VEP annotation, not filtered for quality
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information. The variants have been filtered with Minimum Quality score = 30
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions. It has been filtered with Minimum Quality score = 30 and the DEPTH* field was added
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

-----

*The '''DEPTH''' field in the vcf has six fields, corresponding to the average depth across all samples.

<pre>DEPTH=(DHBFC_1/1, DHBFC_0/1, DHBFC_0/0, DHFFC_1/1, DHFFC_0/1, DHFFC_0/0)</pre>
Depth fold change as described in [https://github.com/brentp/duphold duphold]:
<blockquote>DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
DHFFC: fold-change for the variant depth relative to Flanking regions.</blockquote>

These fields are also in the <code>{prefix}_DUP_DEL_INV_table.tsv</code> file

Population structural variant calling pipeline

2021-10-08T09:42:37Z

Moiti001: added tip on how to define sample name for sample list file

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python - add depth to vcf and create final table
** PyVcf

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example:
<blockquote>sample1,sample1.bam,Pop1 
sample2,sample2.bam,Pop1 
sample3,sample3.bam,Pop2 
sample4,sample4.bam,Pop2 </blockquote>
Tip: use the name of the bam file without the .bam extension as the sample name. Ex: from sample1.bam to sample1
* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==
=== Configuring VEP ===
This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

=== Installing R packages ===
First load R:
<pre>module load R/3.6.2</pre>
Enter the R environment by writing <pre>R</pre> and clicking enter.
Install the packages:
<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>

If you get an error like this:
<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here]

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

Population structural variant calling pipeline

2021-10-08T09:20:57Z

Moiti001: added info about installing R packages

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python - add depth to vcf and create final table
** PyVcf

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example: > sample1,sample1.bam,Pop1 
> sample2,sample2.bam,Pop1 
> sample3,sample3.bam,Pop2 
> sample4,sample4.bam,Pop2 

* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==
=== Configuring VEP ===
This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

=== Installing R packages ===
First load R:
<pre>module load R/3.6.2</pre>
Enter the R environment by writing <pre>R</pre> and clicking enter.
Install the packages:
<pre>list.of.packages <- c("optparse", "data.table", "ggplot2")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)</pre>

If you get an error like this:
<pre>Warning in install.packages(new.packages) :
'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>
Follow the instructions on how to install R packages locally [https://wiki.anunna.wur.nl/index.php/Installing_R_packages_locally here]

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

Population mapping pipeline

2021-09-29T12:55:31Z

Moiti001: added workflow image

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-mapping Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads from several individuals to a reference assembly. It outputs the mapped reads and a qualimap report.

==== Tools used: ====

* Bwa - mapping
* Samtools - processing
* Qualimap - mapping summary

{|
!align="center"| [[File:Population-mapping-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>ASSEMBLY: /path/to/assembly
OUTDIR: /path/to/outdir
PATHS_WITH_FILES:
path1: /path/to/dir</pre>
* ASSEMBLY - path to the assembly file
* OUTDIR - directory where snakemake will run and where the results will be written to. 
If you want the results to be written to this directory (not to a new directory), comment out <code>OUTDIR: /path/to/outdir</code>
* PATHS_WITH_FILES - directory that can contain subdirectories where the '''fq.gz''' reads are located. You can add several paths by adding <code>path2: /path/to/dir</code> under <code>PATHS_WITH_FILES</code>. (The line you add has to have indentation)

The script goes through the subdirectories of the directory you choose under <code>PATHS_WITH_FILES</code> looking for files with '''fq.gz''' extension. 
Example: if <code>path1: /lustre/nobackup/WUR/ABGC/shared/Chicken/Africa/X201SC20031230-Z01-F006_multipath</code>, the subdirectory structure could be:

<pre>/lustre/nobackup/WUR/ABGC/shared/Chicken/Africa/X201SC20031230-Z01-F006_multipath
├── X201SC20031230-Z01-F006_1
│ └── raw_data
│ ├── a109_26_15_1_H
│ │ ├── a109_26_15_1_H_FDSW202597655-1r_HWFFFDSXY_L3_1.fq.gz
│ │ ├── a109_26_15_1_H_FDSW202597655-1r_HWFFFDSXY_L3_2.fq.gz
│ │ └── MD5.txt
│ └── a20_10_16_1_H
│ ├── a20_10_16_1_H_FDSW202597566-1r_HWFFFDSXY_L3_1.fq.gz
│ ├── a20_10_16_1_H_FDSW202597566-1r_HWFFFDSXY_L3_2.fq.gz
│ └── MD5.txt
└── X201SC20031230-Z01-F006_2
└── raw_data
├── a349_Be_17_1_C
│ ├── a349_Be_17_1_C_FDSW202597895-1r_HWFFFDSXY_L3_1.fq.gz
│ ├── a349_Be_17_1_C_FDSW202597895-1r_HWFFFDSXY_L3_2.fq.gz
│ └── MD5.txt
└── a360_Be_05_1_H
├── a360_Be_05_1_H_FDSW202597906-1r_HWFFFDSXY_L3_1.fq.gz
├── a360_Be_05_1_H_FDSW202597906-1r_HWFFFDSXY_L3_2.fq.gz
└── MD5.txt </pre>
== RESULTS ==

* '''<run_date>_files.txt''' dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''processed_reads''' directory with the bam files with the mapped reads for every sample
* '''mapping_stats''' directory containing the qualimap results and a summary of the qualimap results for all samples in <code>sample_quality_summary.tsv</code>
** '''qualimap''' contains qualimap results per sample

File:Population-mapping-workflow.png

2021-09-29T12:51:44Z

Moiti001: workflow for population mapping snakemake pipeline

workflow for population mapping snakemake pipeline

Population mapping pipeline

2021-09-29T12:49:30Z

Moiti001: added info about pipeline

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-mapping Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads from several individuals to a reference assembly. It outputs the mapped reads and a qualimap report.

==== Tools used: ====

* Bwa - mapping
* Samtools - processing
* Qualimap - mapping summary

{|
!align="center"| [[File:https://github.com/CarolinaPB/population-mapping/blob/wur/workflow.png|DAG]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>ASSEMBLY: /path/to/assembly
OUTDIR: /path/to/outdir
PATHS_WITH_FILES:
path1: /path/to/dir</pre>
* ASSEMBLY - path to the assembly file
* OUTDIR - directory where snakemake will run and where the results will be written to. 
If you want the results to be written to this directory (not to a new directory), comment out <code>OUTDIR: /path/to/outdir</code>
* PATHS_WITH_FILES - directory that can contain subdirectories where the '''fq.gz''' reads are located. You can add several paths by adding <code>path2: /path/to/dir</code> under <code>PATHS_WITH_FILES</code>. (The line you add has to have indentation)

The script goes through the subdirectories of the directory you choose under <code>PATHS_WITH_FILES</code> looking for files with '''fq.gz''' extension. 
Example: if <code>path1: /lustre/nobackup/WUR/ABGC/shared/Chicken/Africa/X201SC20031230-Z01-F006_multipath</code>, the subdirectory structure could be:

<pre>/lustre/nobackup/WUR/ABGC/shared/Chicken/Africa/X201SC20031230-Z01-F006_multipath
├── X201SC20031230-Z01-F006_1
│ └── raw_data
│ ├── a109_26_15_1_H
│ │ ├── a109_26_15_1_H_FDSW202597655-1r_HWFFFDSXY_L3_1.fq.gz
│ │ ├── a109_26_15_1_H_FDSW202597655-1r_HWFFFDSXY_L3_2.fq.gz
│ │ └── MD5.txt
│ └── a20_10_16_1_H
│ ├── a20_10_16_1_H_FDSW202597566-1r_HWFFFDSXY_L3_1.fq.gz
│ ├── a20_10_16_1_H_FDSW202597566-1r_HWFFFDSXY_L3_2.fq.gz
│ └── MD5.txt
└── X201SC20031230-Z01-F006_2
└── raw_data
├── a349_Be_17_1_C
│ ├── a349_Be_17_1_C_FDSW202597895-1r_HWFFFDSXY_L3_1.fq.gz
│ ├── a349_Be_17_1_C_FDSW202597895-1r_HWFFFDSXY_L3_2.fq.gz
│ └── MD5.txt
└── a360_Be_05_1_H
├── a360_Be_05_1_H_FDSW202597906-1r_HWFFFDSXY_L3_1.fq.gz
├── a360_Be_05_1_H_FDSW202597906-1r_HWFFFDSXY_L3_2.fq.gz
└── MD5.txt </pre>
== RESULTS ==

* **<run_date>_files.txt** dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''processed_reads''' directory with the bam files with the mapped reads for every sample
* '''mapping_stats''' directory containing the qualimap results and a summary of the qualimap results for all samples in <code>sample_quality_summary.tsv</code>
** '''qualimap''' contains qualimap results per sample

Bioinformatics tips tricks workflows

2021-09-29T12:42:24Z

Moiti001: added page for population mapping pipeline

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]
* [[Population mapping pipeline | Population mapping pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Population structural variant calling pipeline

2021-09-22T11:39:24Z

Moiti001: added info about population SV calling pipeline

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/population-structural-var-calling-smoove/tree/single_run Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to perform structural variant calling in a population using Smoove. It also runs VEP and performs PCA. In addition to the VCF with the SVs, you also get a .tsv file with some summarized information on the SVs: it includes allele frequency per population, as well as VEP annotation and depth fold change as described in [https://github.com/brentp/duphold duphold]: 
> DHBFC: fold-change for the variant depth relative to bins in the genome with similar GC-content. 
> DHFFC: fold-change for the variant depth relative to Flanking regions.

==== Tools used: ====

* Smoove - SV calling
* VEP - determines the effect of the variants
* Plink - perform PCA
* R - plot PCA
* SURVIVOR - basic SV stats
* Python - add depth to vcf and create final table
** PyVcf

{|
!align="center"| [[File:Pop-sv-calling-workflow.png]]
|-
|align="center"| ''Pipeline workflow''
|}

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
SAMPLE_LIST: /path/to/file
REFERENCE: /path/to/assembly
CONTIGS_IGNORE: /path/to/file
SPECIES: <species_name>
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* SAMPLE_LIST - three column csv with the sample name, name of the bam files to use in the second column and the name of the corresponding population on the third column. These bams should all be in the same directory (READS_DIR)
* Example: > sample1,sample1.bam,Pop1 
> sample2,sample2.bam,Pop1 
> sample3,sample3.bam,Pop2 
> sample4,sample4.bam,Pop2 

* REFERENCE - path to the assembly file
* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)
** If you don't want to exclude contigs you'll need to edit the Snakefile to remove this line <code>--excludechroms {params.contigs} \</code>
* SPECIES - species name to be used for VEP
* PREFIX - prefix for the created files

If you want the results to be written to this directory (not to a new directory), comment out or remove

<pre>OUTDIR: /path/to/outdir</pre>
== ADDITIONAL SET UP ==

This pipeline uses VEP in offline mode, which increases performance. In order to use it in this mode, the cache for the species used needs to be installed:
Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with

<pre>/usr/bin/perl /cm/shared/apps/SHARED/ensembl-vep/INSTALL.pl --CACHEDIR /lustre/nobackup/SHARED/cache/ --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

==== Other option: ====

You can install VEP with

<pre>conda install -c bioconda ensembl-vep</pre>
and install the cache with

<pre>vep_install --CACHEDIR <where/to/install/cache> --AUTO c -n --SPECIES <species></pre>
When multiple assemblies are found you need to run it again with <code>--ASSEMBLY <assembly name></code>, where “assembly name” is the name of the assembly you want to use.

In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

== RESULTS ==

* '''<run_date>_files.txt''' Dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''2_merged'''
** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering
* '''5_postprocessing''' directory that contains the final VCF file containing the structural variants found. This file has been annotated with VEP
** {prefix}.smoove.square.vep.vcf.gz - Final VCF
** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP
** {prefix}.nosex, {prefix}.log, {prefix}.eigenvec, {prefix}.eigenval - output files from the PCA
** {prefix}_DUP_DEL_INV_table.tsv - table with the most important information extracted from the VCF. Contains information about the SV, allele frequency for each population, VEP annotation and depth information
** {prefix}_DUP_DEL_INV.vcf - vcf file with annotated duplications, deletions and inversions
** {prefix}_BND.vcf - vcf file with variants annotated with BND
* '''6_metrics''' directory that contains general stats about the number of SVs found
* '''FIGURES''' directory that contains the PCA plot

What you do with the results from this structural variant calling pipeline depends on your research question: a possible next step would be to explore the '''{prefix}_DUP_DEL_INV_table.tsv''' file and look at the largest SVs found (sort by ''SVLEN'') or at a specific effect in the ANNOTATION column, such as “frameshift_variant”.

See [https://m.ensembl.org/info/genome/variation/prediction/predicted_data.html VEP effect descriptions] for a short description of the effects annotated by VEP.

File:Pop-sv-calling-workflow.png

2021-09-22T11:30:39Z

Moiti001: population structural variation calling pipeline workflow

population structural variation calling pipeline workflow

Bioinformatics tips tricks workflows

2021-09-22T11:19:16Z

Moiti001: /* A list of tutorials, workflows, and recipes */

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]
* [[Population structural variant calling pipeline | Population structural variant calling pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Running Snakemake pipelines

2021-07-01T13:38:17Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path 
In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>

You'll need to activate the conda environment again
<pre>conda activate <name-of-pipeline></pre>

Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see if there’s any mistake.

If all looks ok, you can now run your pipeline

<pre>snakemake --profile <name-of-pipeline></pre>
If everything was set up correctly, the jobs should be submitted and you should be able to see the progress of the pipeline in your terminal.

Running Snakemake pipelines

2021-07-01T08:03:10Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path 
In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>
Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see if there’s any mistake.

If all looks ok, you can now run your pipeline

<pre>snakemake --profile <name-of-pipeline></pre>
If everything was set up correctly, the jobs should be submitted and you should be able to see the progress of the pipeline in your terminal.

Running Snakemake pipelines

2021-07-01T06:39:30Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>
Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see if there’s any mistake.

If all looks ok, you can now run your pipeline

<pre>snakemake --profile <name-of-pipeline></pre>
If everything was set up correctly, the jobs should be submitted and you should be able to see the progress of the pipeline in your terminal.

Running Snakemake pipelines

2021-06-30T14:24:33Z

Moiti001: Created page with "Author: Carolina Pita Barros Contact: carolina.pitabarros@wur.nl ABG You can find my pipelines [https://github.com/CarolinaPB/ here] The Snakemake sha..."

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG

 
You can find my pipelines [https://github.com/CarolinaPB/ here]

The Snakemake shared here use modules loaded from the HPC and tools installed with conda.

Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== Clone the repository ==

==== From github ====

Go to the repository’s page, click the green “Code” button and copy the path In your terminal go to where you want to download it to and run

<pre>git clone <path you copied from github></pre>
==== From the the WUR HPC (Anunna) ====

Go to <code>/lustre/nobackup/WUR/ABGC/shared/PIPELINES/</code> and choose which pipeline you want to use.

<pre>cp <pipeline directory> <directory where you want to save it to></pre>
First you’ll need to do some set up. Go to the pipeline’s directory.

== Installation ==

Install <code>conda</code> if you don’t have it

=== Create conda environment ===

<pre>conda create --name <name-of-pipeline> --file requirements.txt</pre>
<blockquote>I recommend giving it the same name as the pipeline
</blockquote>
This environment contains snakemake and the other packages that are needed to run the pipeline.

=== Activate environment ===

<pre>conda activate <name-of-pipeline></pre>
=== To deactivate the environment (if you want to leave the conda environment) ===

<pre>conda deactivate</pre>
== File configuration ==

=== Create HPC config file ===

Necessary for snakemake to prepare and send jobs.

==== Start with creating the directory ====

<pre>mkdir -p ~/.config/snakemake/<name-of-pipeline>
cd ~/.config/snakemake/<name-of-pipeline></pre>
==== Create config.yaml and include the following: ====

<blockquote>My pipelines are configured to work with SLURM
</blockquote>
<pre>jobs: 10
cluster: "sbatch -t 1:0:0 --mem=16000 -c 16 --job-name={rule} --exclude=fat001,fat002,fat101,fat100 --output=logs_slurm/{rule}.out --error=logs_slurm/{rule}.err"

use-conda: true</pre>
<blockquote>Here you should configure the resources you want to use.
</blockquote>
=== Go to the pipeline directory and open config.yaml ===

Configure your paths, but keep the variable names that are already in the config file.

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
If you want the results to be written to this directory (not to a new directory), open the Snakefile and comment out <code>workdir: config["OUTDIR"]</code> and ignore or comment out the <code>OUTDIR: /path/to/output</code> in the config file.

'''Now the setup is complete'''

== How to run the pipeline ==

Since the pipelines can take a while to run, it’s best if you use a [https://linuxize.com/post/how-to-use-linux-screen/ screen session]. By using a screen session, Snakemake stays “active” in the shell while it’s running, there’s no risk of the connection going down and Snakemake stopping.

Start by creating a screen session:

<pre>screen -S <name of session></pre>
Then run

<pre>snakemake -np</pre>
This will show you the steps and commands that will be executed. Check the commands and file names to see i

Bioinformatics tips tricks workflows

2021-06-30T14:19:28Z

Moiti001:

This page is intended as a portal to pages concerning best practices, workflows and pipelines, and other protocols (including scripts).

== A list of tutorials, workflows, and recipes ==
* [[Mapping_reads_with_Mosaik | Mapping Illumina GA2/HiSeq reads to the Sus scrofa genome assembly]]
* [[convert_fastq_to_fasta | A Perl script to convert fastq to fasta file format]]
* [[Mapping Pair-end reads with Stampy]]
* [[making_slices_from_BAM_files | Create slices from a collection of BAM files ]]
* [[Setting_up_Python_virtualenv | Setting up and using a virtual environment for Python3 ]]
* [[ssh_without_password | ssh without password]]
* [[Create_shortcut_log-in_command | Create a shortcut for the ssh log-in command]]
* [[Installing_R_packages_locally | Installing R packages locally]]
* [[command_line_tricks_for_manipulating_fastq | Command-line tricks for manipulating fastq files]]
* [[assemble_mitochondrial_genomes_from_short_read_data | Assemble mitochondrial genomes from whole-genome short-read data]]
* [[1000Bulls_mapping_pipeline_at_ABGC | 1000 Bulls mapping pipeline at ABGC]]
* [[ABGSA | Animal Breeding and Genomics Sequence Archives (ABGSA)]]
* [[Short_read_mapping_pipeline_pig | Pig mapping pipeline at ABGC]]
* [[Extract_noncall_snps_from_soy | Extract a set of pig SNPs not called in a control sample (soybean)]]
* [[calculate_corrected_theta_from_resequencing_data | Calculate nucleotide diversity (theta) corrected for sequencing depth]]
* [[RNA-seq analysis | RNA-seq analysis with Tophat]]
* [[Variant_annotation_tutorial | Variant annotation tutorial]]
* [[issues_asreml | Issues with ASReml]]
* [[Checkpointing | Checkpointing]]
* [[Assembly & Annotation | Assembly and Annotation guidelines (denovo)]]
* [[DE expression | DE expression analysis with tophat2 / cuffdiff]]
* [[JBrowse | JBrowse]]
* [[Running Snakemake pipelines | Running Snakemake pipelines]]
* [[Mapping and variant calling pipeline | Mapping and variant calling pipeline]]

== External links ==
* [http://en.wikipedia.org/wiki/Help:Cheatsheet Help with editing Wiki pages]

Mapping and variant calling pipeline

2021-06-30T14:18:19Z

Moiti001:

Author: Carolina Pita Barros 
Contact: carolina.pitabarros@wur.nl 
ABG 

[https://github.com/CarolinaPB/WUR_mapping-variant-calling Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads to a reference assembly. It outputs the mapped reads, a qualimap report and does variant calling.

==== Tools used: ====

* Bwa - mapping
* Samtools - processing
* Qualimap - mapping summary
* Freebayes - variant calling

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* ASSEMBLY - path to the assembly file
* PREFIX - prefix for the final mapped reads file

If you want the results to be written to this directory (not to a new directory), open Snakefile and comment out

<pre>workdir: config["OUTDIR"]</pre>
== RESULTS ==

* dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''sorted_reads''' directory with the file containing the mapped reads
* '''results''' directory containing the qualimap results
* '''variant_calling''' directory containing the variant calling VCF file

Mapping and variant calling pipeline

2021-06-30T14:17:06Z

Moiti001:

Author: Carolina Pita Barros
Contact: carolina.pitabarros@wur.nl
ABG

[https://github.com/CarolinaPB/WUR_mapping-variant-calling Link to the repository]

== First follow the instructions here: ==

[https://carolinapb.github.io/2021-06-23-how-to-run-my-pipelines/ Step by step guide on how to use my pipelines] 
Click [https://github.com/CarolinaPB/snakemake-template/blob/master/Short%20introduction%20to%20Snakemake.pdf here] for an introduction to Snakemake

== ABOUT ==

This is a pipeline to map short reads to a reference assembly. It outputs the mapped reads, a qualimap report and does variant calling.

==== Tools used: ====

* Bwa - mapping
* Samtools - processing
* Qualimap - mapping summary
* Freebayes - variant calling

=== Edit config.yaml with the paths to your files ===

<pre>OUTDIR: /path/to/output
READS_DIR: /path/to/reads/ # don't add the reads files, just the directory where they are
ASSEMBLY: /path/to/assembly
PREFIX: <output name></pre>
* OUTDIR - directory where snakemake will run and where the results will be written to
* READS_DIR - path to the directory that contains the reads
* ASSEMBLY - path to the assembly file
* PREFIX - prefix for the final mapped reads file

If you want the results to be written to this directory (not to a new directory), open Snakefile and comment out

<pre>workdir: config["OUTDIR"]</pre>
== RESULTS ==

* dated file with an overview of the files used to run the pipeline (for documentation purposes)
* '''sorted_reads''' directory with the file containing the mapped reads
* '''results''' directory containing the qualimap results
* '''variant_calling''' directory containing the variant calling VCF file