Population structural variant calling pipeline - Revision history

Dawes0011 at 09:00, 16 June 2023

2023-06-16T09:00:29Z

← Older revision		Revision as of 11:00, 16 June 2023
Line 73:		Line 73:
	Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with		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>		<pre>/usr/bin/perl /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.		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.

Line 86:		Line 86:
	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.		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>		In the Snakefile, in rule <code>run_vep</code>, replace <code>/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>

	=== Installing R packages ===		=== Installing R packages ===
Line 99:		Line 99:
	If you get an error like this:		If you get an error like this:
	<pre>Warning in install.packages(new.packages) :		<pre>Warning in install.packages(new.packages) :
	'lib = "~~/cm~~/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>		'lib = "/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.		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.

Moiti001: added info about split read support option

2021-10-28T08:50:29Z

added info about split read support option

← Older revision		Revision as of 10:50, 28 October 2021
Line 24:		Line 24:
	* R - plot PCA		* R - plot PCA
	* SURVIVOR - basic SV stats		* SURVIVOR - basic SV stats
	* Python - add depth to vcf and create final table		* Python
	** ~~PyVcf~~		** PyVcf - add depth to vcf and create final table
			** bamgroupreads.py + samblaster - create bam files with split and discordant reads

	{\|		{\|
Line 42:		Line 43:
	SPECIES: <species_name>		SPECIES: <species_name>
	PREFIX: <output name>		PREFIX: <output name>
	NUM_CHRS: <number of chromosomes></pre>		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		* OUTDIR - directory where snakemake will run and where the results will be written to
	* READS_DIR - path to the directory that contains the reads		* READS_DIR - path to the directory that contains the reads
Line 58:		Line 61:
	* PREFIX - prefix for the created files		* PREFIX - prefix for the created files
	* NUM_CHRS - number of chromosomes for your species (necessary for plink). ex: 38		* 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		If you want the results to be written to this directory (not to a new directory), comment out or remove

Moiti001: added new config parameter - num of chrs

2021-10-15T11:49:16Z

added new config parameter - num of chrs

← Older revision		Revision as of 13:49, 15 October 2021
Line 4:		Line 4:

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


	== First follow the instructions here: ==		== First follow the instructions here: ==
Line 40:		Line 41:
	CONTIGS_IGNORE: /path/to/file		CONTIGS_IGNORE: /path/to/file
	SPECIES: <species_name>		SPECIES: <species_name>
	PREFIX: <output name></pre>		PREFIX: <output name>
			NUM_CHRS: <number of chromosomes></pre>
	* OUTDIR - directory where snakemake will run and where the results will be written to		* OUTDIR - directory where snakemake will run and where the results will be written to
	* READS_DIR - path to the directory that contains the reads		* READS_DIR - path to the directory that contains the reads
Line 55:		Line 57:
	* SPECIES - species name to be used for VEP		* SPECIES - species name to be used for VEP
	* PREFIX - prefix for the created files		* 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		If you want the results to be written to this directory (not to a new directory), comment out or remove

Moiti001 at 12:52, 8 October 2021

2021-10-08T12:52:38Z

← Older revision		Revision as of 14:52, 8 October 2021
Line 91:		Line 91:
	<pre>Warning in install.packages(new.packages) :		<pre>Warning in install.packages(new.packages) :
	'lib = "/cm/shared/apps/R/3.6.2/lib64/R/library"' is not writable</pre>		'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]		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 ==		== RESULTS ==
Line 99:		Line 99:
	** {prefix}.smoove-counts.html - shows a summary of the number of reads before and after filtering		** {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		* '''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 - Final VCF - with VEP annotation, not filtered for quality
	** {prefix}.smoove.square.vep.vcf.gz_summary.html - statistics from VEP		** {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}.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_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		** {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		** {prefix}_BND.vcf - vcf file with variants annotated with BND
	* '''6_metrics''' directory that contains general stats about the number of SVs found		* '''6_metrics''' directory that contains general stats about the number of SVs found
Line 111:		Line 111:

	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.		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.<br />
			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

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

2021-10-08T09:42:37Z

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

← Older revision		Revision as of 11:42, 8 October 2021
Line 44:		Line 44:
	* READS_DIR - path to the directory that contains the reads		* 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)		* 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<br />		* Example:
	~~>~~ sample2,sample2.bam,Pop1<br />		<blockquote>sample1,sample1.bam,Pop1<br />
	~~>~~ sample3,sample3.bam,Pop2<br />		sample2,sample2.bam,Pop1<br />
	~~>~~ sample4,sample4.bam,Pop2<br />		sample3,sample3.bam,Pop2<br />
			sample4,sample4.bam,Pop2<br /></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		* REFERENCE - path to the assembly file
	* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)		* CONTIGS_IGNORE - contigs to be excluded from SV calling (usually the small contigs)

Moiti001: added info about installing R packages

2021-10-08T09:20:57Z

added info about installing R packages

← Older revision		Revision as of 11:20, 8 October 2021
Line 59:		Line 59:
	<pre>OUTDIR: /path/to/outdir</pre>		<pre>OUTDIR: /path/to/outdir</pre>
	== ADDITIONAL SET UP ==		== 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:		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		Check if the cache file for your species already exist in <code>/lustre/nobackup/SHARED/cache/</code>. If it doesn’t, create it with
Line 77:		Line 77:

	In the Snakefile, in rule <code>run_vep</code>, replace <code>/cm/shared/apps/SHARED/ensembl-vep/vep</code> with <code>vep</code>		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 ==		== RESULTS ==

Moiti001: added info about population SV calling pipeline

2021-09-22T11:39:24Z

added info about population SV calling pipeline

New page

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.