Merge pull request #10 from TRON-Bioinformatics/bcftools-csq

Bcftools csq

Author: priesgo

Makefile

@@ -28,3 +28,4 @@ test:
 	bash test/scripts/run_test_10.sh
 	bash test/scripts/run_test_11.sh
 	bash test/scripts/run_test_12.sh
+	bash test/scripts/run_test_13.sh

README.md

@@ -36,7 +36,7 @@ Find the help as follows:
  ```
  $ nextflow run tron-bioinformatics/tronflow-vcf-postprocessing --help
  
- TronFlow VCF normalization v${VERSION}
+ TronFlow VCF postprocessing v${VERSION}
 
 Usage:
     nextflow run main.nf --input_vcfs input_vcfs --reference reference.fasta
@@ -67,7 +67,8 @@ Optional input:
     Example input file:
     sample1	primary:3
     sample1	metastasis:/path/to/metastasis.local_clonalities.bed
-    * --reference: path to the FASTA genome reference (indexes expected *.fai, *.dict) [required for normalization]
+    * --reference: absolute path to the FASTA genome reference (indexes expected *.fai, *.dict) [required for normalization and for functional annotation with BCFtools]
+    * --gff: absolute path to a GFF gene annotations file [required for functional annotation with BCFtools, only Ensembl-like GFF files]
     * --vcf-without-ad: indicate when the VCFs to normalize do not have the FORMAT/AD annotation
     * --output: the folder where to publish output
     * --skip_normalization: flag indicating to skip all normalization steps
@@ -269,7 +270,10 @@ No technical annotations are performed if the parameter `--input_bams` is not pa
 ## Functional annotations
 
 The functional annotations provide a biological context for every variant. Such as the overlapping genes or the effect 
-of the variant in a protein. These annotations are provided by SnpEff (Cingolani, 2012).
+of the variant in a protein. These annotations are provided by SnpEff (Cingolani, 2012) or by BCFtools csq (Danecek, 2017).
+Only one of the previous can be used.
+
+### Using SnpEff
 
 The SnpEff available human annotations are:
 * GRCh37.75 
@@ -289,7 +293,15 @@ To provide any additional SnpEff arguments use `--snpeff_args` such as
 `--snpeff_args "-noStats -no-downstream -no-upstream -no-intergenic -no-intron -onlyProtein -hgvs1LetterAa -noShiftHgvs"`, 
 otherwise defaults will be used.
 
-No functional annotations are performed if the parameters `--snpeff_organism` and `--snpeff_datadir` are not passed.
+No SnpEff functional annotations are performed if the parameters `--snpeff_organism` and `--snpeff_datadir` are not passed.
+
+### Using BCFtools csq
+
+BCFtools does not require any previous preparation. It expects two parameters: 
+* `--reference`: absolute path to the FASTA reference genome
+* `--gff`: absolute path to the Ensembl-like GFF annotations (ie: Gencode GFF files do not work https://github.com/samtools/bcftools/issues/1078)
+
+Importantly, BCFtools does use the available phasing information to evaluate all mutations affecting any given transcript together.
 
 
 ## References
@@ -298,3 +310,4 @@ No functional annotations are performed if the parameters `--snpeff_organism` an
 * Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, Li H. Twelve years of SAMtools and BCFtools. Gigascience. 2021 Feb 16;10(2):giab008. doi: 10.1093/gigascience/giab008. PMID: 33590861; PMCID: PMC7931819.
 * Di Tommaso, P., Chatzou, M., Floden, E. W., Barja, P. P., Palumbo, E., & Notredame, C. (2017). Nextflow enables reproducible computational workflows. Nature Biotechnology, 35(4), 316–319. 10.1038/nbt.3820
 * Cingolani P, Platts A, Wang le L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM. (2012) A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3.". Fly (Austin). 2012 Apr-Jun;6(2):80-92. PMID: 22728672
+* Danecek, P., & McCarthy, S. A. (2017). BCFtools/csq: haplotype-aware variant consequences. Bioinformatics, 33(13), 2037–2039. https://doi.org/10.1093/bioinformatics/btx100

main.nf

@@ -6,7 +6,7 @@ include { FILTER_VCF } from './modules/01_filter'
 include { BCFTOOLS_NORM; VT_DECOMPOSE_COMPLEX; REMOVE_DUPLICATES } from './modules/02_normalization'
 include { SUMMARY_VCF; SUMMARY_VCF as SUMMARY_VCF_2 } from './modules/03_summary'
 include { VAFATOR; MULTIALLELIC_FILTER } from './modules/04_vafator'
-include { VARIANT_ANNOTATION } from './modules/05_variant_annotation'
+include { VARIANT_ANNOTATION_SNPEFF; VARIANT_ANNOTATION_BCFTOOLS } from './modules/05_variant_annotation'
 
 params.help= false
 params.input_vcfs = false
@@ -19,6 +19,8 @@ params.input_clonalities = false
 
 
 params.reference = false
+params.gff = false
+
 params.output = "output"
 params.skip_normalization = false
 params.skip_decompose_complex = false
@@ -41,6 +43,10 @@ if ( params.snpeff_organism && ! params.snpeff_datadir) {
   exit 1, "To run snpEff, please, provide your snpEff data folder with --snpeff_datadir"
 }
 
+if (params.snpeff_organism && params.gff) {
+    exit 1, "Please use either SnpEff (--snpeff_organism) or BCFtools csq (--gff), but not both"
+}
+
 if (params.skip_normalization && ! params.input_bams && ! params.snpeff_organism) {
   exit -1, "Neither normalization, VAFator annotation or SnpEff annotation enabled! Nothing to do..."
 }
@@ -123,12 +129,6 @@ workflow {
             .join(input_purities.groupTuple(), remainder: true)
             .join(input_clonalities.groupTuple(), remainder: true)
 
-        //if (params.input_purities) {
-        //    vafator_input = vafator_input.join(input_purities.groupTuple())
-        //}
-        //if (params.input_clonalities) {
-        //    vafator_input = vafator_input.join(input_clonalities.groupTuple())
-        //}
         VAFATOR(vafator_input)
 
         final_vcfs = VAFATOR.out.annotated_vcf
@@ -139,8 +139,12 @@ workflow {
     }
 
     if (params.snpeff_organism) {
-        VARIANT_ANNOTATION(final_vcfs)
-        final_vcfs = VARIANT_ANNOTATION.out.annotated_vcf
+        VARIANT_ANNOTATION_SNPEFF(final_vcfs)
+        final_vcfs = VARIANT_ANNOTATION_SNPEFF.out.annotated_vcf
+    }
+    else if (params.gff) {
+        VARIANT_ANNOTATION_BCFTOOLS(final_vcfs)
+        final_vcfs = VARIANT_ANNOTATION_BCFTOOLS.out.annotated_vcf
     }
 }
 

modules/05_variant_annotation.nf

@@ -6,7 +6,7 @@ params.snpeff_organism = false
 params.snpeff_args = ""
 
 
-process VARIANT_ANNOTATION {
+process VARIANT_ANNOTATION_SNPEFF {
     cpus params.cpus
     memory params.memory
     publishDir "${params.output}/${name}", mode: "copy"
@@ -25,3 +25,25 @@ process VARIANT_ANNOTATION {
     snpEff eff ${datadir_arg} ${params.snpeff_args} -nodownload ${params.snpeff_organism} ${vcf} > ${name}.annotated.vcf
     """
 }
+
+process VARIANT_ANNOTATION_BCFTOOLS {
+    cpus params.cpus
+    memory params.memory
+    publishDir "${params.output}/${name}", mode: "copy"
+
+    conda (params.enable_conda ? "conda-forge::libgcc-ng=10.3.0 conda-forge::gsl=2.7 bioconda::bcftools=1.15.1" : null)
+
+    input:
+        tuple val(name), file(vcf)
+
+    output:
+    tuple val(name), file("${name}.annotated.vcf") , emit: annotated_vcf
+
+    """
+    bcftools csq \\
+        --fasta-ref ${params.reference} \\
+        --gff-annot ${params.gff} ${vcf} \\
+        --output-type v \\
+        --output ${name}.annotated.vcf
+    """
+}

nextflow.config

@@ -40,15 +40,15 @@ manifest {
 }
 
 params.help_message = """
-TronFlow variant normalization v${VERSION}
+TronFlow VCF postprocessing v${VERSION}
 
 Usage:
-    nextflow run main.nf -profile conda --input_vcfs input_vcfs --reference reference.fasta
+    nextflow run main.nf --input_vcfs input_vcfs --reference reference.fasta
 
 
 Input:
     * --input_vcf: the path to a single VCF to normalize (not compatible with --input_files)
-    * --input_vcfs: a tab-separated values file containing in each row the sample name and path to the VCF file (not compatible with --input_vcf)
+    * --input_vcfs: the path to a tab-separated values file containing in each row the sample name  and path to the VCF file (not compatible with --input_vcf)
     The input file does not have header!
     Example input file:
     sample1	/path/to/your/file.vcf
@@ -71,12 +71,13 @@ Optional input:
     Example input file:
     sample1	primary:3
     sample1	metastasis:/path/to/metastasis.local_clonalities.bed
-    * --reference: path to the FASTA genome reference (indexes expected *.fai, *.dict) [required for normalization]
+    * --reference: absolute path to the FASTA genome reference (indexes expected *.fai, *.dict) [required for normalization and for functional annotation with BCFtools]
+    * --gff: absolute path to a GFF gene annotations file [required for functional annotation with BCFtools, only Ensembl-like GFF files]
+    * --vcf-without-ad: indicate when the VCFs to normalize do not have the FORMAT/AD annotation
     * --output: the folder where to publish output
     * --skip_normalization: flag indicating to skip all normalization steps
     * --skip_decompose_complex: flag indicating not to split complex variants (ie: MNVs and combinations of SNVs and indels)
-    * --filter: specify a comma-separated list of filters to apply (e.g.: PASS,.), only variants with these values will be kept. If not provided all varianst are kept
-    * --vcf-without-ad: indicate when the VCFs to normalize do not have the FORMAT/AD annotation
+    * --filter: specify the filter to apply if any (e.g.: PASS), only variants with this value will be kept
     * --input_bams: a tab-separated values file containing in each row the sample name, tumor and normal BAM files for annotation with Vafator
     * --skip_multiallelic_filter: after VAFator annotations if any multiallelic variant is present (ie: two different
     mutations in the same position) only the highest VAF variant is kept unless this flag is passed

test/data/Homo_sapiens.GRCh37.87.minimal.gff3.gz

@@ -0,0 +1,16 @@
+#!/bin/bash
+
+### this test cannot be automated as it relies on SnpEff references which need to be downloaded beforehand
+source test/scripts/assert.sh
+output_folder=test/output/test13
+data_folder=`pwd`/test/data
+
+nextflow main.nf -profile test,conda --output $output_folder \
+  --gff ${data_folder}/Homo_sapiens.GRCh37.87.minimal.gff3.gz \
+  --reference ${data_folder}/ucsc.hg19.minimal.fasta
+
+test -s $output_folder/single_sample/single_sample.annotated.vcf || { echo "Missing test 1 output file!"; exit 1; }
+test -s $output_folder/tumor_normal/tumor_normal.annotated.vcf || { echo "Missing test 1 output file!"; exit 1; }
+
+assert_eq `grep -v '#' $output_folder/single_sample/single_sample.annotated.vcf| grep BCSQ | wc -l | cut -d' ' -f 1` 32 "Wrong number of variants"
+assert_eq `grep -v '#' $output_folder/tumor_normal/tumor_normal.annotated.vcf | grep BCSQ | wc -l | cut -d' ' -f 1` 32 "Wrong number of variants"

test/scripts/run_test_13.sh

@@ -2,7 +2,7 @@
 
 ### this test cannot be automated as it relies on SnpEff references which need to be downloaded beforehand
 source test/scripts/assert.sh
-output_folder=test/output/test8
+output_folder=test/output/test9
 snpeff_datadir=/home/you/snpeff
 
 nextflow main.nf -profile test,conda --output $output_folder --snpeff_organism hg19 --snpeff_datadir $snpeff_datadir