iarcbioinfo/needlestack

Multi-sample somatic variant caller

A multi-sample somatic variant caller

A preprint describing needlestack is available on biorxiv: https://www.biorxiv.org/content/10.1101/639377v1

Contact: [email protected]

Jump over here for a quick start!

Needlestack development is support by the US National Cancer Institute (grant number R21CA175979) and the French Institut national du cancer.

WARNING

Since the last release, we have harmonized our pipelines and the following option names changed:

• bam_folder is now input_bams and accepts text file
• fasta_ref is now ref
• max_DP is now max_dp
• out_folder is now output_folder
• out_vcf is now output_vcf and is now a mandatory argument
• out_annotated_vcf is now output_annotated_vcf
• pairs_file is now tn_pairs
• no_plots is now plots (see the Detailed description section)

Description

Needlestack is an ultra-sensitive multi-sample variant caller for Next Generation Sequencing (NGS) data. It is based on the idea that analysing several samples together can help estimate the distribution of sequencing errors to accurately identify variants. It has been initially developed for somatic variant calling using very deep NGS data from circulating free DNA, but is also applicable to lower coverage data like Whole Exome Sequencing (WES) or even Whole Genome Sequencing (WGS). It is a highly scalable and reproducible pipeline thanks to the use of nextflow and Docker/Singularity container technologies.

Here is a summary of the method:

• At each position and for each candidate variant, we model sequencing errors using a negative binomial regression with a linear link and a zero intercept. The data is extracted from the BAM files using samtools.
• Genetic variants are detected as being outliers from the error model. To avoid these outliers biasing the regression we use a robust estimator for the negative binomial regression (published here with code available here).
• We calculate for each sample a p-value for being a variant (outlier from the regression) that we further transform into q-values to account for multiple testing.

Note that needlestack is not performing a re-assembly of reads to determine haplotypes as some other modern variant callers are doing (Mutect2, Strelka2 etc.). For this reason we recommend using ABRA2, a realigner for NGS data, before running needlestack. We also host a nextflow pipeline for ABRA2.

The detection of a variant can be visualized graphically. Here are two examples of a variant detected, one from a series of samples with a coverage of 120X (left) and one with a very high coverage (50000X, right). Each dot represents a sequenced sample colored according to its Phred scale q-value. The black regression line shows the estimated sequencing-error rate along with the 99% confidence interval (black dotted lines) containing samples. Colored-dotted lines correspond to the limits of regions defined for different significance q-value thresholds. In both case one sample appears as outlier from the regression (in red), and is therefore classified as carrying the given mutation. Note that the position on the left has a high error rate, while the position on the right has a very low error rate, allowing the detection of a variant with a tiny allelic fraction (0.025%).

Dependencies

Needlestack works under most Linux distributions and Apple OS X.

1. This pipeline is based on nextflow. As we host several nextflow pipelines, we have centralized the common information in the IARC-nf repository. Please read it carefully as it contains essential information for the installation, basic usage and configuration of nextflow and our pipelines.

2. External software:

• bedtools
• samtools
• R with Rscript
• Compile the file mpileup2readcounts.cc located here
• Add the previous tools to your path (executables are assumed to be respectively called samtools, Rscript and mpileup2readcounts)
• Optionally, the bioconductor packages required for the alignments plots (see the Detailed description section)

You can avoid installing all the external software by using Docker or Singularity. See the IARC-nf repository for more information. Warning: the docker/singularity containers don't include the bioconductor packages needed to plot alignments.

To use needlestack without nextflow, in addition to the previous tools, download the files in this bin directory and add them to your path. See the Usage section to run needlestack without nextflow.

Input

Type	Description
BAM files	BAM files (called *.bam) grouped in a single folder or described in a text file (one BAM path per line) along with their index files (called *.bam.bai or *.bai). A minimum of 20 BAM files is recommended.
fasta file	A reference fasta file (eventually compressed with bgzip) along with its faidx index (and *.gzi faidx index if compressed).
bed file	Optional input file. Otherwise the variant calling is performed on the whole reference provided.

A sample dataset is available here for testing.

Parameters

Type --help to get the full list of options. All parameters are prefixed with a double-dash like in --help.

• Mandatory

Name	Example value	Description
input_bams	BAM/ or bams.txt	Folder containing the BAM files
ref	fasta_name.fa	Reference fasta file
output_vcf	vcf_name.vcf	Name of the output VCF file generated by needlestack

• Optional

Name	Default value	Description
min_dp *	30	Minimum median coverage to consider a site. In addition, at least 10 samples have to be covered by min_dp.
min_ao *	3	Minimum number of non-ref reads in at least one sample to consider a site
min_af *	0	Minimum allelic fraction in at least one sample to consider a site
nsplit	1	Split the bed file in nsplit pieces and run in parallel
min_qval	50	qvalue threshold in Phred scale to consider a variant
sb_type	SOR	Strand bias measure, either SOR, RVSB or FS
sb_snv	100 or 1000	Strand bias threshold for SNVs. The default of 100 (1000 if sb_type is FS) means no filter
sb_indel	100 or 1000	Strand bias threshold for indels. The default of 100 (1000 if sb_type is FS) means no filter
map_qual	0	Min mapping quality (passed to samtools)
base_qual	13	Min base quality (passed to samtools)
max_dp	50000	Downsample coverage per sample (passed to samtools)
plots	SOMATIC if tn_pairs provided, ALL if not	Creates pdf plots of regressions in the output. To remove plots set --plots to NONE (See Plot options paragraph).
genome_release	-	Reference genome for the annotions on the alignments plots. Examples: Hsapiens.UCSC.hg19, Hsapiens.UCSC.hg18, Hsapiens.UCSC.hg38, Mmusculus.UCSC.mm10. The right terminology is described below (Plot options paragraph). WARNING: this option is mandatory if the --do_alignments flag is used
output_folder	--bam_folder	Output folder, by default equals to the input bam folder
bed	-	BED file containing a list of regions (or positions) where needlestack will perform variant calling
region	-	A region in format CHR:START-END where calling should be done
tn_pairs	-	A tab-delimited file containing two columns (normal and tumor sample names) for each sample in line. This enables matched tumor/normal pair calling features (see below)
sigma_normal	0.1	Sigma parameter for negative binomial modeling of expected germline allelic fraction. We strongly recommend not to change this parameter unless you really know what it means
min_af_extra_rob	0.2	Minimum allelic fraction to exclude a sample at a position for extra-robust regression
min_prop_extra_rob	0.1	Minimum proportion of samples having an allelic fraction to be excluded from extra-robust regression
max_prop_extra_rob	0.8	Maximum proportion of samples having an allelic fraction to be excluded from extra-robust regression
power_min_af		Allelic fraction used to classify genotypes as 0/0 or ./. depending of the power to detect a variant at this fraction (see below)
input_vcf		A VCF file (from GATK) where calling should be done. Needlestack will extract DP and AO from this VCF (DP and AD fields) and annotate it with phred q-value score (FORMAT/QVAL field), error rate (INFO/ERR) and over-dispersion sigma parameter (INFO/SIG). WARNING: by default, only works with split (coreutils) version > 8.13

* Caution: min_dp, min_ao and min_af are position-based filters, which means that needlestack might still identify a variant with an allelic fraction lower than min_af in a sample if at the same position, another sample has an allelic fraction larger than min_af.

By default, if neither --bed nor --region are provided, needlestack runs on the whole reference genome provided, building a bed file from fasta index. If --bed and --region are both provided, it runs on the region only.

• Flags

Flags are parameters without value.

Name	Description
help	Display help
all_SNVs	Output all SNVs, even when no variant is detected
extra_robust_gl	Perform extra-robust regression (useful for common germline SNPs, see below)
no_labels	No label for the outliers on regression plots
no_indels	Do not perform variant calling for insertions and deletions
no_contours	Do not plot q-values contours (for q-value threshold={10,30,50,70,100} by default) and do not plot minimum detectable allelic fraction as a function function of coverage
use_file_name	Use the bam file names as sample names. By default the sample name is extracted from the bam file SM tag
do_alignments	Add alignments plots to the pdf plots of regressions

Usage

1. Optionally download a sample dataset.

   git clone --depth=1 https://github.com/IARCbioinfo/data_test

2. Run the pipeline using nextflow (here using docker).

   cd data_test
   nextflow run iarcbioinfo/needlestack -with-docker \
                    --bed BED/TP53_all.bed \
                    --input_bams BAM/BAM_multiple/ \
                    --ref REF/17.fasta \
                    --output_vcf all_variants.vcf 

   You will find a VCF file called all_variants.vcf in the BAM_multiple/ folder once done.

   The first time it will take more time as the pipeline will be downloaded from github and the docker container from dockerhub.

   Official releases can be found here. There is a corresponding official docker container for each release and one can run a particular version using (for example for v0.3):

   nextflow run iarcbioinfo/needlestack -r v0.3 -with-docker \
   	        --bed BED/TP53_all.bed --bam_folder BAM/BAM_multiple/ --fasta_ref REF/17.fasta --out_vcf all_variants.vcf

It is also possible to run the pipeline without nextflow (not recommended) using the needlestack.sh. Here on the example dataset downloaded above: bash cd data_test/ needlestack.sh --region=17:7572814-7573814 --bam_folder=BAM/BAM_multiple --ref=REF/17.fasta --output_vcf=all_variants.vcf

Output

Type	Description
VCF file	File containing all the variants detected. It contains a header of meta-information and data lines containing information about each variant position
PDF files	For each variant detected a pdf is generated, this pdf file contains the regression plots and alignments plots if the --do_alignments flag is used

Detailed description

Germline, somatic, matched Tumor-Normal pairs calling and contamination

When using matched tumor/normal, Needlestack can classify variants (VCF FORMAT/STATUS field) according to the following table:

For this, one need to provide a tab-delimited file containing two columns with normal and tumor sample names using the --tn_pairs option. The first line of this file is a header with TUMOR and NORMAL keywords. When one normal or one tumor is missing, one can write NA. In this mode, the parameter power_min_af defines the allelic fraction in the tumor one is trying to detect to classify genotypes as ./. or 0/0 depending on the power to detect this allelic fraction. Variants found as somatic in a tumor, but germline in another sample of the series will be flagged as POSSIBLE_CONTAMINATION. We found this particularly important, as needlestack is very sensitive to low allelic fractions, to filter out contamination among samples, coming for example from pooled exome capture.

In other cases (when there is no --tn_pairs parameter defined), genotypes are defined as ./. or 0/0 assuming one is looking for allelic fractions expected for germline variants (negative binomial distribution centered at 0.5 with over-dispersion parameter sigma=sigma_normal, with sigma_normal=0.1 by default). If you are looking for somatic variants without matched-normal and assuming you are interesting to correctly distinguish ./. and 0/0genotypes, you can set the power_min_af parameter to the lowest allelic fraction of somatic variants you are interested with (and your coverage allows you to find). Note that this is by far not the most common situation, and that in most cases you don't have to worry about the power_min_af parameter.

Plot options

Two options to consider:

1. --plots:

   • SOMATIC: produce pdf regression plots (like in the description above) only for somatic variants (option not available in the annotation mode). Default value when using matched tumor/normal.
   • ALL: produce pdf regression plots for all variants. Default value when not using matched tumor/normal or when using annotation mode.
   • NONE: remove pdf regression plots from the output.

2. --do_alignments: To add the alignments plots after the regression plots. If this flag is used, the name of the reference genome (--genome_release option) needs to be provided to choose the correct annotation (See --genome_release option below).

The following bioconductor packages need to be installed for plotting the alignments:

• The Gviz package
• An Annotation package for TxDb objects.
• A Genome wide annotation, it contains mappings between Entrez Gene identifiers and GenBank accession numbers. Examples: the Genome wide annotation package for Human: org.Hs.eg.db and for the mouse: org.Mm.eg.db.

Example: If the hg19 release of the human genome is used, the following packages should be installed: Gviz, TxDb.Hsapiens.UCSC.hg19.knownGene (hg18 and hg38 UCSC version can also be used) and org.Hs.eg.db

These packages exist for other organisms than Human but have not been tested. One can for example generate the alignments plot for mouse data by installing TxDb.Mmusculus.UCSC.mm10.knownGene (mm9 UCSC version can also be used) and org.Mm.eg.db. For the other organisms the packages need to have the same nomenclature as the ones listed above.

The --genome_release option needs to be provided and corresponds to the TxDb annotation package name without its prefix and suffix. For the hg19 release of the human genome, one needs to set --genome_release to Hsapiens.UCSC.hg19. Note that the packages chosen for the annotations are compatible with the UCSC notations since most of the Gviz fonctionalities can handle these notations. The reference genome used for the BAM alignments can be based on GENCODE, UCSC or ENSEMBL genome varieties.

Warning: the docker/singularity containers don't include the bioconductor packages needed to plot alignments.

The alignment plot represents from top to bottom:

• Chromosome representation: a red vertical line shows the position of the variant.
• Genomic axis associated with the alignment.
• BAM alignment(s) (50 bases on both sides of the variant): the variant position is highlighted in red. Note that insertions and deletions are not currently supported by Gviz.
• The reference genome: the variant position is highlighted in red.
• Genome annotation: the yellow blocks represent exons, the variant position is highlighted in red. The annotation is represented only if the variant is not in an intergenic region.
• Zoom-out on the genome annotation: representation of the whole gene with its name on the left side. The annotation is represented only if the variant is not in an intergenic region. A red vertical line shows the position of the variant.
• Genomic axis corresponding to the previous gene annotation. A red vertical line shows the position of the variant.

Common variants

Needlestack is designed to identify rare variants (i.e. only a few samples in your set of samples have a particular variant), because of the robust regression used. Therefore, common SNPs (>10%) or strong somatic hotspots will be missed. The optional extra_robust_gl can overcome partially this issue for common germline mutations: it first discard high allelic fraction (>20%, assuming these are likely true variants) before fitting the regression model when between 10% and 50% of sample have such a high allelic fraction. A flag is written in the VCF INFO/WARN field when this happened (EXTRA_ROBUST_GL). Additionally, when an other allele than the reference allele is the most common, it is taken as the reference allele and a flag is also written in the VCF (INFO/WARN=INV_REF).

Strand bias

For conventional variant callers, GATK WES/WGS recommended values for SOR strand bias are SOR < 4 for SNVs and < 10 for indels. We haven't found this particularly useful, but a more detailed evaluation is necessary. For amplicon based targeted sequencing, RVSB>0.85 seems to reduce some errors. There is no hard filter by default as this is easy to do afterward using bcftools filter command.

Contributions

Name	Email	Description
Matthieu Foll*	[email protected]	Lead developer, contact for support
Tiffany Delhomme*	[email protected]	Main developer, contact for support
Nicolas Alcala	[email protected]	tumor/normal pair mode, power calculation
Aurelie Gabriel	[email protected]	Alignment plot, nextflow pipeline, mpileup2readcounts, bash version

References

• Needlestack: an ultra-sensitive variant caller for multi-sample next generation sequencing data Tiffany M. Delhomme, Patrice H. Avogbe, Aurelie A.G. Gabriel, Nicolas Alcala, Noemie Leblay, Catherine Voegele, Maxime Vallee, Priscilia Chopard, Amelie Chabrier, Behnoush Abedi-Ardekani, Valerie Gaborieau, Ivana Holcatova, Vladimir Janout, Lenka Foretova, Sasa Milosavljevic, David Zaridze, Anush Mukeriya, Elisabeth Brambilla, Paul Brennan, Ghislaine Scelo, Lynnette Fernandez-Cuesta, Graham Byrnes, Florence Le Calvez-Kelm, James D D. McKay, Matthieu Foll bioRxiv 639377; doi: https://doi.org/10.1101/639377. https://www.biorxiv.org/content/10.1101/639377v1

Needlestack has been used in these papers:

• Urinary TERT promoter mutations as non-invasive biomarkers for the comprehensive detection of urothelial cancer Avogbe et al. 2019 EBioMedicine.
• Fernandez-Cuesta L, Perdomo S, Avogbe PH, Leblay N, Delhomme TM, Gaborieau V, Abedi-Ardekani B, Chanudet E, Olivier M, Zaridze D, Mukeria A, Vilensky M, Holcatova I, Polesel J, Simonato L, Canova C, Lagiou P, Brambilla C, Brambilla E, Byrnes G, Scelo G, Le Calvez-Kelm F, Foll M, McKay JD, Brennan P. Identification of Circulating Tumor DNA for the Early Detection of Small-cell Lung Cancer. EBioMedicine. 2016 Aug;10:117-23. doi: 10.1016/j.ebiom.2016.06.032. Epub 2016 Jun 25. PubMed PMID: 27377626; PubMed Central PMCID: PMC5036515.
• Le Calvez-Kelm F, Foll M, Wozniak MB, Delhomme TM, Durand G, Chopard P, Pertesi M, Fabianova E, Adamcakova Z, Holcatova I, Foretova L, Janout V, Vallee MP, Rinaldi S, Brennan P, McKay JD, Byrnes GB, Scelo G. KRAS mutations in blood circulating cell-free DNA: a pancreatic cancer case-control. Oncotarget. 2016 Nov 29;7(48):78827-78840. doi: 10.18632/oncotarget.12386. PubMed PMID: 27705932; PubMed Central PMCID: PMC5346680.

needlestack's Issues

Let min_qval be equal to 0

User must have the choice to keep variants with min_qval=0 (1 in non-Phred), e.g if want to keep all positions.
That is the nextflow script which has to be modify.

Add an option to directly input a region for calling in the command line

Format like chr1:123456-123457. In this case bed file can be omitted and is ignored if present.

error when coverage is null for every bam file

If we try to call a region where the coverage is null for every files, samtools_mpileup will produce an empty file and mpileup2table will not output any sample*.txt whereas they are required outputs.

Error running process > Missing output file(s): 'sample*.txt' expected by process: mpileup2table

We have to stop the pipeline after samtools_mpileup process if all emitted mpileup files are empty.

Add information about the pipeline in the log

See here for example.

Change the way we publish new version

We should keep the master branch fixed between version releases. Then we have a branch for the next version where we make pull requests from other(s) dev branches. Once the new version is ready we merge it to master and delete the "next version" branch. We can then recreate a "next version" branch and so on. This is linked to #68

Change the number of entry in the INFO and FORMAT VCF fields

Number=A means one value per alternate allele. It is more general but in our case we always have only one alternate allele per line, so we can safely put 1 instead. This makes easier to process the VCF file in bioconductor for example.

QVAL is wrongly called GQ for indels

Check if the gzi is present if the ref is gz

Improve R script command line parsing

Use @tdelhomme code example here https://github.com/IARC-bioinfo/R-tricks#building-an-argument-section-for-your-r-script

Move to IARC-bioinfo organisation repo

Check that BAM folder contains bam files

Send an error if not (check bai too) and a warning if less than 20? (refuse to run with less than 10?)

Calling doesn't work when a region contains only T in the reference

This is because when R reads the pileup files here Tin the reference are converted automatically to booleans TRUE by read.table.

samtools mpileup is fixed and we no longer need to add `*` in zero coverage regions

We can remove the | sed 's/ / * */g' from here https://github.com/IARCbioinfo/needlestack/blob/v0.3b/needlestack.nf#L236 since samtools mpileup bug has been fixed:
samtools/samtools@afae50c

It's not yet present in samtools 1.3 so we should wait until the next release to do it.

This will solved an unexpected behaviour in needlestack: because of the pipe, even if samtools returns a non-zero exit code, the command returns a zero exit code (see http://unix.stackexchange.com/questions/14270/get-exit-status-of-process-thats-piped-to-another). One option would be to use pipefail in bash but it's not compatible with all bash versions and also doesn't work will with head -n 1like here https://github.com/IARCbioinfo/needlestack/blob/v0.3b/needlestack.nf#L276. For example:

$ seq 1 100000 | head -1 ; echo $?
1
0
$ set -o pipefail ; seq 1 100000 | head -1 ; echo $?
1
141

The regression is not always robust for common SNPs

Look at this example:

I think we can safely exclude AF>0.25 before fitting the regression to avoid this?

In the absence of a bed file the pipeline should run on the full reference genome

bedtools has a function makewindows that can create a bed from a fasta index file with windows of a given size. For example split the whole genome in 10Mb regions:

bedtools makewindows -g reference.fasta.fai -w 10000000

I could actually also use this function to split the bed instead/in combination with my own R script.

Make a script that would run needlestack without nextflow

It should be easy to simply take the shell portions from nextflow script. Of course we just ignore the first one splitting the bed file. The difficult part could be to parse arguments (and to give them the same names) and to make sure we keep it up-to-date when we update the nextflow script. For arguments parsing, have a look at https://gist.github.com/jehiah/855086

Add version numbers in VCF output

It would be good to write at least the pipeline version, nextflow version and the docker container version.

Add program usage when launched with --help

See here for example

Add logo image

VCF files have to be sorted

See here:

All entries for a specific CHROM should form a contiguous block within the VCF file.
Positions are sorted numerically, in increasing order, within each reference sequence CHROM.

Use something like: https://github.com/vcftools/vcftools/blob/master/src/perl/vcf-sort or https://www.biostars.org/p/133487/

Add alignment plot in the PDF

Look at ggbio or Gviz R packages.

Add an option to use regression parameters previously estimated for SNVs

This could be used to call variants from a single sample if we already estimated the parameters from other data (hoping the batch effect won't be too important). The input format could be a VCF obtained from needlestack with all_SNVs option as it contains the two parameters of the regression.

Choose strand bias filter to apply

It would be better to be able to choose between SOR and RVSB

color points by qvalues in regression plot

for example take colors on a rainbow, don't forget to show the corresponding palette (same as for benchmarking? i.e. stopped at 100)

par(mfrow=c(1,2))
vect=sample(1:10,10)
plot(vect,vect,col=rainbow(length(vect)),pch=19)
plot(vect,vect,col=rainbow(length(vect))[match(vect,sort(vect))],pch=19)

Remove unnecessary intermediate outputs

The only important outputs are: the final VCF and PDF files. All other intermediate files should only be present in nextflow work directory. This will imply:

• removing storeDir directives in split_bed, samtools_mpileup and mpileup2table processes
• deleting the intermediate VCF outputs at the end of collect_vcf_result (see how we delete the empty pdf files).
• renaming the VCF output dir as PDF
• deleting the PDF folder and the all_variants.vcf file as the first step in the pipeline

It will has the advantage that the pipeline will be restarted from the beginning by default. We should note in the manual that the output is erased when it exists already. Users will still have the option to resume an interrupted pipeline by adding the -resume option to the command line (add that to the manual too).

Sometimes large number in VCF files are written in scientific notations

This can be a problem in particular for round chromosome positions.

Add additional filters on our calls

Look what Mutect and GATK VariantAnnotator do.

Correct english typos in readme, help and log

Ask a native speaker

Choose a better name for the pipeline and change file names accordingly

Don't forget the VCF header line

CircleCI docker build image sometimes crashes

Error looks like:

E: Failed to fetch http://httpredir.debian.org/debian/pool/main/d/dbus-python/python3-dbus_1.2.0-2+b3_amd64.deb  Error reading from server. Remote end closed connection

It seems to be a common issue with docker: http://serverfault.com/questions/690639/api-get-error-reading-from-server-under-docker

Add an option to choose output VCF file name (--out_vcf?)

Add optional additional informations into VCF contig header

Usual fields are:

• assembly
• md5 (could be calculated)
• URL

Add contour lines for a set of qvalues in the plot

The idea is to compute, at each position, for a grid (DP,AO) of n points, the corresponding qvalues if we add the point in the observations.
For each point of the grid, its qvalue will correspond to a transformed p-value computed with pnbinom using estimated parameters of the regression including the point.
This will allow us to plot contour lines for a set of qvalues.
From this, we could also compute, for a specific q-value threshold, the distribution of min_af (=AO/DP in the contour line of the threshold) in function of DP.

Change the way chromosome length is calculated

I just discovered that the size of chromosomes/contigs is written in the faidx index of the fasta reference file (see HTSlib). So my script to calculate that is useless (and very slow for nothing). I should rather simply use the faidx.

Add on option to give a VCF as an input to only recall variants it contains

In this case the bed file is not needed. The problem will be non-standard notation of variants: they should be left aligned/normalized first with only one ALT allele per line.

Add more tests in CircleCI

We should at least test each option once. Maybe we don't need to do all tests twice with the stable docker file (just try it in one case).

Verify the user inputs are correct

For parameters:

• sample_names
• do_plots
• all_sites

See here for example

Dockerfile always adds scripts from master branch

Docker image should only contain binaries and not scripts, as nextflow is adding the scripts from the bin directory to the environment in the docker container. So adding scripts in docker image is useless.

Option `all_sites` should rather be called `all_SNVs`

Make the stable docker file more stable

Stable docker file should become the default and the one installing the latest version should be on a dev/beta branch only. We should only update the docker file at each version update. This is linked to #69

Add a zoomed regression plot

The y-axis max value should be twice the max of the 99% confidence line. Also make a zoomed q-value / AF plot with the same corresponding max y-axis. In this example we don't really see how the regression fits the errors:

Pipeline crashed when no variant is found

It seems that when the process R_regression does not produce variants, as it does not produce pdf output (but the errorStrategy 'ignore' makes that acceptable), the empty vcf is not sent to the vcf channel:

output:
     file "${region_tag}.vcf" into vcf
     file '*.pdf' into PDF

This creates an error in the collect_vcf_result process as there is no vcf there.

Add contigs in VCF header

bcftools view wants it to subset samples (see samtools/bcftools#326)

Improve the bed split method

My R script is not splitting regions in the bed file but is rather grouping them to launch in parallel. This is good for exomes as the bed file contains many line (~ one per exon), but the user should also have the option to allow cutting lines of the bed to make calculation faster when they contain only a few large lines. This is related to #39 where for full genomes there is no bed so we will have to cut the chromosomes in chunks.

CircleCI deploy.sh doesn't trigger correctly Docker Hub

add --no_indel option

Let the user choose if calling indels is not necessary.
Add this option in the R script.
To improve accuracy, modifiy the building of pileup file with the perl script.

Cutting bed files with zero length regions was not working properly

plot improved error rate confidence interval

Plot a better estimated confidence interval (0.01, 0.99) according to the data.
At this time, we compute only the confidence interval for the maximum coverage with qnbinom(), and plot a line from 0 to this value (see big blue dotted lines on the following plot).
The issue is that the shape of this interval is not necessary linear.
Now we would like to compute this confidence interval for a sufficient number of points to be sufficiently accurate. Let's take 1000 equidistant points in {0,max(coverage)} if max(coverage)>1000, otherwise all coverage between 0 and the maximum coverage.

Use this opportunity to change colors for error rate and confidence interval from blue to black to erase conflict with the contour for qvalue=10.

Add strand bias estimated using Fisher's Exact Test

Like in GATK

The pipeline randomly crashes with java.nio.file.NoSuchFileException: XXX_empty.pdf

The issue seems to be related to an NFS copy error on our IARC cluster. Sometimes an empty pdf file is not correctly copied at the end of R_regression process, but as nextflow is tracking it, it crashes later on when it tries to delete it.

It's not a bug in our pipeline neither in nextflow.

Nextflow is not checking at the end of a process if files are correctly copied. The last line of the R_regression file .command.run.1 is:

cp -fR *.pdf /path/to/my/folder/BAM/VCF || true

I do see the XXX_empty.pdf in the corresponding nextflow workfolder but it has not been copied into the /path/to/my/folder/BAM/VCF folder, even though the process ended with no error.