Control-FREEC improvements

CHANGELOG.md

@@ -10,6 +10,7 @@ Piellorieppe is one of the main massif in the Sarek National Park.
 
 ### Added - [2.6dev]
 
+- [#51](https://github.com/nf-core/sarek/pull/51) - Script for RankScore annotation
 - [#76](https://github.com/nf-core/sarek/pull/76) - Add `GATK Spark` possibilities to Sarek
 - [#87](https://github.com/nf-core/sarek/pull/87) - Add `GATK BaseRecalibrator` plot to `MultiQC` report
 - [#115](https://github.com/nf-core/sarek/pull/115) - Add [@szilvajuhos](https://github.com/szilvajuhos) abstract for ESHG2020

README.md

@@ -73,6 +73,7 @@ The nf-core/sarek pipeline comes with documentation about the pipeline, found in
     * [Complementary information about ASCAT](docs/ascat.md)
     * [Complementary information about Sentieon](docs/sentieon.md)
     * [Extra documentation on annotation](docs/annotation.md)
+    * [Downstream processing, ranking, and connection to Scout](docs/downstream.md)
 5. [Troubleshooting](https://nf-co.re/usage/troubleshooting)
 
 ## Credits

docs/downstream.md

@@ -0,0 +1,109 @@
+# Downstream processing, ranking and connection to [Scout](https://github.com/Clinical-Genomics/scout)
+
+## RankScore annotations
+
+[Scout](https://github.com/Clinical-Genomics/scout) was developed to help clinicians to visualise relevant
+variations, focusing on rare diseases. However, adding cancer capabilities are on the way, and the first step to add
+somatic mutations is to rank them in the VCF and load them to Scout. The scores have to be added as a `RankScore` INFO
+field, and have to have the format like:
+
+```text
+chr1  514206  rs1247069502  A C . PASS  ECNT=1;dbSNP=rs1247069502 RankScore=ScoreTest:0 GT:AD:AF  0/1:2,6:0.750 0/0:33,1:0.029
+```
+
+where the `ScoreTest` is the sample (normal) ID, and the score of this particular variation is 0. Scout will accept and
+show any sort of scores if the VCF is annotated properly. There have to be an info line in the VCF like:
+
+```text
+##INFO=<ID=RankScore,Number=.,Type=String,Description="The rank score for this variant ... . family_id:rank_score.">
+```
+
+RankScores for Scout considering rare diseases are calculated by [genmod](https://github.com/moonso/genmod), hence we
+are following the same annotation format to accomodate VCF files generated by Sarek.
+
+## Calculating RankScores with [Pathfindr](https://github.com/NBISweden/pathfindr)
+
+There are many approaches to score variants, Pathfindr considers many statistics found in the VCF file, collates
+[snpEff](http://snpeff.sourceforge.net/) and [VEP](https://www.ensembl.org/info/docs/tools/vep/index.html) annotations
+to calculate a single score that is written into a HTML page, also into a CSV file. This CSV can be processed with the
+`scripts/addPFRanksToVCF.py` python3 script. The CSV do not have to be generated by Pathfindr: until you have a
+comma-delimited file with a header with the `ID` in the first column and further columns named as `chr`,`start` and
+`rank_score`, the VCF will be annotated with the provided rank score in the CSV:
+
+```text
+ID,sample,rank_score,rank_terms,LOH,whatever,Consequence,IMPACT,SomethingIrrelevant,chr,start,end,width
+chr1:47248600_A/G,SampleTheBig,2,T1_gene,Y,0.926928,MODIFIER,0,0.853685,chr1,47248600,47248600,1
+chr1:110337540_C/T,SampleTheBig,2,T1_gene,Y,0.920036,MODIFIER,0,0.674548,chr1,110337540,110337540,1
+rs1205765240;rs4067337,SampleTheBig,2,T1_gene,"",0.605296,MODIFIER,0.348,0.209617,chr1,148963716,148963717,2
+chr1:197775910_G/A,SampleTheBig,2, TFBS not_canonical,"",0.813335,MODIFIER,0,0.471831,chr1,197775910,197775910,1
+rs142472845,SampleTheBig,2,T2_gene,"",0.519579,MODIFIER,0.002,0.140411,chr1,248458340,248458340,1
+chr2:15453359_C/T,SampleTheBig,2,T2_gene,"",0.514598,MODIFIER,0,0.00526232,chr2,15453359,15453359,1
+chr2:29413588_G/A,SampleTheBig,2,T1_gene,"",0.00738279,MODIFIER,0,0.994607,chr2,29413588,29413588,1
+chr2:70181162_C/T,SampleTheBig,2, high_impact,"",0.201275,HIGH,0,0.76965,chr2,70181162,70181162,1
+rs1335878992,SampleTheBig,2,T2_gene,"",0.00478691,MODIFIER,0,0.66847,chr2,80397367,80397367,1
+rs56269417,SampleTheBig,2,T2_gene,"",0.76779,MODIFIER,0,0.378918,chr2,98265733,98265733,1
+...
+```
+
+The basis of the annotation is `chr:startposition`: all variants in the VCF with the same coordinates will get the same
+RankScore irrespectively of the ALT allele. If there are multiple scores for a single genomic coordinate in the CSV, the
+very last one in the file will be stored only.
+
+## Usage of `scripts/addPFRanksToVCF.py` : transferring Pathfindr rank_score to Scout
+
+Usage is:
+
+```bash
+$ ./addPFRanksToVCF.py -h
+Usage: addPFRanksToVCF.py [OPTIONS]
+
+Options:
+  -v, --vcf TEXT     VCF file to annotate  [required]
+  -c, --csv TEXT     CSV file with ranked scores  [required]
+  -f, --family TEXT  family ID as expected in scout yaml file (usually sample
+                     name of the normal)  [required]
+  -h, --help         Show this message and exit.
+```
+
+Three parameters have to be provided, `-v` and `-c` are the name of the VCF and CSV files respectively, the `-f`
+parameter is the name of the normal sample in the VCF file. MuTect2 (GATK version < 4.0) VCF has a header line like:
+
+```text
+#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  TUMOR   NORMAL
+```
+
+not showing sample names, only the `TUMOR` and `NORMAL` columns. For Scout, the `NORMAL` part has to be altered to
+reflect actual sample name: this sample name will be also used in the RankScore annotation. By providing the something
+like `-f ThisIsIt` , the header line will be
+
+```text
+#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  TUMOR  ThisIsIt
+```
+
+and annotations will be like `RankScore=ThisIsIt:7` . The `.yaml` file for Scout to load the case has to be prepared in
+a way, that the `family` entry and the `sample_id` in the `samples:` sections also contain this key. Example:
+
+```text
+---
+
+owner:  BTB
+human_genome_build: 38
+# family has to be the same as the NORMAL column in the VCF
+# also have to be the same as the sample_id
+family: 'ScoreTest'
+# family name will appear in the Scout table
+family_name: 'Mutect2Testing'
+samples:
+  - analysis_type: wgs
+    sample_id: ScoreTest
+    tumor_type: medulloblastoma
+    phenotype: affected
+    father: 0
+    mother: 0
+    expected_coverage: 90
+    sex: male
+
+# This VCF contains the RankScore annotations already  
+vcf_cancer: testM2.vcf
+
+```

main.nf

@@ -109,6 +109,7 @@ def helpMessage() {
       --known_indels           [file] knownIndels file
       --known_indels_index     [file] knownIndels index
                                       If none provided, will be generated automatically if a knownIndels file is provided
+      --mappability            [file] Mappability file for Control-FREEC
       --species                 [str] Species for VEP
       --snpeff_db               [str] snpEff Database version
       --vep_cache_version       [int] VEP Cache version
@@ -125,6 +126,13 @@ def helpMessage() {
       --max_multiqc_email_size  [str] Theshold size for MultiQC report to be attached in notification email. If file generated by pipeline exceeds the threshold, it will not be attached (Default: 25MB)
       -name                     [str] Name for the pipeline run. If not specified, Nextflow will automatically generate a random mnemonic
 
+
+      --normal_pileup          [file] Normal pileup file for Control-FREEC
+      --tumor_pileup           [file] Tumor pileup file for Control-FREEC
+      --cf_coeff                [str] Control-FREEC coefficientOfVariation [0.015]
+      --cf_ploidy               [str] Control-FREEC ploidy [2]
+      --cf_window               [str] Control-FREEC window size [not set, as using coefficientOfVariation by default] 
+
     AWSBatch options:
       --awsqueue                [str] The AWSBatch JobQueue that needs to be set when running on AWSBatch
       --awsregion               [str] The AWS Region for your AWSBatch job to run on
@@ -361,6 +369,10 @@ if (!checkParameterList(annotateTools,annoList)) exit 1, 'Unknown tool(s) to ann
 
 // Check parameters
 if ((params.ascat_ploidy && !params.ascat_purity) || (!params.ascat_ploidy && params.ascat_purity)) exit 1, 'Please specify both --ascat_purity and --ascat_ploidy, or none of them'
+if ((params.normal_pileup && !params.tumor_pileup) || (!params.normal_pileup && params.tumor_pileup)) exit 1, 'Please specify both --normal_pileup and --tumor_pileup, or none of them'
+
+// for Control-FREEC window and coefficientOfVariation we need only one of them
+if (params.cf_window && params.cf_coeff) exit 1, 'Please specify either --cf_window OR --cf_coeff, but not both of them'
 
 // Has the run name been specified by the user?
 // This has the bonus effect of catching both -name and --name
@@ -477,6 +489,12 @@ ch_cadd_wg_snvs_tbi = params.cadd_wg_snvs_tbi ? Channel.value(file(params.cadd_w
 ch_pon = params.pon ? Channel.value(file(params.pon)) : "null"
 ch_target_bed = params.target_bed ? Channel.value(file(params.target_bed)) : "null"
 
+// parameters to experiment with Control-FREEC
+cf_window = params.cf_window ? Channel.value(file(params.cf_window)) : "null"
+cf_coeff = params.cf_coeff ? Channel.value(file(params.cf_coeff)) : "null"
+cf_ploidy = params.cf_ploidy ? Channel.value(file(params.cf_ploidy)) : "null"
+
+
 /*
 ================================================================================
                                 PRINTING SUMMARY
@@ -912,7 +930,7 @@ if (step in ['recalibrate', 'variantcalling', 'annotate']) {
 
 (inputBam, inputBamFastQC) = inputBam.into(2)
 
-// Removing inputFile2 wich is null in case of uBAM
+// Removing inputFile2 which is null in case of uBAM
 inputBamFastQC = inputBamFastQC.map {
     idPatient, idSample, idRun, inputFile1, inputFile2 ->
     [idPatient, idSample, idRun, inputFile1]
@@ -2194,8 +2212,7 @@ vcfFreebayesSingle = vcfFreebayesSingle.groupTuple(by: [0,1,2])
                              SOMATIC VARIANT CALLING
 ================================================================================
 */
-
-// Ascat, Control-FREEC
+// Ascat, pileup, pileups with no intervals, recalibrated, and IDs when the pileup is provided by CLI for Control-FREEC
 (bamAscat, bamMpileup, bamMpileupNoInt, bamRecalAll) = bamRecalAll.into(4)
 
 // separate BAM by status
@@ -2945,18 +2962,25 @@ process Mpileup {
         file(fastaFai) from ch_fai
 
     output:
-        set idPatient, idSample, file("${prefix}${idSample}.pileup.gz") into mpileupMerge
+        set idPatient, idSample, file("${prefix}${idSample}.pileup") into mpileupMerge
 
     when: 'controlfreec' in tools || 'mpileup' in tools
 
     script:
     prefix = params.no_intervals ? "" : "${intervalBed.baseName}_"
     intervalsOptions = params.no_intervals ? "" : "-l ${intervalBed}"
+
+    if(params.normal_pileup && params.tumor_pileup)   // we already have pileups
+    """
+    touch ${prefix}${idSample}.pileup 
+    """
+    else    
     """
+    # Control-FREEC reads uncompresses the zipped file TWICE in single-threaded mode.
+    # For the moment we are not using compressed pileups.
     samtools mpileup \
         -f ${fasta} ${bam} \
-        ${intervalsOptions} \
-    | bgzip --threads ${task.cpus} -c > ${prefix}${idSample}.pileup.gz
+        ${intervalsOptions} > ${prefix}${idSample}.pileup # | bgzip --threads ${task.cpus} -c > ${prefix}${idSample}.pileup.gz
     """
 }
 
@@ -2969,8 +2993,9 @@ if (!params.no_intervals) {
 }
 
 // STEP MPILEUP.2 - MERGE
-
 process MergeMpileup {
+    label 'cpus_1'
+
     tag {idSample}
 
     publishDir params.outdir, mode: params.publish_dir_mode, saveAs: { it == "${idSample}.pileup.gz" ? "VariantCalling/${idSample}/mpileup/${it}" : '' }
@@ -2979,19 +3004,25 @@ process MergeMpileup {
         set idPatient, idSample, file(mpileup) from mpileupMerge
 
     output:
-        set idPatient, idSample, file("${idSample}.pileup.gz") into mpileupOut
+        set idPatient, idSample, file("${idSample}.pileup") into mpileupOut
 
     when: !(params.no_intervals) && 'controlfreec' in tools || 'mpileup' in tools
 
     script:
+    if(params.normal_pileup && params.tumor_pileup)
+    """
+    if [[ ${params.normal_pileup} =~ $idSample ]] ; then ln -s ${params.normal_pileup} ${idSample}.pileup ; fi
+    if [[ ${params.tumor_pileup} =~ $idSample ]] ; then ln -s ${params.tumor_pileup} ${idSample}.pileup ; fi
     """
-    for i in `ls -1v *.pileup.gz`;
-        do zcat \$i >> ${idSample}.pileup
+    else
+    """
+    for i in `ls -1v *.pileup`;
+        do cat \$i >> ${idSample}.pileup
     done
 
-    bgzip --threads ${task.cpus} -c ${idSample}.pileup > ${idSample}.pileup.gz
+    #bgzip --threads ${task.cpus} -c ${idSample}.pileup > ${idSample}.pileup.gz
 
-    rm ${idSample}.pileup
+    #rm ${idSample}.pileup
     """
 }
 
@@ -3015,7 +3046,8 @@ mpileupOut = mpileupOut.map {
 // STEP CONTROLFREEC.1 - CONTROLFREEC
 
 process ControlFREEC {
-    label 'memory_singleCPU_2_task'
+    label 'cpus_max'
+    //label 'memory_singleCPU_2_task'
 
     tag {idSampleTumor + "_vs_" + idSampleNormal}
 
@@ -3031,28 +3063,44 @@ process ControlFREEC {
         file(fastaFai) from ch_fai
 
     output:
-        set idPatient, idSampleNormal, idSampleTumor, file("${idSampleTumor}.pileup.gz_CNVs"), file("${idSampleTumor}.pileup.gz_ratio.txt"), file("${idSampleTumor}.pileup.gz_normal_CNVs"), file("${idSampleTumor}.pileup.gz_normal_ratio.txt"), file("${idSampleTumor}.pileup.gz_BAF.txt"), file("${idSampleNormal}.pileup.gz_BAF.txt") into controlFreecViz
-        set file("*.pileup.gz*"), file("${idSampleTumor}_vs_${idSampleNormal}.config.txt") into controlFreecOut
+        set idPatient, idSampleNormal, idSampleTumor, file("${idSampleTumor}.pileup_CNVs"), file("${idSampleTumor}.pileup_ratio.txt"), file("${idSampleTumor}.pileup_normal_CNVs"), file("${idSampleTumor}.pileup_normal_ratio.txt"), file("${idSampleTumor}.pileup_BAF.txt"), file("${idSampleNormal}.pileup_BAF.txt") into controlFreecViz
+        set file("*.pileup*"), file("${idSampleTumor}_vs_${idSampleNormal}.config.txt") into controlFreecOut
 
     when: 'controlfreec' in tools
 
     script:
     config = "${idSampleTumor}_vs_${idSampleNormal}.config.txt"
     gender = genderMap[idPatient]
+    ploidy = params.cf_ploidy ? "${params.cf_ploidy}" : "2"
+    // if we are using coefficientOfVariation, we must delete the window parameter 
+    if(params.cf_coeff) {
+      coeff = "coefficientOfVariation = $params.cf_coeff"
+      window = ""
+    }
+    else if(params.cf_window) {
+      coeff = ""
+      window = "window = ${params.cf_window}"
+    }
+    else {  // default settings
+      coeff = "coefficientOfVariation = 0.015"
+      window = "" // it is "window = 20000" in the default settings, without coefficientOfVariation set, but we do not like it. Note, it is not written in stone
+    }
+    mappability = "gemMappabilityFile = ${params.mappability}"
     """
     touch ${config}
     echo "[general]" >> ${config}
     echo "BedGraphOutput = TRUE" >> ${config}
     echo "chrFiles = \${PWD}/${chrDir.fileName}" >> ${config}
     echo "chrLenFile = \${PWD}/${chrLength.fileName}" >> ${config}
-    echo "coefficientOfVariation = 0.05" >> ${config}
+    echo "${coeff}" >> ${config}
     echo "contaminationAdjustment = TRUE" >> ${config}
-    echo "forceGCcontentNormalization = 0" >> ${config}
+    echo "forceGCcontentNormalization = 1" >> ${config}
     echo "maxThreads = ${task.cpus}" >> ${config}
     echo "minimalSubclonePresence = 20" >> ${config}
-    echo "ploidy = 2,3,4" >> ${config}
+    echo "ploidy = ${ploidy}" >> ${config}
     echo "sex = ${gender}" >> ${config}
-    echo "window = 50000" >> ${config}
+    echo "${window}" >> ${config}
+    echo "${mappability}" >> ${config}
     echo "" >> ${config}
 
     echo "[control]" >> ${config}
@@ -3094,11 +3142,20 @@ process ControlFreecViz {
     when: 'controlfreec' in tools
 
     """
-    cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/assess_significance.R | R --slave --args ${cnvTumor} ${ratioTumor}
-    cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/assess_significance.R | R --slave --args ${cnvNormal} ${ratioNormal}
+    echo "Shaping CNV files to make sure we can assess signifacance"
+    awk 'NF==9{print}' ${cnvTumor} > TUMOR.CNVs
+    awk 'NF==7{print}' ${cnvNormal} > NORMAL.CNVs
+    echo "############### Calculating significance values for TUMOR CNVs #############"
+    cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/assess_significance.R | R --slave --args TUMOR.CNVs ${ratioTumor}
+    echo "############### Calculating significance values for NORMAL CNVs ############"
+    cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/assess_significance.R | R --slave --args NORMAL.CNVs ${ratioNormal}
+    echo "############### Creating graph for TUMOR ratios ###############"
     cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/makeGraph.R | R --slave --args 2 ${ratioTumor} ${bafTumor}
+    echo "############### Creating graph for NORMAL ratios ##############"
     cat /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/makeGraph.R | R --slave --args 2 ${ratioNormal} ${bafNormal}
+    echo "############### Creating BED files for TUMOR ##############"
     perl /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/freec2bed.pl -f ${ratioTumor} > ${idSampleTumor}.bed
+    echo "############### Creating BED files for NORMAL #############"
     perl /opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/freec2bed.pl -f ${ratioNormal} > ${idSampleNormal}.bed
     """
 }
@@ -3694,6 +3751,8 @@ def nfcoreHeader() {
     c_reset  = params.monochrome_logs ? '' : "\033[0m";
     c_white  = params.monochrome_logs ? '' : "\033[0;37m";
     c_yellow = params.monochrome_logs ? '' : "\033[0;33m";
+
+    return "--"
 
     return """    -${c_dim}--------------------------------------------------${c_reset}-
                                             ${c_green},--.${c_black}/${c_green},-.${c_reset}